text:
using RoME
using Statistics
using Manifolds
# , Distributions
using Test
using DistributedFactorGraphs
using TensorCast
import Base: convert

##

@testset "test sampling from BearingRange factor..." begin

##

p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))

fg = initfg()
addVariable!(fg, :x0, Pose2)
addVariable!(fg, :x1, Point2)
addFactor!(fg, [:x0;:x1], p2br, graphinit=false)

meas = sampleFactor(IIF._getCCW(fg, :x0x1f1), 100)
##

# meas = getSample(p2br, 100)
M = getManifold(p2br)
mcords = vee.(Ref(M), Ref(identity_element(M)), meas)

mu = Statistics.mean(mcords)
sigma = Statistics.std(mcords)

@test abs(mu[1]) < 0.1
@test 0.05 < abs(sigma[1]) < 0.2

@test abs(mu[2] - 20.0) < 1.0
@test 0.5 < abs(sigma[2]) < 1.5

##

end

@testset "test BearingRange factor residual function..." begin

##

# dummy variables
fg = initfg()
X0 = addVariable!(fg, :x0, Pose2)
X1 = addVariable!(fg, :x1, Point2)

##

p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))

xi = getPointIdentity(Pose2)
li = zeros(2); li[1] = 20.0;

M = getManifold(p2br)
_zi = [0,20.0]
zi = Manifolds.hat(M, identity_element(M), _zi)

res = calcFactorResidualTemporary( p2br, (Pose2, Point2), zi, (xi, li))
#
# calcFactorResidualTemporary(p2br, zi, (Pose2, xi), (Point2, li))

@show res
@test norm(res) < 1e-14

##

xi = getPointIdentity(Pose2)
li = zeros(2); li[2] = 20.0;
_zi = [pi/2,20.0]
zi = Manifolds.hat(M, identity_element(M), _zi)
# idx = 1
# res = zeros(2)
# p2br(res, fmd, idx, zi, xi, li)

res = calcFactorResidualTemporary( p2br, (Pose2, Point2), zi, (xi, li))

@show res
@test norm( res ) < 1e-14

##

Xi = zeros(3); Xi[3] = pi/2
xi = getPoint(Pose2, Xi)
li = zeros(2); li[2] = 20.0;
_zi = [0.0,20.0]
zi = Manifolds.hat(M, identity_element(M), _zi)

res = calcFactorResidualTemporary( p2br, (Pose2, Point2), zi, (xi, li))

#
@show res
@test norm(res) < 1e-14

##

Xi = zeros(3); Xi[3] = -pi/2
xi = getPoint(Pose2, Xi)
li = zeros(2); li[1] = 20.0;
# zi = ([0.0;pi/2],[0.0;20.0],)
_zi = [pi/2,20.0]
zi = Manifolds.hat(M, identity_element(M), _zi)
# idx = 2
# res = zeros(2)
# p2br(res, fmd, idx, zi, xi, li)

res = calcFactorResidualTemporary( p2br, (Pose2, Point2), zi, (xi, li))

@show res
@test norm(res) < 1e-14

##
x1 = ProductRepr([0.,0], [1. 0; 0 1])
x2 = ProductRepr([0.,0], [0 -1.; 1 0])

#measurement setup 1
meas = (0., 10)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [10.,0]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)
# x2
p = x2
q = [0., 10]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)

meas = (pi/2., 10)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [0.,10]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)
# x2
p = x2
q = [-10., 0]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)

meas = (pi, 10.)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [-10.,0]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)
# x2
p = x2
q = [0., -10]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)

meas = (-pi/2., 10.)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [0.,-10]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)
# x2
p = x2
q = [10., 0]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,0], atol = 1e-9)

##
# testing non zero errors on range
#FIXME BR range sign is broken, needed for gradients

meas = (0., 10)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [11., 0]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test_broken isapprox(res, [0,-1], atol = 1e-9)
# x2
p = x2
q = [0., 11]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test_broken isapprox(res, [0,-1], atol = 1e-9)

meas = (0., 10)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [9., 0]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,1], atol = 1e-9)
# x2
p = x2
q = [0., 9]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res, [0,1], atol = 1e-9)

# on small angles also
s,c = 10 .* sincos(0.001)
meas = (0., 10)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [c, s]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res[1], -0.001, atol = 1e-9)
#FIXME ?
@test isapprox(res[2], 0, atol = 0.1)
# x2
p = x2
q = [s, c]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test isapprox(res[1], 0.001, atol = 1e-9)
#FIXME ?
@test isapprox(res[2], 0, atol = 0.1)

# WIP testing non zero errors
# I don't know if this test is needed or even possible
r2 = 10/sqrt(2)
meas = (0., 10)
f = Pose2Point2BearingRange(Normal(meas[1],1), Normal(meas[2],1))
X = hat(M, Identity(M), [meas[1],meas[2]])
# x1
p = x1
q = [r2, r2]
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test_broken isapprox(res, [-pi/4,0], atol = 1e-9)
# x2
p = x2
res = calcFactorResidualTemporary( f, (Pose2, Point2), X, (p, q))
@test_broken isapprox(res, [pi/4,0], atol = 1e-9)

##
end

@testset "test unimodal bearing range factor, solve for landmark..." begin

##

# Start with an empty graph
N = 1
fg = initfg()

#add pose with partial constraint
addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*Matrix{Float64}(LinearAlgebra.I, 3,3))), graphinit=false)
# force particular initialization
setVal!(fg, :x0, [getPointIdentity(Pose2)])

##----------- sanity check that predictbelief plumbing is doing the right thing
_pts, = predictbelief(fg, :x0, ls(fg, :x0), N=75)
@cast pts[j,i] := DFG.getCoordinates.(Pose2, _pts)[i][j]
@test sum(abs.(Statistics.mean(pts,dims=2)) .< [0.1; 0.1; 0.1]) == 3
@test sum([0.05; 0.05; 0.05] .< Statistics.std(pts,dims=2) .< [0.15; 0.15; 0.15]) == 3
#------------

# Add landmark
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
li = zeros(2); li[1] = 20.0;
setVal!(fg, :l1, [li])

# Add bearing range measurement between pose and landmark
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br, graphinit=false)

# there should be just one (the bearingrange) factor connected to :l1
@test length(ls(fg, :l1)) == 1
# drawGraph(fg, show=true)

# check the forward convolution is working properly
_pts, = predictbelief(fg, :l1, ls(fg, :l1), N=75)
@cast pts[j,i] := _pts[i][j]
@show tp = mean(TranslationGroup(2), _pts)
@warn "weak test tolerance, suspect partial products need to be upgraded first.  Please see likely AMP #41 and IIF #1010 for known issues likely the root cause."
@test isapprox( tp, [20.0; 0.0], atol=5.0 )
@test sum([0.1; 0.1] .< Statistics.std(pts,dims=2) .< [3.0; 3.0]) == 2

# using Gadfly, KernelDensityEstimate, KernelDensityEstimatePlotting
#
# pl = plotKDE(kde!(pts))
# pl.coord = Coord.Cartesian(xmin=-5,xmax=25, ymin=-10.0,ymax=10)
# pl

##

end

@testset "test unimodal bearing range factor, solve for pose..." begin

##

# Start with an empty graph
N = 1
fg = initfg()

# Add landmark
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
addFactor!(fg, [:l1], PriorPoint2(MvNormal([20.0;0.0], Matrix(Diagonal([1.0;1.0].^2)))),  graphinit=false ) # could be IIF.Prior
li = zeros(2); li[1] = 20.0;
setVal!(fg, :l1, [li])

#add pose with partial constraint
addVariable!(fg, :x0, Pose2)
# force particular initialization
setVal!(fg, :x0, [getPointIdentity(Pose2)])

# Add bearing range measurement between pose and landmark
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br, graphinit=false)

# there should be just one (the bearingrange) factor connected to :l1
@test length(ls(fg, :x0)) == 1
# writeGraphPdf(fg)

# check the forward convolution is working properly
_pts, = predictbelief(fg, :x0, ls(fg, :x0), N=75)
p_μ = mean(SpecialEuclidean(2), _pts)

_pts = getCoordinates.(Pose2, _pts)
@cast pts[j,i] := _pts[i][j]

# pts[3,:] .= TU.wrapRad.(pts[3,:])
@show abs.(Statistics.mean(pts,dims=2))
@test sum(abs.(Statistics.mean(pts,dims=2)) .< [2.0; 2.0; 2.0]) == 3
@show Statistics.std(pts,dims=2)
@test sum([0.1; 2.0; 0.01] .< Statistics.std(pts,dims=2) .< [5.0; 10.0; 2.0]) == 3

##

end

@testset "Testing Pose2Point2Bearing Initialization and Packing" begin

##

p2p2b = Pose2Point2Bearing( MvNormal([0.2,0.2,0.2], [1.0 0 0;0 1 0;0 0 1]) )
packed = convert(PackedPose2Point2Bearing, p2p2b)
p2p2bTest = convert(Pose2Point2Bearing, packed)
@test p2p2b.bearing.μ == p2p2bTest.bearing.μ
@test p2p2b.bearing.Σ.mat == p2p2bTest.bearing.Σ.mat

##

end

#=

fg = LightDFG{SolverParams}(solverParams=SolverParams())

# Add the first pose :x0
x0 = addVariable!(fg, :x0, Pose2)

p = getPoint(Pose2, [10; 10; -pi])

prior = addFactor!(fg, [:x0], PriorPose2( MvNormal([0.1,0.1,0.05]), p ) )

for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    addFactor!(fg, [psym;nsym], pp )
end

addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br)

p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x6; :l1], p2br)

smtasks = Task[]
solveTree!(fg; smtasks)
#add bearing range

## ================================================================================================
## Other tests
## ================================================================================================
hist = hists[3]
step = 7
fnc_ = hist[step].f
# the data before step
csmc_ = deepcopy(hist[step].csmc);
csmc_.enableLogging = false;
csmc_.logger =  SimpleLogger(Base.devnull);#NullLogger()

##

fg = initfg()

pRight = 0.9
pWrong = 0.1
pr_noise = [0.01, 0.01, 0.001]
od_noise = [0.2; 0.2; 0.2]
lm_noise = [0.01, 0.01, 0.001]

#x0 prior
addVariable!(fg, :x0, Pose2)
prpo = MvNormal([0,0,-pi], pr_noise)
addFactor!(fg, [:x0], PriorPose2(MvNormal(rand(prpo), pr_noise)))

#l1 and l2
addVariable!(fg, :l1, Pose2, tags=[:LANDMARK])
addVariable!(fg, :l2, Pose2, tags=[:LANDMARK])

#x0 to l1 or l2
p2ln = MvNormal([1, 1, pi], lm_noise)
p2p = Pose2Pose2(MvNormal(rand(p2ln), lm_noise))
addFactor!(fg, [:x0; :l1; :l2], p2p, multihypo = [1.0, pRight, pWrong])

#x0 to x1
addVariable!(fg, :x1, Pose2)
pp = MvNormal([1.0,0,0], od_noise)
addFactor!(fg, [:x0,:x1], Pose2Pose2(MvNormal(rand(pp), od_noise)))

#x1 to l1 or l2
p2ln = MvNormal([0, 1, pi], lm_noise)
p2p = Pose2Pose2(MvNormal(rand(p2ln), lm_noise))
addFactor!(fg, [:x1; :l1; :l2], p2p, multihypo = [1.0, pRight, pWrong])

#x1 to l2 or l1
p2ln = MvNormal([0, -1, pi], lm_noise)
p2p = Pose2Pose2(MvNormal(rand(p2ln), lm_noise))
addFactor!(fg, [:x1; :l2; :l1], p2p, multihypo = [1.0, pRight, pWrong])

#aslo good example without x2
#x1 to x2
addVariable!(fg, :x2, Pose2)
pp = MvNormal([1.0,0,0], od_noise)
addFactor!(fg, [:x1,:x2], Pose2Pose2(MvNormal(rand(pp), od_noise)))

#x2 to l1 or l2
p2ln = MvNormal([-1, 1, pi], lm_noise)
p2p = Pose2Pose2(MvNormal(rand(p2ln), lm_noise))
addFactor!(fg, [:x2; :l1; :l2], p2p, multihypo = [1.0, pRight, pWrong])

#x2 to l2 or l1
p2ln = MvNormal([-1, -1, pi], lm_noise)
p2p = Pose2Pose2(MvNormal(rand(p2ln), lm_noise))
addFactor!(fg, [:x2; :l2; :l1], p2p, multihypo = [1.0, pRight, pWrong])

##

fg.solverParams.inflation=0.1
fg.solverParams.spreadNH=0.1
solveTree!(fg)

plotSLAM2D(fg)

##

# Random.seed!(42) # The answer to reproducable noise

fg = LightDFG(solverParams=SolverParams(graphinit=false, gibbsIters=5, spreadNH=5.0))

pRight = 0.8
pWrong = 0.2
pr_noise = 0.01
od_noise = 0.1
lm_noise = 0.01

# true positions
# x0 at 0
# x1 at 1
# l1 at 1
# l2 at 2

#x0 prior
addVariable!(fg, :x0, ContinuousScalar)
prpo = Normal(0.0, pr_noise)
addFactor!(fg, [:x0], Prior(Normal(rand(prpo), pr_noise)))

#l1 and l2
addVariable!(fg, :l1, ContinuousScalar, tags=[:LANDMARK])
addVariable!(fg, :l2, ContinuousScalar, tags=[:LANDMARK])

#x0 to l1 or l2
p2ln = Normal(1.0, lm_noise)
p2p = LinearRelative(Normal(rand(p2ln), lm_noise))
addFactor!(fg, [:x0; :l1; :l2], p2p, multihypo = [1, pRight, pWrong])
# addFactor!(fg, [:x0; :l1], p2p) #this one used for sanity check

#x0 to x1
addVariable!(fg, :x1, ContinuousScalar)
pp = Normal(1.0, od_noise)
addFactor!(fg, [:x0,:x1], LinearRelative(Normal(rand(pp), od_noise)))

#x1 to l1 or l2
p2ln = Normal(0.0, lm_noise)
p2p = LinearRelative(Normal(rand(p2ln), lm_noise))
addFactor!(fg, [:x1; :l1; :l2], p2p, multihypo = [1, pRight, pWrong])
# addFactor!(fg, [:x1; :l1], p2p) #this one used for sanity check

#x1 to l2 or l1
p2ln = Normal(1.0, lm_noise)
p2p = LinearRelative(Normal(rand(p2ln), lm_noise))
addFactor!(fg, [:x1; :l2; :l1], p2p, multihypo = [1, pRight, pWrong])
# addFactor!(fg, [:x1; :l2], p2p) #this one used for sanity check

#

# prescribe an elimination order to get a single clique
eo = [:l2,:x1,:x0,:l1]
# fg.solverParams.graphinit=true

fg.solverParams.inflation=3.0
fg.solverParams.spreadNH=1.0
fg.solverParams.dbg = true
##
smtasks = Task[]
tree, _, = solveTree!(fg; smtasks, eliminationOrder=eo) #, smtasks=smtasks, recordcliqs=ls(fg));

# hists = fetchCliqHistoryAll!(smtasks)

plotKDE(fg, ls(fg))

##

@test isapprox(DFG.getPPESuggested(fg, :x0)[], 0, atol = 0.2)
@test isapprox(DFG.getPPESuggested(fg, :x1)[], 1, atol = 0.2)
@test isapprox(DFG.getPPESuggested(fg, :l1)[], 1, atol = 0.2)

L2 = getBelief(fg, :l2)
L2_ = manikde!(2 .+ 0.1*randn(size(getPoints(L2),2)), ContinuousScalar)

# test that there is at least a mode present
mmd(L2_, L2, ContinuousScalar)
@test isapprox(DFG.getPPESuggested(fg, :l2)[], 2, atol = 0.2)

=#


text:
using FixedPointNumbers, Test

"""
## Reference
[1] C. S. Turner,  "A Fast Binary Logarithm Algorithm", IEEE Signal
     Processing Mag., pp. 124,140, Sep. 2010.
"""
function log2fix(x::Fixed{T, P}) where {T, P}
    PREC = UInt(P)
    x.i == 0 && return typemin(T) # represents negative infinity

    y = zero(T)
    xi = x.i
    while xi < 1 << PREC
        xi <<= 1
        y -= T(1) << PREC
    end

    while xi >= 2 << PREC
        xi >>= 1
        y += T(1) << PREC
    end

    z = xi
    b = T(1) << (PREC - UInt(1))
    for i = 1:P
        temp = Base.widemul(z, z) >> PREC
        z = T(temp)
        if z >= T(2) << PREC
            z >>= 1
            y += b
        end
        b >>= 1
    end

    return Fixed{T,PREC}(y, nothing)
end

@test log2fix(Fixed{Int, 43}(2^1.24)) ≈ 1.24

ext: jl


text:
@reexport module CanonicalSymplecticStructures

import Base
import LinearAlgebra

export CanonicalSymplecticMatrix

struct CanonicalSymplecticMatrix{T} <: AbstractMatrix{T}
    mid::Int
    function CanonicalSymplecticMatrix{T}(n::Integer) where T
        iseven(n) || throw(ArgumentError("CanonicalSymplecticMatrix must have even size"))
        new{T}(n ÷ 2)
    end
end

CanonicalSymplecticMatrix(n::Integer) = CanonicalSymplecticMatrix{Int}(n)

Base.size(C::CanonicalSymplecticMatrix) = (sz = C.mid * 2; (sz, sz))

function Base.getindex(C::CanonicalSymplecticMatrix{T}, i1::Int, i2::Int) where T
    if i2 - C.mid == i1
        return T(-1)
    elseif i1 - C.mid == i2
        return T(1)
    else
        return T(0)
    end
end

function LinearAlgebra.dot(x::AbstractVector, C::CanonicalSymplecticMatrix, y::AbstractVector)
    m = C.mid
    axes(x, 1) == axes(y, 1) == 1:2m || error()
    s = zero(promote_type(eltype(x), eltype(y)))
    @inbounds for i in 1:m
        s += y[i] * x[m + i]
        s -= y[m + i] * x[i]
    end

    s
end

end  # module


text:
# src/comparison.jl
# Comparisons between Rationals and other Numeral Types

import Base: ==, <, <=

"""```
x::AbstractRatio == y::AbstractRatio
```
Check if the ratio `x` is equal to the ratio `y`.
"""
==(x::AbstractRatio, y::AbstractRatio) = ==(crossterms(x, y)...)

"""```
x::AbstractRatio == y::AbstractRational
```
Check if the ratio `x` is equal to the rational `y`.
"""
==(x::AbstractRatio, y::AbstractRational) = ==(crossterms(x, y)...)

"""```
x::AbstractRational == y::AbstractRational
```
Check if the rational `x` is equal to the rational `y`.
"""
function ==(x::AbstractRational, y::AbstractRational)
    return (denominator(x) == denominator(y)) & (numerator(x) == numerator(y))
end

"""```
x::AbstractRatio == y::Integer
```
Check if the ratio `x` is equal to the integer `y`.
"""
==(x::AbstractRatio, y::Integer) = ==(crossterms(x, y)...)

"""```
x::AbstractRational == y::Integer
```
Check if the rational `x` is equal to the integer `y`.
"""
==(x::AbstractRational, y::Integer) = isone(denominator(x)) & (numerator(x) == y)

"""```
x::AbstractRational == y::AbstractFloat
```
Check if the rational `x` is equal to the float `y`.
"""
function ==(x::AbstractRational, y::AbstractFloat)
    nx, dx = numden(x)
    if isfinite(y)
        return isone(count_ones(dx)) & (y * dx == nx)
    else
        return y == nx / dx
    end
end

"""```
x::AbstractRatio == z::Complex
```
Check if the ratio `x` is equal to the complex number `z`.
"""
==(x::AbstractRatio, z::Complex) = isreal(z) & (x == real(z))

"""```
x == y::AbstractRatio
    === (y == x)
```
Check if `x` is equal to the ratio `y`.
"""
==(x, y::AbstractRatio) = y == x

"""```
x::AbstractRatio < y::AbstractRatio
```
Check if the ratio `x` is less than the ratio `y`.
"""
function <(x::AbstractRatio, y::AbstractRatio)
    nx, dx = numden(x)
    ny, dy = numden(y)
    return dx == dy ? nx < ny : widemul(nx, dy) < widemul(ny, dx)
end

"""```
x::AbstractRatio < y
```
Check if the ratio `x` is less than `y`.
"""
<(x::AbstractRatio, y) = numerator(x) < widemul(denominator(x), y)

"""```
x < y::AbstractRatio
```
Check if `x` is less than the ratio `y`.
"""
<(x, y::AbstractRatio) = widemul(x, denominator(y)) < numerator(y)

"""```
x::AbstractRatio <= y::AbstractRatio
```
Check if the ratio `x` is less than or equal to the ratio `y`.
"""
function <=(x::AbstractRational, y::AbstractRational)
    nx, dx = numden(x)
    ny, dy = numden(y)
    return dx == dy ? nx <= ny : widemul(nx, dy) <= widemul(ny, dx)
end

"""```
x::AbstractRatio <= y
```
Check if the ratio `x` is less than or equal to `y`.
"""
<=(x::AbstractRatio, y) = numerator(x) <= widemul(denominator(x), y)

"""```
x <= y::AbstractRatio
```
Check if `x` is less than or equal to a ratio `y`.
"""
<=(x, y::AbstractRatio) = widemul(x, denominator(y)) <= numerator(y)


text:
# switch statement

macro switch(value, block)
  @gensym value_sym

  :(for i in 1:1
      let $(esc(value_sym)) = $value
          $(esc(switch(value_sym, block)))
      end
    end)
end

@metafunction switch(value, begin :L{:case}(x); :?{xs}; :L{:case}(y); *{ys} end) begin
  next = switch(value, quote :case($y); $(ys...) end)
  switch_clause(value, x, xs, next)
end

@metafunction switch(value, begin :L{:case}(x); *{xs} end) begin
  switch_clause(value, x, xs)
end

switch_clause(value, x, xs, next=nothing) = quote
  if $value == $x
    $(xs...)
  else
    $next
  end
end

r = 0
a = 2

@switch a begin
  :case(1);
    r = 1

  :case(2);
    r = 2

end

if a in [1,2]
  @test a != 1 || r == 1
  @test a != 2 || r == 2
end


text:
#Standard Julia Capabilities

using RDatasets
describe(dataset("datasets", "iris"))


text:
# example uniaxial extension
push!( LOAD_PATH, "../src")

using Materialis

# reference material parameters
pRef             = 2
nu               = 0.4

largo = 1
ladoSecc = .2

Lx = largo; Ly = ladoSecc ; Lz = Ly;

nodalCoords = [ 0    0    0 ;
                0    0   Lz ;
                0   Ly   Lz ;
                0   Ly    0 ;
                Lx   0    0 ;
                Lx   0   Lz ;
                Lx  Ly   Lz ;
                Lx  Ly    0 ]

# the connectivity matrix is given by a vector of vectors.
elemNodalConnec =[[ 1 4 2 6 ] ,
                  [ 6 2 3 4 ] ,
                  [ 4 3 6 7 ] ,
                  [ 4 1 5 6 ] ,
                  [ 4 6 5 8 ] ,
                  [ 4 7 6 8 ]  ]

nLoadSteps       = 1
levelerror       = 0.0

myGrid = gridData( 100, 3, 3, .1*largo, .9*largo, ladoSecc*.1, ladoSecc*.5 , ladoSecc*.1, ladoSecc*.5 )
myROI  = ROIData( .4, .6, ladoSecc*0.25, ladoSecc*0.35, ladoSecc*0.25, ladoSecc*0.35    )

plotParamsVector = [ 3 ]
tolDeltaP        = 1e-4 # convergence criterion
tolDeltaP = 1e-5


text:
function noise_threshold(x,th = 0.4)
    if x<=th
        return 0
    else
        return 1
    end
end
"""
Pixels with no economic activity may show some light due to background noise. These pixels could be in forests, oceans, deserts etc. The ```background_noise_mask``` function generates a background moise mask such that those pixels which are considered dark are marked as 0 and those considered lit are marked as 1. The function uses the datacubes of radiance and clouds to generate annual image of the last year the data. The function considers all the pixels below a provided threshold as dark and remaining to be lit.
```julia
background_noise_mask(radiance_datacube, clouds_datacube)
```
"""
function background_noise_mask(radiance_datacube, clouds_datacube, th = 0.4)
    # This function may be obsolete because Payne Institute is providing annual images for each year.
    last_year_rad      = radiance_datacube[:, :, (size(radiance_datacube)[3]-11):size(radiance_datacube)[3]]
    last_year_cloud   = clouds_datacube[:, :, (size(radiance_datacube)[3]-11):size(radiance_datacube)[3]]
    average_lastyear = copy(radiance_datacube[:, :, 1])
    for i in 1:size(last_year_rad)[1]
        for j in 1:size(last_year_rad)[2]
            average_lastyear[i,j] = weighted_mean(last_year_rad[i, j, :], last_year_cloud[i, j, :])
        end
    end
    mask = noise_threshold.(average_lastyear, th)
    return mask
end


text:
dir = "/Users/jacobquinn/csv-comparison/data/"

for d in (
    "rows_100_cols_20_na_false",
    "rows_100_cols_20_na_true",
    "rows_100_cols_200_na_false",
    "rows_100_cols_200_na_true",
    "rows_10000_cols_20_na_false",
    "rows_10000_cols_20_na_true",
    "rows_10000_cols_200_na_false",
    "rows_10000_cols_200_na_true",
    "rows_1000000_cols_20_na_false",
    "rows_1000000_cols_20_na_true",
    "rows_1000000_cols_200_na_false",
    "rows_1000000_cols_200_na_true",
), f in (
    "mixed_data.csv",
    "uniform_data_catstring.csv",
    "uniform_data_float64.csv",
    "uniform_data_shortfloat64.csv",
    "mixed_data_shortfloat64.csv",
    "uniform_data_datetime.csv",
    "uniform_data_int64.csv",
    "uniform_data_string.csv",
)
    file = joinpath(dir, d, f)
    println("file=$file")
    GC.enable(false)
    @time CSV.read(file)
    GC.enable(true)
    GC.gc(); GC.gc()
    GC.enable(false)
    @time TextParse.csvread(file)
    GC.enable(true)
    GC.gc(); GC.gc()
end

file = "/Users/jacobquinn/csv-comparison/data/rows_1000000_cols_20_na_false/uniform_data_int64.csv"
@time CSV.read(file)


text:
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# PROGRAM NAME: main.jl
# DESCRIPTION:  This program reads the most recent data from
#               the CDC and plots for data visualization purposes
# AUTHOR:       Isaac Weintraub Jan 2021
# VERSION	0.1
# REQ.PACKAGES:	CSV
#
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# using Pkg
# cd(@__DIR__)
# Pkg.activate(".")
# Pkg.instantiate()

;cd("/Users/weintraub0/Documents/Projects/COVID19/COVID-19")
;run(`git pull`)    # Get the latest data
;cd("/Users/weintraub0/Documents/Projects/COVID19/")
using CSV           # Need to read the CSV data
using Plots         # Allows the genation of plots]
using PlotlyJS
plotlyjs()

using BenchmarkTools
using DataFrames
using DelimitedFiles

#ENV["MPLBACKEND"]="agg"
#pygui(true)
# Focusing on how to read data into julia
# What are the common data structures

#using Pandas
include("covidPlot2.jl")

# # Read the CSV Data (ARRAY)
# # This is the method that you would want to use if you don't want to use dataframes (slower than dataframes)
# P_popStates,H_popStates = readdlm("USA.csv",',',;header=true)
# P_casesUSA,H_casesUSA = readdlm("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv",',',;header=true)
# P_deathUSA,H_deathUSA = readdlm("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv",',',;header=true)

# Read the CSV Data (DATAFRAME)
# This is the method you should use if you want to use dataframes (faster than the array method)
DFpopStates = CSV.read("USA.csv", DataFrame)    # Read the population data from the united states
DFcasesUSA  = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv", DataFrame)
DFdeathUSA  = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv", DataFrame)

# DFpopStates[!,:Province] # Prints the list of the provinces in the dataframe
# DFpopStates.Province     # Another way to list the provinces in the dataframe
# names(DFpopStates)       # prints the names for the columns in the dataframe
# describe(DFpopStates)    # prints a summary of the data in compact way

function getPopulation(P_df,Name::String)
# This fuction returns the populatino of a Provice (String) to a provided DataFrame
      loc = findfirst(P_df.Province .== Name)
      !isnothing(loc) && return P_df.Population[loc]
      error("Error: Province Not Found")
end

function getIndex(DF,Name::String)
      loc = findall(DFcasesUSA.Province_State .== Name)
      !isnothing(loc) && return DFcasesUSA.Province_State .== Name
      error("Error: Province Not Found")
end

numProvinces = length(DFpopStates.Province);          # Get the total number of provinces

p1 = nothing
for name in DFpopStates.Province             # Run a for loop over the states
      global p1
      population = getPopulation(DFpopStates,name)   # Population
      indexCasesUSA = getIndex(DFcasesUSA,name)      # Index (Cases)
      indexDeathUSA = getIndex(DFdeathUSA,name)      # Index (Deaths)
      cases = DFcasesUSA[indexCasesUSA,12:end]       # cases
      death = DFdeathUSA[indexDeathUSA,12:end]       # deaths
      casesProvince = sum(Array(cases),dims=1)       # cases  (Array)
      deathProvince = sum(Array(death),dims=1)       # deaths (Array)
      println("$name Cases: $(casesProvince[end]) Deaths: $(deathProvince[end])")   #     Message of the current state
      p1 = plot()
      hspan!(p1,[0,1e2], color  = :red, alpha   = 0.6, lab = nothing);
      hspan!(p1,[0,1e1], color  = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e0], color  = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e-1], color = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e-2], color = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e-3], color = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e-4], color = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e-5], color = :white, alpha = 0.1, lab = nothing);
      hspan!(p1,[0,1e-6], color = :white, alpha = 0.1, lab = nothing);
      plot!(casesProvince[1,2:end]./population.*100,
            lab="Cases",
            #yscale=:log10,
            xscale=:identity,
            color=:black)
      plot!(deathProvince[1,2:end]./population.*100,
            lab="Deaths",
            #yscale=:log10,
            xscale=:identity,
            legend=:bottomright,
            color=:blue,
            yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
            yaxis="Percent of Population")
      title!(name)
      annotate!((size(casesProvince)[2]+5,
                  casesProvince[1,end]./population.*100,
                  "$(casesProvince[end])"))
      annotate!((size(deathProvince)[2]+5,
                  deathProvince[1,end]./population.*100,
                  "$(deathProvince[end])"))

      ylims!((1e-6,100))
      xlims!((1,90+length(casesProvince)))
      plot!(yscale=:log10)
      savefig(p1,"Figures/UnitedStates/$name.png")
      gui()
end

# USA PLOT

# data_popStates.Province .== "Ohio"
# n = size(data_popStates.Province,1); # Number to search through

# state = "Ohio"
# println(state)
# casesDF = data_casesUSA[data_casesUSA.Province_State .== state,:]   # cases
# deathDF = data_deathUSA[data_deathUSA.Province_State .== state,:]   # deaths
# popDF   = data_popStates[data_popStates.Province .== state,:]
# cases = casesDF[:,12:end]
# death = deathDF[:,12:end]
# popState   = popDF[1,2]
# casesArray = convert(Array, cases[1:end,:])
# casesState = sum(casesArray,dims=1)
# deathArray = convert(Array, cases[1:end,:])
# deathState = sum(deathArray,dims=1)

# cases = casesState
# death = deathState
# pop = popState
# name = state
# p1 = plot()
# hspan!(p1,[0,1e2], color  = :red, alpha   = 0.6, labels = false);
# hspan!(p1,[0,1e1], color  = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e0], color  = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e-1], color = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e-2], color = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e-3], color = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e-4], color = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e-5], color = :white, alpha = 0.1, labels = false);
# hspan!(p1,[0,1e-6], color = :white, alpha = 0.1, labels = false);
# scatter!(1:size(cases)[2],
#       100*cases[1,:]./pop,lab=name*" COVID Cases",
#       scale=:log10,
#       yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#       yaxis="Percent of Country", xaxis="Days since Jan 22",
#       marker =:o,
#       color =:green)

# scatter!(1:size(death)[2],
#       100*death[1,|:]./pop,lab=name*" COVID Deaths",
#       scale=:log10,
#       yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#       yaxis="Percent of Country", xaxis="Days since Jan 22",
#       marker =:x,
#       color =:black)

#  plot!(xscale=:identity,xticks=0:100:500, legend=:bottomright)
#  labelCases = cases[1,end];
#  labelDeath = death[1,end];
#  annotate!(size(cases)[2]+5, 200*cases[1,end]./pop, ha="left",va="bottom","$labelCases")
#  annotate!(size(cases)[2]+5, 200*death[1,end]./pop, ha="left",va="bottom","$labelDeath")
#  xlims!((0,500))
#  ylims!((1e-6,100))
#  png(p1,"Figures/UnitedStates/$name.png")

# for state in data_popStates.Province
#     println(state)
#     casesDF = data_casesUSA[data_casesUSA.Province_State .== state,:]   # cases
#     deathDF = data_deathUSA[data_deathUSA.Province_State .== state,:]   # deaths
#     popDF   = data_popStates[data_popStates.Province .== state,:]
#     cases = casesDF[:,12:end]
#     death = deathDF[:,12:end]
#     popState   = popDF[1,2]
#     casesArray = convert(Array, cases[1:end,:])
#     casesState = sum(casesArray,dims=1)
#     deathArray = convert(Array, cases[1:end,:])
#     deathState = sum(deathArray,dims=1)
#
#     cases = casesState
#     death = deathState
#     pop = popState
#     p1 = plot()
#     hspan!(p1,[0,1e2], color  = :red, alpha   = 0.6, labels = false);
#     hspan!(p1,[0,1e1], color  = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e0], color  = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-1], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-2], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-3], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-4], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-5], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-6], color = :white, alpha = 0.1, labels = false);
#     scatter!(1:size(cases)[2],
#           100*cases[1,:]./pop,lab=name*" COVID Cases",
#           scale=:log10,
#           yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#           yaxis="Percent of Country", xaxis="Days since Jan 22",
#           marker =:o,
#           color =:green)
#
#     scatter!(1:size(cases)[2],
#           100*death[1,:]./pop,lab=name*" COVID Deaths",
#           scale=:log10,
#           yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#           yaxis="Percent of Country", xaxis="Days since Jan 22",
#           marker =:x,
#           color =:black)
#      plot!(xscale=:identity,xticks=0:100:500, legend=:bottomright)
#      labelCases = cases[1,end];
#      labelDeath = death[1,end];
#      annotate!(size(cases)[2]+5, 200*cases[1,end]./pop, ha="left",va="bottom","$labelCases")
#      annotate!(size(cases)[2]+5, 200*death[1,end]./pop, ha="left",va="bottom","$labelDeath")
#      xlims!((0,500))
#      ylims!((1e-6,100))
#      png(p1,"Figures/UnitedStates/$name.png")
#
#
#     covidPlot2(casesState,deathState,popState,state)
# end

#data_casesUSA   = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv",header=true)       # cases of COVID 19 in the USA
#data_deathUSA   = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv",header=true)          # deaths in the USA from COVID 19
#data_casesGlobe = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv",header=true)   # casess of COVID 19 around the globe
#data_deathGlobe = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv",header=true)      # deaths in the Globe from COVID 19
#data_recovGlobe = CSV.read("COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv",header=true)   # recovered global COVID 19 Cases
#data_populGlobe = CSV.read("population/data/population.csv",header=true) # get the population data
# data_casesUSAstateList = data_casesUSA[:,7]       # Gather a list of all the states
# countryList = unique(data_populGlobe[:,1]) # Remove Duplicates
#
# n = size(countryList)[1]
# global casesC
# global deathC
# global recovC
# @showprogress 0.01 for i = 1:n              # show a progress meter and for loop through the countries
#     countryName = countryList[i]
#         year = maximum(data_populGlobe[data_populGlobe[:,1].==countryName,3])
#         select = data_populGlobe[:,1].==countryName
#         index = data_populGlobe[select,3].==year
#         selectPopulations = data_populGlobe[select,4]
#     countryPopulation = selectPopulations[index][1]
#     if countryName == "United States"
#         countryName = "US"
#     end
#     if sum(data_casesGlobe[:,2] .== countryName) > 0
#         #print(countryName*"\n")
#         index = data_casesGlobe[:,2] .== countryName
#
# 		cases = data_casesGlobe[index,5:end]
# 		death = data_deathGlobe[index,5:end]
# 		#recov = data_recovGlobe[index,5:end]
#
#         casesMatrix = convert(Array, cases[1:end,:])
#         deathMatrix = convert(Array, death[1:end,:])
# 		#recovMatrix = convert(Array, recov[1:end,:])
#
# 		countryCases = sum(casesMatrix,dims=1)
# 		countryDeath = sum(deathMatrix,dims=1)
# 		#countryRecov = sum(recovMatrix,dims=1)
# 		#covidPlot1(countryCases,countryPopulation,countryName) 	# Plot the Covid Data
# 		covidPlot2(countryCases,countryDeath,countryPopulation,countryName)
#     end
# end
#
#
# n = size(data_popStates)[1]                  # the number of provinces and states of the USA
# global USACases
# global USApop
# global USADeaths
# @showprogress 0.01 for i = 1:n               # show a progress meter and for loop through the states
#     stateName        = data_popStates[i,1]          # get the name of the current state (STRING)
#     statePopulation  = data_popStates[i,2]          # get the population of the current state (FLOAT)
#     index            = data_casesUSA[:,7] .== stateName   # get the subset of indexies whicih corespond to data for that state
#     stateCases       = data_casesUSA[index,12:end]  # parse the subset of the data for the state of interest
#     stateDeaths      = data_deathUSA[index,12:end]             # Gather the deaths
#     stateCaseMatrix  = convert(Array, stateCases[1:end,:])              # Conver the data to a matrix for manipulation reasons
#     stateDeathMatrix = convert(Array, stateDeaths[1:end,:])    # Convert the deat data to a matrix for manipulation reasons
#     stateCases       = sum(stateCaseMatrix,dims=1)              # summing is easier as a matrix than as a data frame NEED THE TOTAL NUMBER OF CASES in each state
#     stateDeathCases  = sum(stateDeathMatrix,dims=1)    # summing of the deaths
#     global USACases
#     global USApop
#     global USADeaths
#     if i == 1
#         USACases  = stateCases
#         USADeaths = stateDeathCases
#         USApop    = statePopulation
#     elseif i > 2
#         USACases  = USACases .+ stateCases
#         USADeaths = USADeaths .+ (stateDeathCases)
#         USApop    = USApop .+ statePopulation
#     end
#     # Plot the current state or provices
#     p1 = plot()
#     hspan!(p1,[0,1e2], color  = :red, alpha   = 0.6, labels = false);
#     hspan!(p1,[0,1e1], color  = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e0], color  = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-1], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-2], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-3], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-4], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-5], color = :white, alpha = 0.1, labels = false);
#     hspan!(p1,[0,1e-6], color = :white, alpha = 0.1, labels = false);
#     scatter!(1:size(stateCases)[2],
#         100*stateCases[1,:]./statePopulation,lab=stateName*" COVID Cases",
#         scale=:log10,
#         yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#         yaxis="Percent of State", xaxis="Days since Jan 22",
#         marker =:o,
#         color =:red)
#     scatter!(1:size(stateDeathCases)[2],
#         100*stateDeathCases[1,:]./statePopulation,lab=stateName*" COVID Deaths",
#         scale=:log10,
#         yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#         yaxis="Percent of State", xaxis="Days since Jan 22",
#         marker =:x,
#         color=:black)
#     plot!(xscale=:identity,xticks=0:10:100)
#     xlims!((0,100))
#     ylims!((1e-6,100))
#     png(p1,"StateFigures/$stateName.png")
# end
# # PLOT THE UNITED STATES PLOT FOR ALL THE CASES IN THE USA USING STATE/PROVINCE DATA
# p2 = plot()
# hspan!(p2,[0,1e2], color  = :red, alpha   = 0.6, labels = false);
# hspan!(p2,[0,1e1], color  = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e0], color  = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e-1], color = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e-2], color = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e-3], color = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e-4], color = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e-5], color = :white, alpha = 0.1, labels = false);
# hspan!(p2,[0,1e-6], color = :white, alpha = 0.1, labels = false);
# scatter!(1:size(USACases)[2],100*USACases[1,:]./USApop,lab="USA",
#      scale=:log10,
#      yticks = [1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,1e1,1e2],
#      yaxis="Percent of USA Pop", xaxis="Days since Jan 22",
#      marker=:o)
# plot!(xscale=:identity,xticks=0:10:100)
# xlims!((0,100))
# ylims!((1e-6,100))
# png(p2,"StateFigures/USA.png")
# print("Done")


text:
# Exercise 1

using LinearAlgebra, Statistics

function lyapunov(A,Σ; tolerance = 1.0E-8, maxiter = 500, S0 = Σ*Σ')

    #initialise loop
    normdiff = Inf
    iter = 1
    S_old = S0
    Q = Σ * Σ'

    while normdiff > tolerance && iter ≤ maxiter

        S_new = A * S_old * A' + Q

        #continue looping
        normdiff = norm(S_old - S_new)
        S_old = S_new
        iter = iter +1
    end
    return (S_old,iter)
end

# test function
A = [0.8  -0.2; -0.1  0.7]
Σ = [0.5 0.4; 0.4 0.6]

@show maximum(abs, eigvals(A))

our_solution, iter = lyapunov(A, Σ)

###############
# Exercise 2
###############
using Plots

Θ = [0.8, 0.9, 0.98] #parameters
γ = 1
σ = 1
T = 150 #horizon
p = plot() #naming plot to add
y0 = 1
τ = 5 # moving average - interval

for θ in Θ
    y = zeros(T+1)
    y[1] = y0
    for i ∈ 1:T
        y[i+1] = γ + θ * y[i] + σ * randn()
    end
    plot!(p, y, label = "θ = $θ")

    # moving average
    x = zeros(T+1)
    for i ∈ τ:(T+1)
        x[i]= sum(y[(i-τ+1):i])/τ
    end

    plot!(p, x, label = "rolling mean θ = $θ" )
end
p

# simulate "nsamp" times

function sim_randomwalk_drift(Θ,γ,σ,T)

    y_sim = zeros(nsamp,length(Θ))
    for j in 1:nsamp

        y = zeros(T+1,length(Θ))
        y[1,:] = ones(eltype(Θ), length(Θ))

        for i in 1:T
            y[i+1,:] = γ .+ Θ .*y[i,:] + σ.*randn(3)
        end

        y_sim[j,:] = y[T+1,:]
    end
    return y_sim
end

y_sim = sim_randomwalk_drift(Θ,γ,σ,T)

# find mean and variance of ensemble averages

@show (mean_θ1, mean_θ2, mean_θ3) = (mean(y_sim[:,1]),mean(y_sim[:,2]),mean(y_sim[:,3]))
@show (var_θ1, var_θ2, var_θ3) = (var(y_sim[:,1]),var(y_sim[:,2]),var(y_sim[:,3]))

# plot histogram of T+1 values

h=histogram()
for (i,θ) in enumerate(Θ)
    histogram!(h, y_sim[:,i],α =0.5, label = "θ = $θ",bins =20)
end
h

##############
# Exercise 3
##############

# simulate observations

N = 50 # number of observations

x_1, x_2 = (randn(N), randn(N))

a, b, c, d, σ = (0.1, 0.2, 0.5, 1.0, 0.1) #parameters

# simulate data

function simdata(a, b, c, d, σ, N; M = 20)
    y = zeros(N,M)
    for i in 1:N
        y[i,:] = a * x_1[i] + b * x_1[i]^2 + c * x_2[i] + d .+ σ.*randn(M)
    end
    return y
end

Y = simdata(a, b, c, d, σ, N)

# perform OLS

X = [x_1 x_1.^2 x_2 ones(N)]
β = zeros(M,size(X)[2])

for m in 1:M
    β[m,:] = inv(X'*X)*X'*Y[:,m]
end
β

# plot histograms
Θ = ["a=0.1" "b=0.2" "c=0.5" "d=1.0"]
h = histogram()

for (i,θ) in enumerate(Θ)
    histogram!(h, β[:,i],α = 0.8, label = θ, bins = 10)
end
h

############################
# Exercise 4 - use NL solve
############################

using NLsolve

A = [0.8  -0.2; -0.1  0.7]
Σ = [0.5 0.4; 0.4 0.6]
Q = Σ * Σ'
S0 = Σ * Σ'

f(S) = A * S * A' + Q

f(S0)
sol = fixedpoint(f,S0)
sol

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
# Variables

"""
    struct AddVariableNotAllowed <: NotAllowedError
        message::String # Human-friendly explanation why the attribute cannot be set
    end

An error indicating that variables cannot be added to the model.
"""
struct AddVariableNotAllowed <: NotAllowedError
    message::String # Human-friendly explanation why the attribute cannot be set
end
AddVariableNotAllowed() = AddVariableNotAllowed("")

operation_name(::AddVariableNotAllowed) = "Adding variables"

"""
    add_variables(model::ModelLike, n::Int)::Vector{VariableIndex}

Add `n` scalar variables to the model, returning a vector of variable indices.

A [`AddVariableNotAllowed`](@ref) error is thrown if adding variables cannot be
done in the current state of the model `model`.
"""
add_variables(model::ModelLike, n) = throw(AddVariableNotAllowed())

"""
    add_variable(model::ModelLike)::VariableIndex

Add a scalar variable to the model, returning a variable index.

A [`AddVariableNotAllowed`](@ref) error is thrown if adding variables cannot be
done in the current state of the model `model`.
"""
add_variable(model::ModelLike) = throw(AddVariableNotAllowed())

"""
    add_constrained_variable(
        model::ModelLike,
        set::AbstractScalarSet
    )::Tuple{MOI.VariableIndex,
             MOI.ConstraintIndex{MOI.SingleVariable, typeof(set)}}

Add to `model` a scalar variable constrained to belong to `set`, returning the
index of the variable created and the index of the constraint constraining the
variable to belong to `set`.

By default, this function falls back to creating a free variable with
[`add_variable`](@ref) and then constraining it to belong to `set` with
[`add_constraint`](@ref).
"""
function add_constrained_variable(model::ModelLike, set::AbstractScalarSet)
    variable = add_variable(model)
    constraint = add_constraint(model, SingleVariable(variable), set)
    return variable, constraint
end

"""
    add_constrained_variables(
        model::ModelLike,
        sets:AbstractVector{<:AbstractScalarSet}
    )::Tuple{Vector{MOI.VariableIndex},
             Vector{MOI.ConstraintIndex{MOI.SingleVariable, eltype(sets)}}}

Add to `model` scalar variables constrained to belong to `sets`, returning the
indices of the variables created and the indices of the constraints constraining
the variables to belong to each set in `sets`. That is, if it returns `variables`
and `constraints`, `constraints[i]` is the index of the constraint constraining
`variable[i]` to belong to `sets[i]`.

By default, this function falls back to calling
[`add_constrained_variable`](@ref) on each set.
"""
function add_constrained_variables(model::ModelLike, sets::AbstractVector{<:AbstractScalarSet})
    variables = Vector{VariableIndex}(undef, length(sets))
    constraints = Vector{ConstraintIndex{SingleVariable, eltype(sets)}}(undef, length(sets))
    for (i, set) in enumerate(sets)
        variables[i], constraints[i] = add_constrained_variable(model, set)
    end
    return variables, constraints
end

"""
    add_constrained_variables(
        model::ModelLike,
        set::AbstractVectorSet
    )::Tuple{Vector{MOI.VariableIndex},
             MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}}

Add to `model` a vector of variables constrained to belong to `set`, returning
the indices of the variables created and the index of the constraint
constraining the vector of variables to belong to `set`.

By default, this function falls back to creating free variables with
[`add_variables`](@ref) and then constraining it to belong to `set` with
[`add_constraint`](@ref).
"""
function add_constrained_variables(model::ModelLike, set::AbstractVectorSet)
    variables = add_variables(model, dimension(set))
    constraint = add_constraint(model, VectorOfVariables(variables), set)
    return variables, constraint
end


text:
using SchumakerSpline
tol = 10*eps()

x = [1,2,3,4,5,6,7,8,9,10,11,12]
y = log.(x) + sqrt.(x)

typeof(x) == Vector{Int}

spline = Schumaker(x,y)
for i in 1:length(x)
    abs(evaluate(spline, x[i]) - y[2]) < tol
end
# Evaluation with a Float64.
evaluate(spline, 11.5)

# Testing second derivatives
xArray = range(1, stop=6, length=1000)
second_derivatives = evaluate.(spline, xArray,2)
maximum(second_derivatives) < tol

# Testing Integrals
analytic_integral(lhs,rhs) = rhs*log(rhs) - rhs + (2/3) * rhs^(3/2) - ( lhs*log(lhs) - lhs + (2/3) * lhs^(3/2) )
lhs = 2.0
rhs = 2.5
numerical_integral = evaluate_integral(spline, lhs,rhs)
abs(analytic_integral(lhs,rhs) - numerical_integral) < 0.01

lhs = 2.1
rhs = 2.11
numerical_integral = evaluate_integral(spline, lhs,rhs)
abs(analytic_integral(lhs,rhs) - numerical_integral) < 0.01

lhs = 1
rhs = 4
numerical_integral = evaluate_integral(spline, lhs,rhs)
abs(analytic_integral(lhs,rhs) - numerical_integral) < 0.1


text:
using QuantumCumulants
using SymbolicUtils
using Test

@testset "parameters" begin

hf = FockSpace(:cavity)
ha = NLevelSpace(:atom,(:g,:e))
h = hf⊗ha

a = Destroy(h,:a)

@cnumbers p q
# @test p*q isa NumberTerm
@test p*q*a isa QuantumCumulants.QTerm

σ = Transition(h,:σ,:g,:e)
σee = Transition(h,:σ,:e,:e)

# JC
@cnumbers ωc ωa g
H = ωc*a'*a + ωa*σ'*σ + g*(a'*σ + σ'*a)

da = commutator(1im*H,a)
@test iszero(simplify(da - ((-1im*g)*σ + (-1im*ωc)*a)))
ds = commutator(1im*H,σ)
@test iszero(simplify(ds - ((-1im*g)*a + (-1im*ωa)*σ + (2im*g)*a*σee)))
dn = commutator(1im*H,a'*a)
@test iszero(simplify(dn - ((-1im*g)*a'*σ + (1im*g)*a*σ')))

he = meanfield([a,σ,a'*a],H)
@test iszero(simplify(he.operator_equations[1].rhs - da))
@test iszero(simplify(he.operator_equations[2].rhs - ds))
@test iszero(simplify(he.operator_equations[3].rhs - dn))

# Lossy JC
@cnumbers κ γ
J = [a,σ]
he_diss = meanfield([a,σ,σ'*σ],H,J;rates=[κ,γ])

@test iszero(simplify(he_diss.operator_equations[1].rhs - ((-1im*ωc - 0.5κ)*a + (-1im*g)*σ); expand=true))
@test iszero(simplify(he_diss.operator_equations[2].rhs - ((-1im*g)*a + (-1im*ωa - 0.5γ)*σ + (2im*g)*a*σee); expand=true))
@test iszero(simplify(he_diss.operator_equations[3].rhs - ((-γ)*σee + (1im*g)*a'*σ + (-1im*g)*a*σ')))

# Single-atom laser
@cnumbers ν
J = [a,σ,σ']
he_laser = meanfield([a'*a,σ'*σ,a'*σ],H,J;rates=[κ,γ,ν])

@test iszero(simplify(he_laser.operator_equations[1].rhs - ((-κ)*a'*a + (-1im*g)*a'*σ + (1im*g)*a*σ')))
@test iszero(simplify(he_laser.operator_equations[2].rhs - (ν + (-ν - γ)*σee + (1im*g)*a'*σ + (-1im*g)*a*σ'); expand=true))
# @test iszero(simplify(he_laser.operator_equations[3].rhs - ((1im*g)*σee + (-1im*g)*a'*a + (1im*ωc + -1im*ωa - 0.5*(κ + γ + ν))*a'*σ + (2im*g)*a'*a*σee); expand=true))

# Test sorting of longer term
@cnumbers λ Γ N
b = Destroy(h, :b)
yy = 1im*σ*b*(N-1)*λ*Γ
# @test isequal(simplify(yy), QuantumCumulants.QMul(im*Γ*λ*(N-1), [b, σ]))

# Test numbers in commutator
@cnumbers Δ
h = NLevelSpace(:a, 2)
s(i,j) = Transition(h, :s, i, j)
H = Δ*s(2,2) - Δ*s(1,1)
H = simplify(H)
@test iszero(simplify(meanfield(s(1,2), H).operator_equations[1].rhs  + 2im*Δ*s(1,2)))

end # testset


text:
using LinearAlgebra, Random, Plots; pyplot()
Random.seed!(1)

N = 10^4
edges = [(1,2), (1,3), (2,4), (1,4), (3,4)]
L = maximum(maximum.(edges))
source, dest = 1, L

function adjMatrix(edges,L)
    R = zeros(Int, L, L)
    for e in edges
        R[ e[1], e[2] ], R[ e[2], e[1] ] = 1, 1
    end
    R
end

pathExists(R, source, destination) = sign.((I+R)^L)[source,destination]
randNet(p) = randsubseq(edges,1-p)

relEst(p) = sum([pathExists(adjMatrix(randNet(p),L),source,dest) for _ in 1:N])/N
relAnalytic(p) = 1-p^3*(p-2)^2

pGrid = 0:0.05:1
scatter(pGrid, relEst.(pGrid),
	c=:blue, ms=5, msw=0,label="Monte Carlo")
plot!(pGrid, relAnalytic.(pGrid),
	c=:red, label="Analytic", xlims=(0,1.05), ylims=(0,1.05),
	xlabel="p", ylabel="Reliability")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
@testset "Parsing" begin

    @test isa(from_string("a", base=false), Constant)
    @test isa(from_string("X"), Variable)
    @test from_string("a").name == "a"
    @test isa(from_string("_"), Variable)
    @test isa(from_string("'UppercaseConstant'", base=false), Constant)

    s = from_string("t(a,b,c)", base=false)
    @test isa(s, Structure)
    @test s.functor.arity == 3
    @test all(isa(a, Constant) for a in s.arguments)

    l = from_string("pred(a,X)")
    @test isa(l, Literal)
    @test l.predicate.arity == 2
    @test isa(l.arguments[begin], Constant)
    @test isa(l.arguments[2], Variable)

    l2 = from_string("[1,2,3,4,5]")
    @test isa(l2, List)
    @test all([isa(a, Integer) for a in l2.elements])

    l3 = from_string("[H|T]")
    @test isa(l3, LPair)

    l4 = from_string("p([a,b,c], s(1,2))")
    @test isa(l4, Literal)
    @test isa(l4.arguments[begin], List)
    @test isa(l4.arguments[2], Structure)

    cl = from_string("h(X,Y) :- p(X,Z), q(Z,Y), r(Z)")
    @test isa(cl, Clause)

end

@testset "macro parsing" begin
    program = @prolog [
        a,
        p(a, b),
        p(X),
        q(X,Y),
        p(1),
        r(a, s(b,c)),
        r(c, s(c,X)),
        p(X) <= q(X,Y) & z(Y),
        p2(X) <= p(X) & !z(X),
        p3(X|Y, Z) <= r(X,Z) & z(Y)
    ]

    clauses = program.clauses

    @test isa(clauses[begin], Proposition)
    @test isa(clauses[2], Literal)
    @test is_ground(clauses[2])
    @test isa(clauses[3], Literal)
    @test isa(clauses[4], Literal)
    @test isa(clauses[5], Literal)
    @test isa(clauses[6], Literal)
    @test isa(clauses[7], Literal)
    @test isa(clauses[8], Clause)
    @test isa(clauses[9], Clause)
    @test isa(clauses[10], Clause)

end

@testset "problog macro parsing" begin
    @test isa((@problog 0.5::p | 0.5::q), AnnotatedDisjunction)
    @test isa((@problog 0.5::p(X) | 0.5::q(X)), AnnotatedDisjunction)
    @test isa((@problog 0.3::p(X,a) | 0.3::p(X,b) | 0.3::p(X,c)), AnnotatedDisjunction)
end


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

@doc raw"""LoadBalancerStatus represents the status of a load-balancer.

    IoK8sApiCoreV1LoadBalancerStatus(;
        ingress=nothing,
    )

    - ingress::Vector{IoK8sApiCoreV1LoadBalancerIngress} : Ingress is a list containing ingress points for the load-balancer. Traffic intended for the service should be sent to these ingress points.
"""
mutable struct IoK8sApiCoreV1LoadBalancerStatus <: SwaggerModel
    ingress::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1LoadBalancerIngress} } # spec name: ingress

    function IoK8sApiCoreV1LoadBalancerStatus(;ingress=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1LoadBalancerStatus, Symbol("ingress"), ingress)
        setfield!(o, Symbol("ingress"), ingress)
        o
    end
end # type IoK8sApiCoreV1LoadBalancerStatus

const _property_map_IoK8sApiCoreV1LoadBalancerStatus = Dict{Symbol,Symbol}(Symbol("ingress")=>Symbol("ingress"))
const _property_types_IoK8sApiCoreV1LoadBalancerStatus = Dict{Symbol,String}(Symbol("ingress")=>"Vector{IoK8sApiCoreV1LoadBalancerIngress}")
Base.propertynames(::Type{ IoK8sApiCoreV1LoadBalancerStatus }) = collect(keys(_property_map_IoK8sApiCoreV1LoadBalancerStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1LoadBalancerStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1LoadBalancerStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1LoadBalancerStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1LoadBalancerStatus[property_name]

function check_required(o::IoK8sApiCoreV1LoadBalancerStatus)
    true
end

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


text:
using Printf

"""
Represents the time of day.

fields: hour, minute, second
"""
struct MyTime
    hour
    minute
    second
end

# Exercise 16-1
function printtime(t::MyTime)
    @printf("%02d:%02d:%02d", t.hour, t.minute, t.second)
end

t = MyTime(23, 10, 58)
# printtime(t)

# Exercise 16-2
function isafter(t1::MyTime, t2::MyTime)
    return (t1.hour, t1.minute, t1.second) > (t2.hour, t2.minute, t2.second)
end

# Exercise 16-3
function increment!(time, seconds)
    time.second += seconds
    dmin, time.second = divrem(time.second, 60)
    time.minute += dmin
    dhour, time.minute = divrem(time.minute, 60)
    time.hour += dhour
    time.hour %= 24
    nothing
end

# Exercise 16-4
function increment(time, seconds)
    nsecond = time.second + seconds
    dmin, nsecond = divrem(nsecond, 60)
    nminute = time.minute + dmin
    dhour, nminute = divrem(nminute, 60)
    nhour = time.hour + dhour
    nhour %= 24
    return MyTime(nhour, nminute, nsecond)
end

# Exercise 16-5
function timetoint(time)
    minutes = time.hour * 60 + time.minute
    seconds = minutes * 60 + time.second
end

function inttotime(seconds)
    (minutes, second) = divrem(seconds, 60)
    hour, minute = divrem(minutes, 60)
    MyTime(hour, minute, second)
end

function increment(time, seconds)
    return inttotime(timetoint(time) + seconds)
end

# Exercise 16-6
function multime(time::MyTime, x::Number)
    return inttotime(timetoint(time) * x)
end

# println(multime(t, 3))

function timepermile(time::MyTime, dist::Number)
    return multime(time, 1/dist)
end

# println(timepermile(t, 3))

# Exercise 16-7
using Dates

td = today()
println(td)
println(dayofweek(td))
println(dayname(td))

function age(bday)
    td = today()
    td - bday
    dyear = year(td) - year(bday)
    dmonth = month(td) - month(bday)
    dday = day(td) - day(bday)
    if dmonth < 0
        dyear -= 1
    elseif dmonth == 0 && dday < 0
        dyear -= 1
    end
    println("Age: $dyear")

    next_bday = bday + Year(dyear + 1)
    println("Next birthday: $next_bday")
    nextdiff = DateTime(next_bday) - now()
    timeleft = canonicalize(Dates.CompoundPeriod(nextdiff))
    println("Time to next birthday: $timeleft")
end

# age(Date(1995, 11, 02))

function doubleday(d1, d2)
    @assert(d2 > d1, "d2 should be later than d1")
    return d2 - d1 + d2
end

d1 = Date(2000, 1, 2)
d2 = Date(2005, 6, 24)
# println(doubleday(d1, d2))

function n_day(d1, d2, n)
    @assert(d2 > d1, "d2 should be later than d1")
    diff = d2 - d1
    age1, age2 = diff.value, 0
    while age2 * n != age1
        if age2 * n > age1
            error("Impossible combination of inputs")
        end
        age1 += 1
        age2 += 1
    end
    return d1 + Day(age1)
end

# println(n_day(d1, d2, 2))
# println(n_day(d1, d2, 3))
# println(n_day(d1, d2, 6))


text:
using ExponentialSmoothing
using Test


text:
# disk-phantom.jl

using ImagePhantoms: disk_phantom_params, Ellipse, phantom, Object2d
#using MIRTjim: jim
using Test: @test, @testset, @test_throws, @inferred

@testset "disk" begin
    params = @inferred disk_phantom_params( ; rhead = () -> rand(100:105))
    @test params isa Matrix{Float32}

    ob = Ellipse(params)
    @test ob isa Vector{<:Object2d}
#=
    x = LinRange(-1,1,200)*130
    y = LinRange(-1,1,202)*130
    tmp = phantom(x, y, ob)
    jim(tmp)
=#
end


text:
@info "Testing PointCloud..."
@testset "PointCloud" begin
    for T in [Float32, Float64], _N in [true, false]
        points = rand(T, 3, 16, 4)
        _normals = rand(T, 3, 16, 4)
        if _N
            normals = _normals
        else
            normals = nothing
        end
        rep = PointCloud(points, normals)
        if normals == nothing
            @test rep.normals isa Nothing
        else
            @test size(rep.normals) == size(_normals)
            @test rep.normals isa Array{Float32,3}
        end
        @test npoints(rep) == size(points, 2)
        @test size(rep.points) == size(points)
        @test rep.points isa Array{Float32,3}

        for i = 1:size(points, 3)
            @test all(isapprox.(rep[i], points[:, :, i]))
        end

        # using other contructor
        rep1 = PointCloud(points = points, normals = normals)
        rep2 = PointCloud(rep)
        @test rep.points == rep1.points
        @test rep.points == rep2.points
        @test rep.normals == rep1.normals
        @test rep.normals == rep2.normals
    end
end

@info "Testing VoxelGrid..."
@testset "VoxelGrid" begin
    res = 4
    voxels = rand(Float32, res, res, res, 2)
    v = VoxelGrid(voxels)
    @test Flux3D._assert_voxel(v)
    @test v isa VoxelGrid{Float32}
    @test v[1] == voxels[:, :, :, 1]
    @test v[2] == voxels[:, :, :, 2]
    @test v.voxels isa Array{Float32,4}
    @test size(v.voxels) == (res, res, res, 2)

    voxels = rand(Float32, res, res, res)
    v2 = VoxelGrid(voxels)
    @test v2 isa VoxelGrid{Float32}
    @test v2[1] == voxels[:, :, :, 1]
    @test size(v2.voxels) == (res, res, res, 1)
end

@info "Testing TriMesh..."
@testset "TriMesh" begin

    T, R = (Float32, Int64)
    verts1 = [
        0.1 0.3 0.5
        0.5 0.2 0.1
        0.6 0.8 0.7
    ]
    verts2 = [
        0.1 0.3 0.3
        0.6 0.7 0.8
        0.2 0.3 0.4
        0.1 0.5 0.3
    ]
    verts3 = [
        0.7 0.3 0.6
        0.2 0.4 0.8
        0.9 0.5 0.2
        0.2 0.3 0.4
        0.9 0.3 0.8
    ]

    verts_list = [T.(verts1'), T.(verts2'), T.(verts3')]

    faces1 = [1 2 3]
    faces2 = [
        1 2 3
        2 3 4
    ]
    faces3 = [
        2 3 1
        1 2 4
        3 4 2
        5 4 3
        5 1 2
        5 4 2
        5 3 2
    ]

    faces_list = [R.(faces1'), R.(faces2'), R.(faces3')]
    m = TriMesh(verts_list, faces_list)

    # IO Tests
    @testset "IO" begin
        mktempdir() do tmpdir
            # FIXME: MeshIO 2dm save/load breaking #62, #63
            # for ext in ["obj", "off", "stl", "ply", "2dm"]
            for ext in ["obj", "off", "stl", "ply"]
                save_trimesh(joinpath(tmpdir, "test.$(ext)"), m)
                m_loaded = load_trimesh(joinpath(tmpdir, "test.$(ext)"))
                @test all(isapprox.(get_verts_packed(m_loaded), verts_list[1]))
                @test get_faces_packed(m_loaded) == faces_list[1]
            end
        end
    end

    @test all(verts_list .== get_verts_list(m))
    @test cat(verts_list...; dims = 2) == get_verts_packed(m)
    _padded = get_verts_padded(m)
    for (i, v) in enumerate(verts_list)
        @test _padded[:, 1:size(v, 2), i] == v
        @test all(_padded[:, size(v, 2)+1:end, i] .== 0)
    end

    @test all(faces_list .== get_faces_list(m))
    _packed = get_faces_packed(m)
    _cur_idx = 1
    _offset = 0
    for (i, f) in enumerate(faces_list)
        @test _packed[:, _cur_idx:_cur_idx+size(f, 2)-1] == (f .+ _offset)
        _cur_idx += size(f, 2)
        _offset += size(verts_list[i], 2)
    end
    _padded = get_faces_padded(m)
    for (i, v) in enumerate(faces_list)
        @test _padded[:, 1:size(v, 2), i] == v
        @test all(_padded[:, size(v, 2)+1:end, i] .== 0)
    end

    _packed = get_faces_packed(m)
    e12 = cat(_packed[1, :], _packed[2, :], dims = 2)
    e23 = cat(_packed[2, :], _packed[3, :], dims = 2)
    e31 = cat(_packed[3, :], _packed[1, :], dims = 2)
    edges = cat(e12, e23, e31, dims = 1)
    edges = sort(edges, dims = 2)
    edges = sortslices(edges; dims = 1)
    edges = unique(edges; dims = 1)
    @test edges == get_edges_packed(m)

    _dict = get_edges_to_key(m)
    for i = 1:size(edges, 1)
        @test i == _dict[Tuple(edges[i, :])]
    end

    _f_edges = get_faces_to_edges_packed(m)
    _packed = get_faces_packed(m)
    for i = 1:size(_f_edges, 1)
        @test edges[_f_edges[i, 1], :] == sort(_packed[[2, 3], i])
        @test edges[_f_edges[i, 2], :] == sort(_packed[[1, 3], i])
        @test edges[_f_edges[i, 3], :] == sort(_packed[[1, 2], i])
    end

    _edges = get_edges_packed(m)
    V = size(get_verts_packed(m), 2)
    L = zeros(V, V)
    deg = zeros(V, V)
    for i = 1:size(_edges, 1)
        L[_edges[i, 1], _edges[i, 2]] = 1
        L[_edges[i, 2], _edges[i, 1]] = 1
    end
    deg = sum(L; dims = 2)
    deg = map(x -> (x > 0 ? 1 / x : x), deg)
    for i = 1:V, j = 1:V
        if i == j
            L[i, j] = -1
        elseif L[i, j] == 1
            L[i, j] = deg[i]
        end
    end
    @test all(isapprox.(L, Array(get_laplacian_packed(m)), rtol = 1e-5, atol = 1e-5))

    # verts_normals, faces_normals, faces_areas
    verts1 = [
        0.1 0.3 0.0
        0.5 0.2 0.0
        0.6 0.8 0.0
        0.0 0.3 0.2
        0.0 0.2 0.5
        0.0 0.8 0.7
        0.5 0.0 0.2
        0.6 0.0 0.5
        0.8 0.0 0.7
        0.0 0.0 0.0
        0.0 0.0 0.0
        0.0 0.0 0.0
    ]

    verts2 = [
        0.1 0.3 0.0
        0.5 0.2 0.0
        0.6 0.8 0.0
        0.0 0.3 0.2
        0.0 0.2 0.5
        0.0 0.8 0.7
    ]

    verts2 = [
        0.1 0.3 0.0
        0.5 0.2 0.0
        0.0 0.3 0.2
        0.0 0.2 0.5
        0.0 0.8 0.7
    ]

    faces1 = [
        1 2 3
        4 5 6
        7 8 9
        10 11 12
    ]

    faces2 = [
        1 2 3
        3 4 5
    ]

    _normal1 = [
        0.0 0.0 1.0
        0.0 0.0 1.0
        0.0 0.0 1.0
        -1.0 0.0 0.0
        -1.0 0.0 0.0
        -1.0 0.0 0.0
        0.0 1.0 0.0
        0.0 1.0 0.0
        0.0 1.0 0.0
        0.0 0.0 0.0
        0.0 0.0 0.0
        0.0 0.0 0.0
    ]

    _normal2 = [
        -0.2408 -0.9631 -0.1204
        -0.2408 -0.9631 -0.1204
        -0.9389 -0.3414 -0.0427
        -1.0000 0.0000 0.0000
        -1.0000 0.0000 0.0000
    ]

    _fnormal1 = [
        -0.0 0.0 1.0
        -1.0 0.0 0.0
        0.0 1.0 0.0
        0.0 0.0 0.0
    ]

    _fnormal2 = [
        -0.2408 -0.9631 -0.1204
        -1.0000 0.0000 0.0000
    ]

    _farea1 = [0.125, 0.1, 0.02, 0.0]
    _farea2 = [0.0415, 0.1]

    _normal1 = _normal1'
    _normal2 = _normal2'
    _fnormal1 = _fnormal1'
    _fnormal2 = _fnormal2'
    _farea1 = _farea1'
    _farea2 = _farea2'

    # verts_normals, faces_normals, faces_areas
    for T in [Float64, Float32], R in [Int64, UInt32]
        verts = [T.(verts1'), T.(verts2')]
        faces = [R.(faces1'), R.(faces2')]
        m = TriMesh(verts, faces)

        @test compute_verts_normals_packed(m) isa Array{T,2}
        @test compute_verts_normals_padded(m) isa Array{T,3}
        @test compute_verts_normals_list(m) isa Array{Array{T,2},1}
        @test all(
            isapprox.(
                cat(_normal1, _normal2; dims = 2),
                compute_verts_normals_packed(m),
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                [_normal1, _normal2],
                compute_verts_normals_list(m),
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                _normal1,
                compute_verts_normals_padded(m)[:, 1:size(_normal1, 2), 1],
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                _normal2,
                compute_verts_normals_padded(m)[:, 1:size(_normal2, 2), 2],
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(compute_verts_normals_padded(m)[:, size(_normal1, 2)+1:end, 1] .== 0)
        @test all(compute_verts_normals_padded(m)[:, size(_normal2, 2)+1:end, 2] .== 0)

        @test compute_faces_normals_packed(m) isa Array{T,2}
        @test compute_faces_normals_padded(m) isa Array{T,3}
        @test compute_faces_normals_list(m) isa Array{Array{T,2},1}
        @test all(
            isapprox.(
                cat(_fnormal1, _fnormal2; dims = 2),
                compute_faces_normals_packed(m),
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                [_fnormal1, _fnormal2],
                compute_faces_normals_list(m),
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                _fnormal1,
                compute_faces_normals_padded(m)[:, 1:size(_fnormal1, 2), 1],
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                _fnormal2,
                compute_faces_normals_padded(m)[:, 1:size(_fnormal2, 2), 2],
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(compute_faces_normals_padded(m)[:, size(_fnormal1, 2)+1:end, 1] .== 0)
        @test all(compute_faces_normals_padded(m)[:, size(_fnormal2, 2)+1:end, 2] .== 0)

        @test compute_faces_areas_packed(m) isa Array{T,1}
        @test compute_faces_areas_padded(m) isa Array{T,3}
        @test compute_faces_areas_list(m) isa Array{Array{T,2},1}
        @test all(
            isapprox.(
                reshape(cat(_farea1, _farea2; dims = 2), :),
                compute_faces_areas_packed(m),
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                [_farea1, _farea2],
                compute_faces_areas_list(m),
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                _farea1,
                compute_faces_areas_padded(m)[:, 1:size(_farea1, 2), 1],
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(
            isapprox.(
                _farea2,
                compute_faces_areas_padded(m)[:, 1:size(_farea2, 2), 2],
                rtol = 1e-4,
                atol = 1e-4,
            ),
        )
        @test all(compute_faces_areas_padded(m)[:, size(_farea1, 2)+1:end, 1] .== 0)
        @test all(compute_faces_areas_padded(m)[:, size(_farea2, 2)+1:end, 2] .== 0)
    end
end

@info "Testing Rep Utils..."
@testset "Rep Utils" begin
    for T in [Float64, Float32, Int64, UInt32]
        a1 = [
            1 2 3
            4 5 6
            7 8 9
            10 11 12
        ]

        a2 = [
            13 14 15
            16 17 18
        ]

        a3 = [
            19 20 21
            22 23 24
            25 26 27
        ]

        a1 = T.(a1')
        a2 = T.(a2')
        a3 = T.(a3')
        _list = [a1, a2, a3]

        _packed = [
            1 2 3
            4 5 6
            7 8 9
            10 11 12
            13 14 15
            16 17 18
            19 20 21
            22 23 24
            25 26 27
        ]
        _packed = T.(_packed')

        _padded = zeros(3, 4, 3)
        _padded[:, 1:size(a1, 2), 1] = a1
        _padded[:, 1:size(a2, 2), 2] = a2
        _padded[:, 1:size(a3, 2), 3] = a3
        _padded = T.(_padded)

        items_len, packed_first_idx, packed_to_list_idx = Flux3D._auxiliary_mesh(_list)

        @test items_len == [4, 2, 3]
        @test packed_first_idx == [1, 5, 7]
        @test all(packed_to_list_idx[1:4] .== 1)
        @test all(packed_to_list_idx[5:6] .== 2)
        @test all(packed_to_list_idx[7:9] .== 3)

        @test Flux3D._list_to_padded(_list, 0) == _padded
        @test Flux3D._list_to_padded(_list, 0) isa Array{T,3}
        @test gradient(x -> 0.5 .* sum(Flux3D._list_to_padded(x, 0) .^ 2), _list)[1] ==
              _list

        @test Flux3D._list_to_packed(_list) == _packed
        @test Flux3D._list_to_packed(_list) isa Array{T,2}
        @test gradient(x -> 0.5 .* sum(Flux3D._list_to_packed(x) .^ 2), _list)[1] ==
              Tuple(_list)

        @test Flux3D._packed_to_padded(_packed, items_len, 0) == _padded
        @test Flux3D._packed_to_padded(_packed, items_len, 0) isa Array{T,3}
        @test gradient(
            x -> 0.5 .* sum(Flux3D._packed_to_padded(x, items_len, 0) .^ 2),
            _packed,
        )[1] == _packed

        @test Flux3D._packed_to_list(_packed, items_len) == _list
        @test Flux3D._packed_to_list(_packed, items_len) isa Array{Array{T,2},1}
        gs = gradient(_packed) do x
            loss = 0
            for x in Flux3D._packed_to_list(x, items_len)
                loss += sum(x .^ 2)
            end
            return 0.5 * loss
        end
        @test gs[1] == _packed

        @test Flux3D._padded_to_list(_padded, items_len) == _list
        @test Flux3D._padded_to_list(_padded, items_len) isa Array{Array{T,2},1}
        gs = gradient(_padded) do x
            loss = 0
            for x in Flux3D._padded_to_list(x, items_len)
                loss += sum(x .^ 2)
            end
            return 0.5 * loss
        end
        @test gs[1] == _padded

        @test Flux3D._padded_to_packed(_padded, items_len) == _packed
        @test Flux3D._padded_to_packed(_padded, items_len) isa Array{T,2}
        @test gradient(
            x -> sum(0.5 .* sum(Flux3D._padded_to_packed(x, items_len) .^ 2)),
            _padded,
        )[1] == _padded
    end
end


text:
struct OptFunc{M<:DescentMethod,F<:Base.Callable}<:Function
    method::M
    fdf::F
end

(optfn::OptFunc)(x, α, d) = callfn!(optfn.method, optfn.fdf, x, α, d)

"""
    optimize!(M::Wrapper, fdf, x0)

Find an optimizer for `fdf`, starting with the initial approximation `x0`.
`fdf(x, g)` must return a tuple (f(x), ∇f(x)) and, if `g` is mutable, overwrite
it with the gradient.
"""
function optimize!(M::Wrapper, fdf, x0; reset=true, log_stream=nothing)
    return __optim_impl!(M, fdf, x0, log_stream; reset)
end

function __optim_impl!(M::Wrapper, fdf, x0, log_stream::IO; reset=true)
    logged_opt = Logger(M, log_stream)
    optfn! = OptFunc(logged_opt, fdf)

    init!(logged_opt, optfn!, x0; reset=reset)
    while !stopcond(logged_opt)
        step!(logged_opt, optfn!)
    end
    return convstat(logged_opt)
end

function __optim_impl!(M::Wrapper, fdf, x0, log_stream::AbstractString; reset=true)
    open(log_stream, "w") do log_io
        return __optim_impl!(M, fdf, x0, log_io; reset)
    end
end

function __optim_impl!(M::Wrapper, fdf, x0, ::Nothing; reset=true)
    optfn! = OptFunc(M, fdf)

    init!(M, optfn!, x0; reset=reset)
    while !stopcond(M)
        step!(M, optfn!)
    end
    return convstat(M)
end

function __optim_impl!(M::Logger, fdf, x0, ::Nothing; reset=true)
    optfn! = OptFunc(M, fdf)

    init!(M, optfn!, x0; reset=reset)
    while !stopcond(M)
        step!(M, optfn!)
    end
    return convstat(M)
end

function __optim_impl!(M::Logger, fdf, x0, log_stream; reset)
    throw(ArgumentError("Cannot wrap a logger into another logger"))
end
"""
    optimize!(
        M::CoreMethod, fdf, x0;
        gtol=1e-6,
        maxiter=100,
        maxcalls=nothing,
        reset=true,
        constrain_step=nothing,
        log_stream=nothing
    )

Find an optimizer for `fdf`, starting with the initial approximation `x0`.

# Arguments:
- `M::CoreMethod`: the core method to use for optimization
- `fdf(x, g)::Function`: function to optimize. It must return a tuple (f(x), ∇f(x)) and,
    if `g` is mutable, overwrite
    it with the gradient.
- `x0`: initial approximation

# Keywords:
- `gtol::Real`: stop optimization when the gradient norm is less
- `maxiter::Integer`: force stop optimization after this number of iterations
    (use `nothing` or a negative value to not constrain iteration number)
- `maxcalls::Integer`: force stop optimization after this number of function calls
    (use `nothing` or a negative value to not constrain call number)
- `reset=true`: a value to pass as a keyword argument to the optimizer `init!` method
- `constrain_step(x0, d)`: a function to constrain step from `x0` in the direction `d`.
    It must return a real-numbered value `α` such that `x0 + αd` is the maximum allowed step
- `log_stream::Union{IO,AbstractString}`: IO stream or a file name to log the optimization
    process
"""
function optimize!(
    M::CoreMethod, fdf, x0;
    gtol = convert(float(eltype(x0)), 1e-6),
    maxiter = 100,
    maxcalls = nothing,
    reset = true,
    constrain_step = nothing,
    log_stream = nothing
)
    if !isnothing(gtol) && gtol > 0
        M = StopByGradient(M, gtol)
    end
    if isnothing(maxiter) || maxiter < 0
        M = LimitIters(M)
    else
        M = LimitIters(M, maxiter)
    end
    if isnothing(maxcalls) || maxcalls < 0
        M = LimitCalls(M)
    else
        M = LimitCalls(M, maxcalls)
    end
    if !isnothing(constrain_step)
        M = ConstrainStepSize(M, constrain_step)
    end
    optimize!(M, fdf, x0; reset=reset, log_stream=log_stream)
end

"""
    DescentMethods.solver(
        M::CoreMethod;
        gtol = 1e-6,
        maxiter = 100,
        maxcalls = nothing,
        constrain_step=nothing,
        log_stream=nothing
    )

Return the wrapper object for a chosen method to solve an optimization problem with given
    parameters. For the description of keywords, see [`optimize!`](@ref)
"""
function solver(
    M::CoreMethod;
    gtol=convert(eltype(argumentvec(M)), 1e-6),
    maxiter=100,
    maxcalls=nothing,
    constrain_step=nothing,
    log_stream=nothing
)
    if !isnothing(gtol) && gtol > 0
        M = StopByGradient(M, gtol)
    end
    if isnothing(maxiter) || maxiter < 0
        M = LimitIters(M)
    else
        M = LimitIters(M, maxiter)
    end
    if isnothing(maxcalls) || maxcalls < 0
        M = LimitCalls(M)
    else
        M = LimitCalls(M, maxcalls)
    end
    if !isnothing(constrain_step)
        M = ConstrainStepSize(M, constrain_step)
    end
    if !isnothing(log_stream)
        M = Logger(M, log_stream)
    end
    return M
end

"""
    DescentMethods.solver(M::DataType, x; gtol = convert(eltype(x0), 1e-6), maxiter = 100, maxcalls = nothing, constrain_step)

Return the wrapper object for a chosen method to solve an optimization problem with given
    parameters compatible with the dimensions of `x`.
"""
function solver(M::Type{<:CoreMethod}, x::AbstractVector; kw...)
    return solver(M(x); kw...)
end


text:
const N = 256
const N2 = N * div(N - 1, 2)
const step = 0.01
const xval_table = zeros(Float64, N)
const tsin_table = zeros(Float64, N)
const tcos_table = zeros(Float64, N)
const ttan_table = zeros(Float64, N)
const rsin_cache = Dict{Float64, Float64}()
const rcos_cache = Dict{Float64, Float64}()
const rtan_cache = Dict{Float64, Float64}()

function rho(x, y, rhocache, i, n)
    if n < 0
        return 0.0
    elseif n == 0
        return y[i+1]
    end
    idx = (N - 1 - n) * div(N - n, 2) + i
    if !haskey(rhocache, idx)
        rhocache[idx] = (x[i+1] - x[i + n+1]) / (rho(x, y, rhocache, i, n - 1) -
            rho(x, y, rhocache, i + 1, n - 1)) + rho(x, y, rhocache, i + 1, n - 2)
    end
    rhocache[idx]
end

function thiele(x, y, r, xin, n)
    if n > N - 1
        return 1.0
    end
    rho(x, y, r, 0, n) - rho(x, y, r, 0, n - 2) + (xin - x[n + 1]) / thiele(x, y, r, xin, n + 1)
end

function thiele_tables()
    for i in 1:N
        xval_table[i] = (i-1) * step
        tsin_table[i] = sin(xval_table[i])
        tcos_table[i] = cos(xval_table[i])
        ttan_table[i] = tsin_table[i] / tcos_table[i]
    end
    println(6 * thiele(tsin_table, xval_table, rsin_cache, 0.5, 0))
    println(3 * thiele(tcos_table, xval_table, rcos_cache, 0.5, 0))
    println(4 * thiele(ttan_table, xval_table, rtan_cache, 1.0, 0))
end

thiele_tables()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export normalize

function normalize(data,method)
  if method==:gauss
    μ=mean(data,2)
    σ=std(data,2)
  elseif method==:center
    ⇧=maximum(data,2);⇩=minimum(data,2)
    μ=(⇧+⇩)/2;σ=(⇧-⇩)/2
  elseif method==:positive
    ⇧=maximum(data,2);⇩=minimum(data,2)
    μ=⇩;σ=⇧-⇩
  end
  return (data.-μ)./σ
end


text:
export PolynomialHomotopyConfig, config

"""
    config(H::AbstractHomotopy, [x])

Create a homotopy corresponding to the homotopy type.
For a `AbstractPolynomialHomotopy` this is just `PolynomialHomotopyConfig`.
"""
config(H::AbstractPolynomialHomotopy) = PolynomialHomotopyConfig(H)
config(H::AbstractPolynomialHomotopy, x) = PolynomialHomotopyConfig(H, x)

"""
    PolynomialHomotopyConfig(H::AbstractPolynomialHomotopy{T}, [x::AbstractVector{S}])

A data structure with which `H` and it's derivatives can be evaluated efficiently.
Note that `x` is only used to determine the
output type of `H(x)`.

    PolynomialHomotopyConfig(H::AbstractPolynomialHomotopy{T}, [S])

Instead of a vector `x` a type can also be given directly.
"""
mutable struct PolynomialHomotopyConfig{T} <: AbstractHomotopyConfig{T}
    start::FP.JacobianConfig{T}
    target::FP.JacobianConfig{T}

    result_start::FP.JacobianDiffResult{T, Vector{T}, Matrix{T}}
    result_target::FP.JacobianDiffResult{T, Vector{T}, Matrix{T}}
end

function PolynomialHomotopyConfig(H::AbstractPolynomialHomotopy{T}) where T
    start = FP.JacobianConfig(H.start)
    target = FP.JacobianConfig(H.target)
    result_start = FP.JacobianDiffResult(start)
    result_target = FP.JacobianDiffResult(target)
    PolynomialHomotopyConfig(start, target, result_start, result_target)
end

function PolynomialHomotopyConfig(H::AbstractPolynomialHomotopy, ::AbstractVector{T}) where T
    start = FP.JacobianConfig(H.start, T)
    target = FP.JacobianConfig(H.target, T)
    result_start = FP.JacobianDiffResult(start)
    result_target = FP.JacobianDiffResult(target)
    PolynomialHomotopyConfig(start, target, result_start, result_target)
end

function evaluate_start_target!(
    cfg::PolynomialHomotopyConfig,
    H::AbstractPolynomialHomotopy,
    x::Vector,
    precomputed=false)

    FP.evaluate!(cfg.result_start.value, H.start, x, cfg.start, precomputed)
    FP.evaluate!(cfg.result_target.value, H.target, x, cfg.target, precomputed)
end

function jacobian_start_target!(
    cfg::PolynomialHomotopyConfig,
    H::AbstractPolynomialHomotopy,
    x::Vector,
    precomputed=false)

    FP.jacobian!(cfg.result_start.jacobian, H.start, x, cfg.start, precomputed)
    FP.jacobian!(cfg.result_target.jacobian, H.target, x, cfg.target, precomputed)
end

function evaluate_and_jacobian_start_target!(
    cfg::PolynomialHomotopyConfig,
    H::AbstractPolynomialHomotopy,
    x::Vector)
    FP.jacobian!(cfg.result_start, H.start, x, cfg.start)
    FP.jacobian!(cfg.result_target, H.target, x, cfg.target)
end

value_start(cfg::PolynomialHomotopyConfig) = cfg.result_start.value
value_target(cfg::PolynomialHomotopyConfig) = cfg.result_target.value

jacobian_start(cfg::PolynomialHomotopyConfig) = cfg.result_start.jacobian
jacobian_target(cfg::PolynomialHomotopyConfig) = cfg.result_target.jacobian


text:
let TDetection = AbstractObjectDetection{T} where T <: Real
    abstract type AbstractObjectDetector{TDetection} end

    detect(detector::AbstractObjectDetector{TDetection}, image)::Vector{TDetection} = error("No implementation for $(typeof(detector))")
end


text:
module TestExamples1DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "MHD Multicomponent" begin

  @trixi_testset "elixir_mhdmulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [8.16048158e-02, 5.46791194e-02, 5.46791194e-02, 1.54509265e-01, 4.13046273e-17, 5.47637521e-02, 5.47637521e-02, 8.37156486e-03, 1.67431297e-02, 3.34862594e-02],
      linf = [1.81982581e-01, 9.13611439e-02, 9.13611439e-02, 4.23831370e-01, 1.11022302e-16, 9.93731761e-02, 9.93731761e-02, 1.57164285e-02, 3.14328569e-02, 6.28657139e-02])
  end

  @trixi_testset "elixir_mhdmulti_ec.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [0.08153372259925547, 0.05464109003345891, 0.05464109003345891, 0.1540576724164453, 4.130462730494e-17, 0.054734930802131036, 0.054734930802131036, 0.008391254781284321, 0.016782509562568642, 0.033565019125137284],
      linf = [0.17492544007323832, 0.09029632168248182, 0.09029632168248182, 0.40798609353896564, 1.1102230246251565e-16, 0.09872923637833075, 0.09872923637833075, 0.01609818847160674, 0.03219637694321348, 0.06439275388642696],
      volume_flux = flux_derigs_etal)
  end

  @trixi_testset "elixir_mhdmulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
      l2   = [7.96871667e-02, 5.39815611e-02, 5.39815611e-02, 1.50151690e-01, 4.13046273e-17, 5.36296787e-02, 5.36296787e-02, 8.27911910e-03, 1.65582382e-02, 3.31164764e-02],
      linf = [1.41195919e-01, 7.82015593e-02, 7.82015593e-02, 3.39107175e-01, 1.11022302e-16, 7.00030105e-02, 7.00030105e-02, 1.49425583e-02, 2.98851165e-02, 5.97702330e-02])
  end

  @trixi_testset "elixir_mhdmulti_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_convergence.jl"),
      l2   = [1.84042774e-05, 3.38781451e-04, 3.38781451e-04, 6.70988781e-05, 4.13046273e-17, 3.34549620e-04, 3.34549620e-04, 2.64858640e-05, 5.29717279e-05, 1.05943456e-04],
      linf = [3.94780688e-05, 1.24292799e-03, 1.24292799e-03, 1.78659026e-04, 1.11022302e-16, 1.24361621e-03, 1.24361621e-03, 5.91550778e-05, 1.18310156e-04, 2.36620311e-04])
    end

  @trixi_testset "elixir_mhdmulti_briowu_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_briowu_shock_tube.jl"),
      l2   = [0.193310260592731, 0.35678108870884484, 0.0, 0.3635484869991876, 5.471141986014606e-16, 0.3668917193096584, 0.0, 0.056532537881382225, 0.11306507576276445],
      linf = [0.43382690098048815, 1.1031774119568432, 0.0, 1.049668973035275, 2.9976021664879227e-15, 1.5126193899444629, 0.0, 0.19658354010817752, 0.39316708021635505])
    end

end

end # module


text:
filename() = "day16.txt"

function read_data()
    data = String(read(filename()))
    blocks = split(data, "\n\n")

    function parse_rule(line)
        name, x = split(line, ": ")
        r1, r2 = split(x, " or ")
        r1 = range(parse.(Int, split(r1, "-"))...; step = 1)
        r2 = range(parse.(Int, split(r2, "-"))...; step = 1)
        return (name = name, ranges = [r1, r2])
    end

    function parse_ticket(line)
        return parse.(Int, split(line, ","))
    end

    return (
        rules = parse_rule.(split(blocks[1], "\n")),
        my_ticket = parse_ticket(split(blocks[2], "\n")[2]),
        nearby_tickets = parse_ticket.(split(blocks[3], "\n")[2:end])
    )
end

# ----------------------------------------------------------------------
# Part 1
# ----------------------------------------------------------------------

"""
Returns `true` if the value does not appear in any of the rules
"""
function error_val(rules, value)
    return !any(in.(value, rules))
end

"""
Return the sum of bad field values from the `ticket`
(which is an array of field values)
"""
function error_sum(rules, ticket)
    bad_values = filter(v -> error_val(rules, v), ticket)
    return sum(bad_values)
end

function part1(input)
    rule_ranges = collect(Iterators.flatten([r.ranges for r in input.rules]))
    return sum(error_sum(rule_ranges, t) for t in input.nearby_tickets)
end

# ----------------------------------------------------------------------
# Part 2
# ----------------------------------------------------------------------

"""
Return an array of rules that can be applied perfectly to all values.
but ignore the ones in the `excludes` dict.
"""
function candidate_fields(input, values, excludes)
    return [r for r in input.rules
        if all(v ∈ r.ranges[1] || v ∈ r.ranges[2] for v in values) &&
            !haskey(excludes, r)]
end

"Gather values for a specific field from the provided tickets."
field_values(tickets, pos) = [t[pos] for t in tickets]

"Return number of fields in the input."
field_count(input) = length(input.my_ticket)

"Return all rules from the input as a flat array."
all_rules(input) = collect(Iterators.flatten([r.ranges for r in input.rules]))

"""
Return good nearby tickets, which is defined as not having any errors in
any of the fields.
"""
function good_nearby_tickets(input)
    rr = all_rules(input)
    [t for t in input.nearby_tickets if !any(error_val(rr, v) for v in t)]
end

"""
Figure out how each field is mapped to which rule using dynamic programming.

The algorithm uses a Dict to keep track of fields that can be mapped.
When the number of candidate rules for that field is determined to be 1,
then we can map that field right away. Otherwise, we put that field back
into a queue and revisit that later.
"""
function part2_search(input)
    tickets = good_nearby_tickets(input)
    found = Dict{Any,Int}()  # rule -> pos mapping
    len = field_count(input)
    queue = collect(1:len)
    cnt = 0
    while length(found) < len
        pos = popfirst!(queue)
        cf = candidate_fields(input, field_values(tickets, pos), found)
        if length(cf) == 1
            # yay! we found a mapping!
            found[cf[1]] = pos
        elseif length(cf) > 1
            # Push it back to queue and try again later
            push!(queue, pos)
        else
            error("not possible $pos")
        end
        cnt += 1
        cnt > 10000 && break  # just a circuit breaker
    end
    @show length(found) len   # integrity check
    return found
end

"""
Find mapping of fields and rules. Then, look up my ticket and multiply
the values for fields that start with the word `departure`.
"""
function part2()
    input = read_data()
    dct = part2_search(input)
    positions = [v for (k,v) in dct if startswith(k.name, "departure")]
    return prod(input.my_ticket[positions])
end


text:
module AutoMLPipeline

greet() = print("Hello World!")
export fit!, transform!, fit_transform!

include("abstracttypes.jl")
using .AbsTypes
export Machine, Computer, Workflow, Learner, Transformer

include("utils.jl")
using .Utils
export getiris, getprofb

include("baselines.jl")
using .Baselines
export Baseline, Identity

include("basefilters.jl")
using .BaseFilters
export OneHotEncoder, Imputer

include("featureselector.jl")
using .FeatureSelectors
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator

include("skpreprocessor.jl")
using .SKPreprocessors
export SKPreprocessor, skpreprocessors

include("sklearners.jl")
using .SKLearners
export SKLearner, sklearners

include("decisiontree.jl")
using .DecisionTreeLearners
export PrunedTree, RandomForest, Adaboost

include("ensemble.jl")
using .EnsembleMethods
export VoteEnsemble, StackEnsemble, BestLearner

include("crossvalidator.jl")
using .CrossValidators
export crossvalidate

include("skcrossvalidator.jl")
using .SKCrossValidators
export crossvalidate

include("naremover.jl")
using .NARemovers
export NARemover

include("pipelines.jl")
using .Pipelines
export @pipeline, @pipelinex, @pipelinez
export Pipeline, ComboPipeline

end # module


text:
module MatrixCore

include("drivers.jl")

export Driver, HumidityDriver,
    IMUDriver, UVDriver, PressureDriver,
    GPIODriver, EverloopDriver,HumidityDriver,
    printmessage,
    connecterror, connectdata,
    connectkeepalive, configure,
    close, start

end # module


text:
mutable struct Conv1 <: NN
    w::PType4 #weights
    b::PType4 #biases
    ms::KnetMoment #batch norm moments
    bn_p::PType1 #batch norm parameters
    stride::Int64
    padding::Int64
    training::Bool
end

function (c::Conv1)(x::HGType)
    bn = batchnorm(x,c.ms,c.bn_p,training=c.training)
    relu_out = relu.(bn)
    output1 = conv4(c.w,relu_out,stride=c.stride,padding=c.padding)
    output2 = output1 .+ c.b

    myfree(bn); myfree(relu_out); myfree(output1)

    output2
end

function Conv1(in_dim::Int,out_dim::Int,k_s,stride,padding,atype=KnetArray)

    bn_p = bnparams(Float32,in_dim)
    bn_p = Param(convert(atype{Float32,1},bn_p))
    ms = bnmoments()

    w = xavier_normal(Float32,k_s,k_s,in_dim,out_dim)
    b = xavier_normal(Float32,1,1,out_dim,1)
    w = Param(convert(atype,w))
    b = Param(convert(atype,b))
    Conv1(w,b,ms,bn_p,stride,padding,true)
end

function set_testing(f::Conv1,training=false)
    f.training=training
    nothing
end

mutable struct Conv0 <: NN
    w::PType4 #weights
    b::PType4 #biases
    stride::Int64
    padding::Int64
end

function (c::Conv0)(x::HGType)
    out1=conv4(c.w,x,stride=c.stride,padding=c.padding)
    out2 = out1 .+ c.b

    myfree(out1)

    out2
end

function Conv0(in_dim::Int,out_dim::Int,k_s,stride,padding,atype=KnetArray)
    w = xavier_normal(Float32,k_s,k_s,in_dim,out_dim)
    b = xavier_normal(Float32,1,1,out_dim,1)
    w = Param(convert(atype,w))
    b = Param(convert(atype,b))
    Conv0(w,b,stride,padding)
end

struct Residual <: NN
    c1::Conv1
    c2::Conv1
    c3::Conv1
end

function (r::Residual)(x::HGType)
    c1=r.c1(x)
    c2=r.c2(c1)
    output = r.c3(c2) .+ x

    myfree(c1); myfree(c2)

    output
end

function Residual(in_dim::Int,out_dim::Int,atype=KnetArray)
    c1=Conv1(in_dim,div(out_dim,2),1,1,0,atype)
    c2=Conv1(div(out_dim,2),div(out_dim,2),3,1,1,atype)
    c3=Conv1(div(out_dim,2),out_dim,1,1,0,atype)
    Residual(c1,c2,c3)
end

function set_testing(f::Residual,training=false)
    set_testing(f.c1,training)
    set_testing(f.c2,training)
    set_testing(f.c3,training)
    nothing
end

mutable struct Residual_skip <: NN
    w::PType4
    b::PType4
    c1::Conv1
    c2::Conv1
    c3::Conv1
end

function (r::Residual_skip)(x::HGType)
    residual1 = conv4(r.w,x,stride=1)
    residual2 = residual1 .+ r.b
    c1=r.c1(x)
    c2=r.c2(c1)
    output = r.c3(c2) .+ residual2

    myfree(c1); myfree(c2); myfree(residual1); myfree(residual2);

    output
end

function Residual_skip(in_dim::Int,out_dim::Int,atype=KnetArray)
    c1=Conv1(in_dim,div(out_dim,2),1,1,0,atype)
    c2=Conv1(div(out_dim,2),div(out_dim,2),3,1,1,atype)
    c3=Conv1(div(out_dim,2),out_dim,1,1,0,atype)

    w = xavier_normal(Float32,1,1,in_dim,out_dim)
    b = xavier_normal(Float32,1,1,out_dim,1)
    w=Param(convert(atype,w))
    b=Param(convert(atype,b))
    Residual_skip(w,b,c1,c2,c3)
end

function set_testing(f::Residual_skip,training=false)
    set_testing(f.c1,training)
    set_testing(f.c2,training)
    set_testing(f.c3,training)
    nothing
end

struct Pool
end
(p::Pool)(x::HGType) = pool(x)

struct Unpool
end
(u::Unpool)(x::HGType) = unpool(x)


text:
module ContinuousOptimization

using LearningStrategies
const LS = LearningStrategies

#=

First, we need LearningStrategy types to do individual components of Optim.optimize:

function optimize{T, M<:Optimizer}(d, initial_x::Array{T}, method::M, options::OptimizationOptions)

-- LS.TimeLimit
    t0 = time() # Initial time stamp used to control early stopping by options.time_limit

-- ignore?
    if length(initial_x) == 1 && typeof(method) <: NelderMead
        error("Use optimize(f, scalar, scalar) for 1D problems")
    end

-- TODO
    state = initial_state(method, options, d, initial_x)

-- LS.ShowStatus, LS.Tracer, LS.IterFunction
    tr = OptimizationTrace{typeof(method)}()
    tracing = options.store_trace || options.show_trace || options.extended_trace || options.callback != nothing
    stopped, stopped_by_callback, stopped_by_time_limit = false, false, false

- LS.GradConverged
    x_converged, f_converged = false, false
    g_converged = if typeof(method) <: NelderMead
        nmobjective(state.f_simplex, state.m, state.n) < options.g_tol
    elseif  typeof(method) <: ParticleSwarm || typeof(method) <: SimulatedAnnealing
        g_converged = false
    else
        vecnorm(state.g, Inf) < options.g_tol
    end

    converged = g_converged
    iteration = 0

-- part of LS.ShowStatus's pre_hook
    options.show_trace && print_header(method)
    trace!(tr, state, iteration, method, options)

-- ignore... part of MetaLearner
    while !converged && !stopped && iteration < options.iterations

-- LS.MaxIter
        iteration += 1

-- LS.GradientLearner
        update_state!(d, state, method) && break # it returns true if it's forced by something in update! to stop (eg dx_dg == 0.0 in BFGS)
        update_g!(d, state, method)

-- LS.FunctionConverged, LS.GradConverged, LS.ParamsConverged
        x_converged, f_converged,
        g_converged, converged = assess_convergence(state, options)
        # We don't use the Hessian for anything if we have declared convergence,
        # so we might as well not make the (expensive) update if converged == true
        !converged && update_h!(d, state, method)

-- LS.ShowStatus, LS.Tracer
        # If tracing, update trace with trace!. If a callback is provided, it
        # should have boolean return value that controls the variable stopped_by_callback.
        # This allows for early stopping controlled by the callback.
        if tracing
            stopped_by_callback = trace!(tr, state, iteration, method, options)
        end

-- LS.TimeLimit
        # Check time_limit; if none is provided it is NaN and the comparison
        # will always return false.
        stopped_by_time_limit = time()-t0 > options.time_limit ? true : false

-- ignore... part of MetaLearner
        # Combine the two, so see if the stopped flag should be changed to true
        # and stop the while loop
        stopped = stopped_by_callback || stopped_by_time_limit ? true : false
    end # while

    after_while!(d, state, method, options)

-- TODO: how do we nicely summarize the sub-learners??
    return MultivariateOptimizationResults(state.method_string,
                                            initial_x,
                                            state.x,
                                            Float64(state.f_x),
                                            iteration,
                                            iteration == options.iterations,
                                            x_converged,
                                            options.x_tol,
                                            f_converged,
                                            options.f_tol,
                                            g_converged,
                                            options.g_tol,
                                            tr,
                                            state.f_calls,
                                            state.g_calls,
                                            state.h_calls)
end

Here's a rough outline of how we'd do this with MetaLearner:
-----------------------------------------------------------

# keep vector of: ‖θ_true - θ‖
tracer = Tracer(Float64, (model,i) -> norm(θ - params(model)))

# build the MetaLearner
learner = make_learner(
    GradientLearner(...),
    TimeLimit(60), # stop after 60 seconds
    MaxIter(1000), # stop after 1000 iterations
    ShowStatus(100), # show a status update before iterating and every 100 iterations
    Converged(params, tol=1e-6), # similar to x_converged for the function-case
    Converged(output, tol=1e-6), # similar to f_converged for the function-case
    Converged(grad, tol=1e-6, every=10), # similar to g_converged for the function-case
                                         # note: we can also only check every ith iteration
)

# learn is like optimize.
learn!(model, learner)

# note: for the function minimization case, it "iterates" over obs == nothing, since the
#   `x` in f(x) is treated as learnable parameters in Transformations.Differentiable, NOT input

------------------------------------------------------------
=#

end # module


text:
@doc raw"""
AshkinTeller model with energy
``E = -\sum_{ij} J^\sigma_{ij} \sigma_i \sigma_j + J^\tau_{ij} \tau_i \tau_j + K_{ij} \sigma_i \sigma_j \tau_i \tau_j``,
where ``\sigma_i`` and ``\tau_i`` takes value of 1 (up spin) or -1 (down spin).
"""
mutable struct AshkinTeller <: Model
    lat :: Lattice
    spins :: Matrix{Int}
    rng :: Random.MersenneTwister

    function AshkinTeller(lat::Lattice, rng::Random.AbstractRNG)
        model = new()
        model.lat = lat
        model.rng = rng
        model.spins = rand(model.rng, [1,-1], 2, numsites(lat))
        return model
    end
end

AshkinTeller(lat::Lattice) = AshkinTeller(lat, Random.seed!(Random.MersenneTwister(0)))
AshkinTeller(lat::Lattice, seed) = AshkinTeller(lat, Random.seed!(Random.MersenneTwister(0), seed...))

@doc doc"""
    AshkinTeller(param)

Generates `AshkinTeller` using `param["Lattice"]` and `param["Seed"]` (if defined).
Each spin will be initialized randomly and independently.
"""
function AshkinTeller(param::Parameter)
    lat = generatelattice(param)
    if "Seed" in keys(param)
        return AshkinTeller(lat, param["Seed"])
    else
        return AshkinTeller(lat)
    end
end

include("update.jl")
include("estimator.jl")
include("postproc.jl")