text:
#= ======================
   LINEAR REGRESSOR (OLS)
   ====================== =#

"""
$SIGNATURES

Standard linear regression model.

## Parameters

* `fit_intercept` (Bool): whether to fit the intercept or not.
* `solver`: type of solver to use (if `nothing` the default is used). The
            solver is Cholesky by default but can be Conjugate-Gradient as
            well. See `?Analytical` for more information.

"""
@with_kw_noshow mutable struct LinearRegressor <: MMI.Deterministic
    fit_intercept::Bool    = true
    solver::Option{Solver} = nothing
end

glr(m::LinearRegressor) = LinearRegression(fit_intercept=m.fit_intercept)

descr(::Type{LinearRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2``."

#= ===============
   RIDGE REGRESSOR
   =============== =#

"""
$SIGNATURES

Ridge regression model with objective function

``|Xθ - y|₂²/2 + λ|θ|₂²/2``

## Parameters

* `lambda` (Real): strength of the L2 regularisation.
* `fit_intercept` (Bool): whether to fit the intercept or not.
* `penalize_intercept` (Bool): whether to penalize the intercept.
* `solver`: type of solver to use (if `nothing` the default is used). The
            solver is Cholesky by default but can be Conjugate-Gradient as
            well. See `?Analytical` for more information.
"""
@with_kw_noshow mutable struct RidgeRegressor <: MMI.Deterministic
    lambda::Real             = 1.0
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::RidgeRegressor) =
    RidgeRegression(m.lambda,
                    fit_intercept=m.fit_intercept,
                    penalize_intercept=m.penalize_intercept)

descr(::Type{RidgeRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2 + λ|θ|₂²/2``."

#= ===============
   LASSO REGRESSOR
   =============== =#

"""
$SIGNATURES

Lasso regression model with objective function

``|Xθ - y|₂²/2 + λ|θ|₁``

## Parameters

* `lambda` (Real): strength of the L1 regularisation.
* `fit_intercept` (Bool): whether to fit the intercept or not.
* `penalize_intercept` (Bool): whether to penalize the intercept.
* `solver`: type of solver to use (if `nothing` the default is used). Either
            `FISTA` or `ISTA` can be used (proximal methods, with/without
            acceleration).
"""
@with_kw_noshow mutable struct LassoRegressor <: MMI.Deterministic
    lambda::Real             = 1.0
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::LassoRegressor) =
    LassoRegression(m.lambda,
                    fit_intercept=m.fit_intercept,
                    penalize_intercept=m.penalize_intercept)

descr(::Type{LassoRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2 + λ|θ|₁``."

#= =====================
   ELASTIC NET REGRESSOR
   ===================== =#

"""
$SIGNATURES

Elastic net regression model with objective function

``|Xθ - y|₂²/2 + λ|θ|₂²/2 + γ|θ|₁``

## Parameters

* `lambda` (Real): strength of the L2 regularisation.
* `gamma` (Real): strength of the L1 regularisation.
* `fit_intercept` (Bool): whether to fit the intercept or not.
* `penalize_intercept` (Bool): whether to penalize the intercept.
* `solver`: type of solver to use (if `nothing` the default is used). Either
            `FISTA` or `ISTA` can be used (proximal methods, with/without
            acceleration).
"""
@with_kw_noshow mutable struct ElasticNetRegressor <: MMI.Deterministic
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::ElasticNetRegressor) =
    ElasticNetRegression(m.lambda, m.gamma,
                         fit_intercept=m.fit_intercept,
                         penalize_intercept=m.penalize_intercept)

descr(::Type{ElasticNetRegressor}) = "Regression with objective function ``|Xθ - y|₂²/2 + λ|θ|₂²/2 + γ|θ|₁``."

#= ==========================
   ROBUST REGRESSOR (General)
   ========================== =#

"""
$SIGNATURES

Robust regression model with objective function

``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁``

where `ρ` is a robust loss function (e.g. the Huber function).

## Parameters

* `rho` (RobustRho): the type of robust loss to use (see `HuberRho`,
                     `TalwarRho`, ...)
* `penalty` (Symbol or String): the penalty to use, either `:l2`, `:l1`, `:en`
                                (elastic net) or `:none`. (Default: `:l2`)
* `lambda` (Real): strength of the regulariser if `penalty` is `:l2` or `:l1`.
                   Strength of the L2 regulariser if `penalty` is `:en`.
* `gamma` (Real): strength of the L1 regulariser if `penalty` is `:en`.
* `fit_intercept` (Bool): whether to fit an intercept (Default: `true`)
* `penalize_intercept` (Bool): whether to penalize intercept (Default: `false`)
* `solver` (Solver): type of solver to use, default if `nothing`.
"""
@with_kw_noshow mutable struct RobustRegressor <: MMI.Deterministic
    rho::RobustRho           = HuberRho(0.1)
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::RobustRegressor) =
    RobustRegression(m.rho, m.lambda, m.gamma;
                     penalty=Symbol(m.penalty),
                     fit_intercept=m.fit_intercept,
                     penalize_intercept=m.penalize_intercept)

descr(::Type{RobustRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` for a given robust `ρ`."

#= ===============
   HUBER REGRESSOR
   =============== =#

"""
$SIGNATURES

Huber Regression, see `RobustRegressor`, it's the same but with the robust loss
set to `HuberRho`.  The parameters are the same apart from `delta` which
parametrises the `HuberRho` function (radius of the ball within which the loss
is a quadratic loss).
"""
@with_kw_noshow mutable struct HuberRegressor <: MMI.Deterministic
    delta::Real              = 0.5
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::HuberRegressor) =
    HuberRegression(m.delta, m.lambda, m.gamma;
                    penalty=Symbol(m.penalty),
                    fit_intercept=m.fit_intercept,
                    penalize_intercept=m.penalize_intercept)

descr(::Type{HuberRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where `ρ` is the Huber Loss."

#= ==================
   QUANTILE REGRESSOR
   ================== =#

"""
$SIGNATURES

Quantile Regression, see `RobustRegressor`, it's the same but with the robust
loss set to `QuantileRho`.  The parameters are the same apart from `delta`
which parametrises the `QuantileRho` function (indicating the  quantile to use
with default `0.5` for the median regression).
"""
@with_kw_noshow mutable struct QuantileRegressor <: MMI.Deterministic
    delta::Real              = 0.5
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::QuantileRegressor) =
    QuantileRegression(m.delta, m.lambda, m.gamma;
                       penalty=Symbol(m.penalty),
                       fit_intercept=m.fit_intercept,
                       penalize_intercept=m.penalize_intercept)

descr(::Type{QuantileRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where `ρ` is the Quantile Loss."

#= ==================================
   LEAST ABSOLUTE DEVIATION REGRESSOR
   ================================== =#

"""
$SIGNATURES

Least Absolute Deviation regression with with objective function

``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁``

where `ρ` is the absolute loss.

See also `RobustRegressor`.
"""
@with_kw_noshow mutable struct LADRegressor <: MMI.Deterministic
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::LADRegressor) =
    LADRegression(m.lambda, m.gamma;
                  penalty=Symbol(m.penalty),
                  fit_intercept=m.fit_intercept,
                  penalize_intercept=m.penalize_intercept)

descr(::Type{LADRegressor}) = "Robust regression with objective ``∑ρ(Xθ - y) + λ|θ|₂² + γ|θ|₁`` where `ρ` is the Absolute Loss."


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

using Test
using Multivectors
using LinearAlgebra

include("Blades_runtests.jl")
include("KVectors_runtests.jl")

module PG2
  using Multivectors, LinearAlgebra
  @generate_basis("++0")

  const eₒ_index = 3
  const eₒ = e₃
  const e₂ₒ = e₂₃
  const e₁ₒ = e₁₃

  meet(a, b) = a∧b
  point(x::Real, y::Real) = x*e₂ₒ + y*e₁ₒ + one(x)*e₁₂
  direction(x::Real, y::Real) = x*e₂ₒ + y*e₁ₒ
  point( v::KVector{T} ) where T = point(coords(v)[1:2]...)
  line(a,b) = a∨b

  """ project point to euclidean (1eₒ) plane """
  projr(a) = a/(dual(1.0e₃)⋅a)

  function circumcentre( a, b, c )
    ab = line(a, b)
    bc = line(b, c)
    ab_mid = 0.5*(a+b)
    bc_mid = 0.5*(b+c)
    abperp = ab⋅ab_mid
    bcperp = bc⋅bc_mid

    projr(abperp∧bcperp)
  end

end

using .PG2
@testset "PG2" begin
  @test PG2.circumcentre( PG2.point(0.0, 0.0), PG2.point(1.0,0.0), PG2.point(0.0,1.0) ) == PG2.point(0.5, 0.5)
end

module HomogeneousG2
  using Multivectors, LinearAlgebra
  @generate_basis("+++")

  const eₒ_index = 3
  const eₒ = e₃

  ishomog( v::KVector{T} ) where T = filter(x->x==eₒ_index, subspace.(v)) |> isempty
  point( a, b ) = a*e₁+b*e₂ + one(a)*e₃
  point( v::KVector{T} ) where T = ishomog(v) ? v : v + one(T)eₒ
  line( p, q ) = p∧q
  plane( p, q, r ) = p∧q∧r
  unitize( p::KVector{T} ) where T = (α = one(T)eₒ⋅p; iszero(α) ? p : p*inv(alpha))

  function meet(a, b)
    J = inv(1.0e₁₂₃)
    ((a*J)∧(b*J))*J
  end

  """ project to euclidean (1eₒ) plane """
  projr(a) = a/norm(1.0eₒ⋅a)

  function circumcentre( a, b, c )
    mid_ab = 0.5*(a+b)
    mid_bc = 0.5*(b+c)
    ab = (a∧b)
    bc = (b∧c)
    rab = 1.0eₒ⋅ab
    rbc = 1.0eₒ⋅bc
    tri = rab∧rbc
    abperp = rab*inv(tri)
    bcperp = rbc*inv(tri)
    centre_ab = mid_ab∧abperp
    centre_bc = mid_bc∧bcperp
    meet(centre_ab, centre_bc)
  end
end

using .HomogeneousG2
@testset "HomogeneousG2" begin
  pointh = HomogeneousG2.point
  @test HomogeneousG2.circumcentre( pointh(0.0, 0.0), pointh(1.0,0.0), pointh(0.0,1.0) ) == pointh(0.5, 0.5)
end

module CGA
  using Multivectors, LinearAlgebra

  @generate_basis("++++-",false,true)
  const eo = 0.5*(1.0e₄ + 1.0e₅)
  const e∞ = 1.0e₅-1.0e₄

  point( p::K ) where {T, K<:Union{KVector{T,1}, Blade{T,1}}} = one(T)*eo + p + (one(T)/2)*p*p*e∞
  locate( a, M ) = (a⋅M)/M
  circumcentre( a, b, c) = locate(a, (a∧b + b∧c + c∧a)/(a∧b∧c∧e∞))
end

using .CGA

@testset "CGA3D" begin
  e₁, e₂, e₃ = alle( CGA, 5 )[1:3]
  using .CGA: eo, e∞, point, circumcentre
  @test e∞⋅e∞ == 0
  @test eo⋅e∞ == -1.0

  # test circumcentre of triangle

  a = point(0.0e₁+0.0e₂)
  b = point(1.0e₁)
  c = point(1.0e₂)

  C = (a∧b + b∧c + c∧a)/(a∧b∧c∧e∞)
  ccentre = (a⋅C)/C
  @test coords(grade(ccentre,1))[1:2] == [0.5,0.5]
  a = point(0.0e₁+0.0e₂)
  b = point(1.0e₁)
  c = point((sqrt(3.0)/2.0)e₂)
  ccentre = grade(circumcentre(a,b,c), 1)

  @test length(a.B) >= length(prune(a).B)
  @test iszero(prune(a-a))
  @test norm(ccentre - grade(a,1)) ≈ norm(ccentre - grade(b,1)) ≈ norm(ccentre - grade(c,1))

  i = 1.0e₁+2.0e₂; j = 3.0(e₁∧e₃) + 4.0(e₂∧e₃); k = KVector(5.0CGA.e₂₃₄)
  @test i+j == Multivector(i) + Multivector(j) == j+i == Multivector(i) + j == i + Multivector(j)
  @test i+j+k == i+(j+k) == (i+k)+j

  v = i*1.1+j+k+1.2
  @test Multivectors.Δⱼ(v, 1) == involute(v)
  @test Multivectors.Δⱼ(v, 2) == ~v
  @test prune(shirokov_inv(i+j+k)*(i+j+k)) ≈ 1.0
end

module PG3
  using Multivectors
  @generate_basis("+++0")
end
using .PG3

@testset "PGA3D" begin
  e₁, e₂, e₃, e₄ = alle( PG3, 4)[1:4]
  e₁₂ = PG3.e₁₂; e₃₄ = PG3.e₃₄; e₁₂₃₄ = PG3.e₁₂₃₄;
  e₁₂₃ = PG3.e₁₂₃; e₂₃₄ = PG3.e₂₃₄; e₁₃₄ = PG3.e₁₃₄; e₁₂₄ = PG3.e₁₂₄
  𝐼 = e₁₂₃₄
  @test KVector(dual(1𝐼)) == dual(KVector(1𝐼))
  @test Multivector(KVector(dual(1𝐼))) == Multivector(dual(KVector(1𝐼))) == dual(Multivector(1𝐼))

  point(x, y, z) = -Float64(x)*e₂₃₄ + Float64(y)*e₁₃₄ - Float64(z)*e₁₂₄ + 1.0*e₁₂₃
  dir(x, y, z) = -Float64(x)*e₂₃₄ + Float64(y)*e₁₃₄ - Float64(z)*e₁₂₄
  line(p, q) = p∨q
  Line{T} = Union{KVector{T, 2}, Blade{T, 2}}

  # dual quaternion for translation
  # constructed as ratio between point (0,0,0) and (x,y,z)
  translator(x,y,z) = -0.5*dir(x,y,z)*1.0e₁₂₃+1.0
  # rotation constructed with ratio between line a+b and line a.  so that (a+b)/a * a *~(a+b)/a = b
  # a rotation from a to b. lines should intersect at origin or we'll get a translation
  rotator(a::Line, b::Line) = normalize(a+b)/a

  q = rotator(line(point(0,0,0), point(1,0,0)),
              line(point(0,0,0), point(1/sqrt(2), 1/sqrt(2), 0)))

  # rotate (1,0,0.1) by 45 degrees around z axis twice
  rxz = normalize(q*q*point(1,0,0.1)*reverse(q)*reverse(q))
  rxz_xyz = coords(dual(rxz)[1])
  rxz_xyz = rxz_xyz .* (1.0/rxz_xyz[4])
  @test isapprox(rxz_xyz[1], 0.0; atol = 1e-10)
  @test isapprox(rxz_xyz[2], 1.0; atol = 1e-10)
  @test isapprox(rxz_xyz[3], 0.1; atol = 1e-10)

  tx = translator(1,0,0)
  @test tx*point(1,0,10)*reverse(tx) == point(2,0,10)
  @test (tx*tx)*point(1,0,10)*reverse(tx*tx) == point(3,0,10)

  # dual quaternion is composition of translator and rotator
  x = point(1,0,0)
  dq = tx*q  # rotate then translate
  qd = q*tx  # translate then rotate
  dqx = dq*x*reverse(dq)
  # check x coord after transformation
  @test (dqx∧1.0e₁)[4][1] |> scalar ≈ 1/sqrt(2) + 1.0
  qdx = qd*x*reverse(qd)
  @test (qdx∧1.0e₁)[4][1] |> scalar ≈ 2/sqrt(2)

  # z coord should not have changed
  qdxz = qd*point(1.0,0,0.1)*reverse(qd)
  dqxz = dq*point(1.0,0,0.1)*reverse(dq)
  @test (qdxz∧1.0e₃)[4][1] |> scalar ≈ 0.1
  @test (dqxz∧1.0e₃)[4][1] |> scalar ≈ 0.1

  a = e₁(1.0); b = e₂(2.0); c = e₃(3.0); d = e₄(4.0)
  @test typeof(a+b*c) == typeof(Multivector{Float64,2}())
  B1 = KVector(a)
  B2 = a*b + b*c
  B3 = a*b*c + b*c*d
  B4 = KVector(a*b*c*d)
  M13 = B1+B3
  @test grade(M13) == 3
  @test grades(M13) == [1,3]
  @test grade(M13, 2) == 0.0
  @test grade(M13, 1) == B1
  @test grade(M13, 3) == B3
  @test grade(M13, 0) == 0.0
  M013 = M13+42.0
  @test grade(M013, 0) == 42.0
  B = a+b
  @test B*B == a*a+a*b+b*a+b*b
  @test grade((2e₂+3e₃)*(2e₂+3e₃), 0) == 2*2+3*3

  # degenerate wedge
  D = 1.0 + 2.0e₄
  E = 2.0 + 3.0(e₁∧e₄)
  @test grade(D∧E, 3) == grade(D*E, 3)
  F = 1.0 + 2.0e₃
  G = 2.0 + 3.0(e₁∧e₃)
  @test grade(F∧G, 3) == grade(F*G, 3)
  @test grade(D∧G, 3) == grade(D*G, 3)

  A = 2e₃ +3e₂
  B = 5e₄+6e₃
  # if we had an e₁ (degen metric) component this would not be true.
  @test A⋅!B == !(A∧B)

  A₀ = 2e₁ + 3e₂
  B₀ = 5e₁ + 6e₃

  # interesting though, this is ok.
  𝐼 = 1e₁₂₃₄
  @test A₀⋅(B₀⋅𝐼) == (A₀∧B₀)⋅𝐼 == (A₀∧B₀)*𝐼 == A₀⋅(B₀*𝐼)# == A₀⋅(B₀*dual(1))

  C = 1e₁₂+4e₃₄
  @test (A∧B)⋅C == A⋅(B⋅C)

  # inner product tests
  # conjugate symmetry? a property an innerproduct is supposed to have.  probably not
  #  @test A⋅C == conj(C⋅A)
  # linearity in first argument
  @test (A+B)⋅C == (A⋅C)+(B⋅C)

  @test 2*A×B == (A*B-B*A)
  A = 2-A
  B = 2*B + 3(e₂∧e₃)
  @test 2A×B == (A*B-B*A)

  v = 1.0+2.0e₁₂ + 1.1(e₁∧e₃)
  @test isapprox(norm(v*inv(v)), 1.0; atol=1e-4)
end

module G4
  using Multivectors
  @generate_basis("++++")
end
using .G4

@testset "conversion" begin
  e₁, e₂, e₃, e₄  = alle( G4, 4)[1:4]
  @test  1.0e₁ + 2e₁ == 1e₁ + 2.0e₁ == 1.0e₁ + 2.0e₁
  @test  1.0e₁ + 2e₂ == 1e₁ + 2.0e₂ == 1.0e₁ + 2.0e₂
end

@testset "contraction" begin

  e₁, e₂, e₃, e₄  = alle( G4, 4)[1:4]
  e₁₃ = G4.e₁₃

  rand1() = rand(1)[1]
  a = rand1()∧e₁ + rand()∧e₃
  # test lcontraction result is orthogonal to original
  B = rand()*(e₁∧e₂)
  @test a⋅(a⋅B) == 0.0
  B = B + rand()*(e₃∧e₄)
  @test a⋅(a⋅B) == 0.0

  α = Multivector(3.0)
  β = 3.0

  A = 1.0 + 3.0G4.e₁₃
  @test lcontraction(α, β) == lcontraction(β, α) == lcontraction(α, α) == 9.0
  @test rcontraction(α, β) == rcontraction(β, α) == rcontraction(α, α) == 9.0
  @test lcontraction(α, A) == α*A
  @test lcontraction(A, α) == α*A[0]
  @test rcontraction(α, A) == α*A[0]
  @test rcontraction(A, α) == A*α
  @test ∙(α, A) == ∙(A, α) == 0.0
  @test ∙(a, a) == grade(a*a, 0)
  @test A*a == A[0]*a + A[2]*a
  @test lcontraction(a,A) == a⋅A == grade(a, 1)*A[2]

  C = A + a
  D = a + 2.0*pseudoscalar(e₁)
  @test scalarprod(A∧C, D) == scalarprod(A, lcontraction(C,D))
  @test scalarprod(D, C∧A) == scalarprod(rcontraction(D,C), A)
  @test_broken rcontraction(C,D) == C[1]*D[1] + C[2]*D[1]  #!me not sure about the new right contraction, comes from mix of Lengyal and Browne
  @test_broken rcontraction(D,C) == D*C[0] + D[1]*C[1] + D[4]*C[1] + D[4]*C[2]
  @test lcontraction(C,D) == C[0]*D + C[1]*D[1] + C[1]*D[4] + C[2]*D[4]
  @test_broken reverse(lcontraction(A,C)) == rcontraction(reverse(C), reverse(A))
  @test_broken reverse(lcontraction(C,D)) == rcontraction(reverse(D), reverse(C))
  @test_broken reverse(lcontraction(A,D)) == rcontraction(reverse(D), reverse(A))
  @test_broken reverse(lcontraction(D,A)) == rcontraction(reverse(A), reverse(D))
  @test grades(A) == [0,2]
  show(A)
  @test first(A) == 1.0
  @test [i for i in A] == map(i->i, A) |>collect == (i for i in A) |> collect
  @test length(A) == 2
  @test mapreduce(identity, +, A) == A
  @test mapreduce(conj, +, A) == conj(A)
  @test 1.0e₁ + A[2] + A[0] == A + 1.0e₁
  @test iszero(grade(A,10))
  @test grades(1.0) == [0]
  @test grade(1.0, 0) == 1.0
  @test grade(1.0, 1) == 0.0
  @test grades(A[2]) == [2]
  @test prune(dual(A)*pseudoscalar(A[2])) == A
  @test A/2.0 == A*0.5
  @test A[2] == A-A[0]
  @test mapreduce(LinearAlgebra.norm_sqr, +, A) == scalarprod(A,reverse(A)) == 10.0
  @test normalize(A) == A/sqrt(scalarprod(A,reverse(A)))

  for i in A
    show(i)
  end

  a = 1.0 + 1.5(e₁) + 2.0(e₁∧e₂)
  b = 1.0 - 10.0e₁ - 10.0e₄ + 2.5(e₁∧e₂∧e₄) + 3.0(e₁∧e₂∧e₃∧e₄)
  @test a⋅b != grade(a*b, grade(b)-grade(a)) # only works for homogeneous k-vectors
  @test a∧b != grade(a*b, grade(b)+grade(a)) # only works for homogeneous k-vectors

  a = 1.0e₁+2.0e₂+3.0e₃
  B = b
  C = 1.0e₁ + 2.0G4.e₁₂ + 3.0G4.e₁₂₃ + 1.1e₂ + 2.2G4.e₂₄
  @test a⋅(B∧C) == (a⋅B)∧C + involute(B)∧(a⋅C)
  @test_broken rcontraction(B, a) == -a⋅involute(B)
  @test lcontraction(a, B) == 0.5*(a*B - involute(B)*a)
  @test_broken rcontraction(B, a) == 0.5*(B*a - a*involute(B))

  # less common multivector products
  A,B = Multivector(2.0e₁+1.0e₁₃), Multivector(2.0e₁+3.0e₂)
  @test length(A×B) > 0
  @test length(B⋅A) == 2
  @test length(A⋅B) == 1
  @test length(B∧A) == 2
  @test length(A∧B) == 2
  @test 2.0⋅A == 2.0*A
  @test 2.0∧A == 2.0*A
  @test A⋅2.0 == 0.0

  @test (1e₁ ∈ 1G4.e₁₂) == true
  @test (1e₁ ∈ 1G4.e₂₃+1G4.e₁₂) == false
  @test (1+1e₁ ∈ 1G4.e₂₃+1G4.e₁₂) == false

  @test isapprox(newton_inv(10.0+2.0G4.e₂₄+3.0G4.e₂₃+0.1e₂)*(10.0+2.0G4.e₂₄+3.0G4.e₂₃+0.1e₂), Multivector(1.0); atol = 0.001)
end

module G3
  using Multivectors
  @generate_basis("+++",false,true,true)
end
using .G3

@testset "Outermorphism" begin
  e₁, e₂, e₃, e₁₂, e₁₃, e₂₃, e₁₂₃ = alle( G3, 3)
  b = 3.0e₁₃
  k = 1.0e₁ + 2.0e₂
  M = 42.0 + b + k
  scale = [2.0 0 0; 0 0.75 0; 0 0 10.0]
  @test outermorphism(scale, b) == 10.0*2.0*b
  @test prune(outermorphism(scale, M)[1]) == 2.0e₁ + 0.75*2.0e₂
end

@testset "Quaternion" begin

  e₁, e₂, e₃, e₁₂, e₁₃, e₂₃, e₁₂₃ = alle( G3, 3)
#==
  Quaternion multiplication
×	i	j	k
i	−1	k	−j
j	−k	−1	i
k	j	−i	−1
==#

  # quaternion basis
  𝑖 = e₂₃; 𝑗 = e₁₃; 𝑘 = e₁₂
  @test 1𝑖*1𝑗*1𝑘 == -1

  @test (1𝑖*1𝑖, 1𝑖*1𝑗, 1𝑖*1𝑘) == (-1, 1𝑘, -1𝑗)
  @test (1𝑗*1𝑖, 1𝑗*1𝑗, 1𝑗*1𝑘) == (-1𝑘, -1, 1𝑖)
  @test (1𝑘*1𝑖, 1𝑘*1𝑗, 1𝑘*1𝑘) == (1𝑗, -1𝑖, -1)

  𝑖 = e₁₂; 𝑗 = e₂₃; 𝑘 = e₁₃
  @test 1𝑖*1𝑗*1𝑘 == -1
  @test (1𝑖*1𝑖, 1𝑖*1𝑗, 1𝑖*1𝑘) == (-1, 1𝑘, -1𝑗)
  @test (1𝑗*1𝑖, 1𝑗*1𝑗, 1𝑗*1𝑘) == (-1𝑘, -1, 1𝑖)
  @test (1𝑘*1𝑖, 1𝑘*1𝑗, 1𝑘*1𝑘) == (1𝑗, -1𝑖, -1)

  half45 = normalize(1.0e₁ + normalize(1.0e₁+1.0e₂))
  q = half45/1.0e₁

  # Transform a 1-vector with the sandwich product.
  v = reverse(q)*(1.0e₁+1.0e₂+1.0e₃)*q

  v′ = grade(v, 1) |> prune∘sort_basis
  @test v′⋅1.0e₃ == 1.0
  @test v′⋅1.0e₁ ≈ sqrt(2.0)

  # Rotors can be constructed using half-angle of trig functions, like quaternions.

  @test cos(π/8) - sin(π/8)*1.0e₁₂ == q

  @test isapprox(grade(log(exp((π/10)*𝑖)), 2), π/10*𝑖; atol = 0.0001)
  @test isapprox(grade(log(exp((π/10)*𝑖)), 0), 0.0; atol = 0.0001)

end

@testset "Barycentric" begin

  e₁, e₂, e₃  = alle( G3, 3)[1:3]

  a = 0.0e₁+0.0e₂; b = 1.0e₁ + 0.0e₂; c = 0.0e₁ + 1.0e₂;  # a simple right angle triangle

  A = (b-a)∧(c-a)  # twice the area of the triangle. we don't worry about the factor of 2

  # Make a function to calculate barycentric coords as the ratio of the area of a triangle made with a point `p` and an edge over original triangle.  i.e. the barycentric coord for vertex `a` is the ratio Δpbc/Δabc

  barycoords(p) = ((c-b)∧(p-b)/A, (a-c)∧(p-c)/A, (b-a)∧(p-a)/A)  # a tuple of coords

  # Notice how the code very directly represents the geometric relationship.  The body of the function is also coordinate free ( we never index into the points or vertices ).

  @test barycoords(0.0e₁)[1] == 1.0
  @test barycoords(1.0e₁)[2] == 1.0
  @test barycoords(0.5e₁+0.5e₂) == (0.0, 0.5, 0.5)

  baryscalars(p) = map(k->grade(k, 0), barycoords(p))

  @test baryscalars(0.1e₁+0.25e₂+10.0e₃) == (0.65, 0.1, 0.25)
  @test baryscalars(1.0e₁) == (0.0, 1.0, 0.0)

  d = 1.0e₃
  V = A∧d

  # tetrahedron
  barycoords4(p) = ((c-b)∧(p-b)∧(d-c)/V,
                    (a-c)∧(p-c)∧(d-a)/V,
                    (b-a)∧(p-a)∧(d-a)/V,
                    (b-a)∧(p-a)∧(a-c)/V)

  @test barycoords4(0.25e₁+0.25e₂+0.25e₃) == (0.25, 0.25, 0.25, 0.25)
end

@testset "Cayley" begin
  e₁, e₂, e₃  = alle( G3, 3)[1:3]

  γ₁ = matrix_representation(e₁)
  γ₂ = matrix_representation(e₂)
  γ₁₂ = matrix_representation(G3.e₁₂)

  @test γ₂ ≈ -γ₁₂*γ₁

  @test cayley_table(e₁)[2,4] == 1G3.e₁₃
end

ext: jl


text:
module MiscTest

import NaNMath

using Compat
using Base.Test
using ForwardDiff

include(joinpath(dirname(@__FILE__), "utils.jl"))

##########################
# Nested Differentiation #
##########################

# README example #
#----------------#

x = rand(5)

f = x -> sum(sin, x) + prod(tan, x) * sum(sqrt, x)
g = x -> ForwardDiff.gradient(f, x)
j = x -> ForwardDiff.jacobian(g, x)

@test_approx_eq ForwardDiff.hessian(f, x) j(x)

# higher-order derivatives #
#--------------------------#

function tensor(f, x)
    n = length(x)
    out = ForwardDiff.jacobian(y -> ForwardDiff.hessian(f, y), x)
    return reshape(out, n, n, n)
end

test_tensor_output = reshape([240.0  -400.0     0.0;
                             -400.0     0.0     0.0;
                                0.0     0.0     0.0;
                             -400.0     0.0     0.0;
                                0.0   480.0  -400.0;
                                0.0  -400.0     0.0;
                                0.0     0.0     0.0;
                                0.0  -400.0     0.0;
                                0.0     0.0     0.0], 3, 3, 3)

@test_approx_eq tensor(rosenbrock, [0.1, 0.2, 0.3]) test_tensor_output

test_nested_jacobian_output = [-sin(1)  0.0     0.0;
                               -0.0    -0.0    -0.0;
                               -0.0    -0.0    -0.0;
                                0.0     0.0     0.0;
                               -0.0    -sin(2) -0.0;
                               -0.0    -0.0    -0.0;
                                0.0     0.0     0.0;
                               -0.0    -0.0    -0.0;
                               -0.0    -0.0    -sin(3)]

@test_approx_eq ForwardDiff.jacobian(x -> ForwardDiff.jacobian(sin, x), [1, 2, 3]) test_nested_jacobian_output

# Issue #59 example #
#-------------------#

x = rand(2)

f = x -> sin(x)/3 * cos(x)/2
df = x -> ForwardDiff.derivative(f, x)
testdf = x -> (((cos(x)^2)/3) - (sin(x)^2)/3) / 2
f2 = x -> df(x[1]) * f(x[2])
testf2 = x -> testdf(x[1]) * f(x[2])

@test_approx_eq ForwardDiff.gradient(f2, x) ForwardDiff.gradient(testf2, x)

######################################
# Higher-Dimensional Differentiation #
######################################

x = rand(5, 5)

@test_approx_eq ForwardDiff.jacobian(inv, x) -kron(inv(x'), inv(x))

#########################################
# Differentiation with non-Array inputs #
#########################################

x = rand(5,5)

# Sparse
f = x -> sum(sin, x) + prod(tan, x) * sum(sqrt, x)
gfx = ForwardDiff.gradient(f, x)
@test_approx_eq gfx ForwardDiff.gradient(f, sparse(x))

# Views
jinvx = ForwardDiff.jacobian(inv, x)
@test_approx_eq jinvx ForwardDiff.jacobian(inv, Compat.view(x, 1:5, 1:5))

########################
# Conversion/Promotion #
########################

# target function returns a literal (Issue #71) #
#-----------------------------------------------#

@test ForwardDiff.derivative(x->zero(x), rand()) == ForwardDiff.derivative(x->1.0, rand())
@test ForwardDiff.gradient(x->zero(x[1]), [rand()]) == ForwardDiff.gradient(x->1.0, [rand()])
@test ForwardDiff.hessian(x->zero(x[1]), [rand()]) == ForwardDiff.hessian(x->1.0, [rand()])
@test ForwardDiff.jacobian(x->[zero(x[1])], [rand()]) == ForwardDiff.jacobian(x->[1.0], [rand()])

# arithmetic element-wise functions #
#-----------------------------------#

N = 4
a = ones(N)
jac0 = reshape(vcat([[zeros(N*(i-1)); a; zeros(N^2-N*i)] for i = 1:N]...), N^2, N)

for op in (-, +, .-, .+, ./, .*)
    f = x -> [op(x[1], a); op(x[2], a); op(x[3], a); op(x[4], a)]
    jac = ForwardDiff.jacobian(f, a)
    @test_approx_eq jac0 jac
end

# NaNs #
#------#

@test ForwardDiff.partials(NaNMath.pow(ForwardDiff.Dual(-2.0,1.0),ForwardDiff.Dual(2.0,0.0)),1) == -4.0

end


text:
using DiffEqSensitivity, OrdinaryDiffEq, Zygote, LinearAlgebra, FiniteDiff, Test
A  = [1.0*im  2.0; 3.0 4.0]
u0 = [1.0 0.0*im; 0.0 1.0]
tspan = (0.0, 1.0)

function f(u,p,t)
    (A*u)*(p[1]*t + p[2]*t^2 + p[3]*t^3 + p[4]*t^4)
end

p = [1.5 + im, 1.0, 3.0, 1.0]
prob = ODEProblem{false}(f,u0,tspan,p)

utarget = [0.0*im 1.0; 1.0 0.0]

function loss_adjoint(p)
    ufinal = last(solve(prob, Tsit5(), p=p, abstol=1e-12, reltol=1e-12, sensealg = InterpolatingAdjoint()))
    loss = 1 - abs(tr(ufinal*utarget')/2)^2
    return loss
end

grad1 = Zygote.gradient(loss_adjoint,Complex{Float64}[1.5, 1.0, 3.0, 1.0])[1]
grad2 = FiniteDiff.finite_difference_gradient(loss_adjoint,Complex{Float64}[1.5, 1.0, 3.0, 1.0])
@test grad1 ≈ grad2

function rhs(u, p, t)
    p .* u
end

function loss_fun(sol)
    final_u = sol[:, end]
    err = sum(abs.(final_u))
    return err
end

function inner_loop(prob, p, loss_fun; sensealg = InterpolatingAdjoint())
    sol = solve(prob, Tsit5(), p=p,  saveat=0.1; sensealg)
    err = loss_fun(sol)
    return err
end

tspan = (0.0, 1.0)
p = [1.0]
u0=[1.0, 2.0]
prob = ODEProblem(rhs, u0, tspan, p)
grads = Zygote.gradient((p)->inner_loop(prob, p, loss_fun), p)[1]

u0=[1.0 + 2.0*im, 2.0 + 1.0*im]
prob = ODEProblem(rhs, u0, tspan, p)
dp1 = Zygote.gradient((p)->inner_loop(prob, p, loss_fun), p)[1]
dp2 = Zygote.gradient((p)->inner_loop(prob, p, loss_fun; sensealg = QuadratureAdjoint()), p)[1]
dp3 = Zygote.gradient((p)->inner_loop(prob, p, loss_fun; sensealg = BacksolveAdjoint()), p)[1]
@test dp1 ≈ dp2 ≈ dp3
@test eltype(dp1) <: Float64


text:
"""
    Model for testing steady state solver with non-linear equations
"""
module E7

using ModelBaseEcon

model = Model()
model.substitutions = true

@parameters model begin
    delta = 0.1000000000000000
    p_dlc_ss = 0.0040000000000000
    p_dlinv_ss = 0.0040000000000000
    p_growth = 0.0040000000000000
end

@variables model begin
    dlc; dlinv; dly; lc; linv;
    lk; ly;
end

@shocks model begin
    dlc_shk; dlinv_shk;
end

@autoexogenize model begin
    lc = dlc_shk
    linv = dlinv_shk
end

@equations model begin
   dlc[t]=(1-0.2-0.2)*p_dlc_ss+0.2*dlc[t-1]+0.2*dlc[t+1]+dlc_shk[t]
   dlinv[t]=(1-0.5)*p_dlinv_ss+0.1*dlinv[t-2]+0.1*dlinv[t-1]+0.1*dlinv[t+1]+0.1*dlinv[t+2]+0.1*dlinv[t+3]+dlinv_shk[t]
   lc[t]=lc[t-1]+dlc[t]
   linv[t]=linv[t-1]+dlinv[t]
   ly[t]=log(exp(lc[t])+exp(linv[t]))
   dly[t]=ly[t]-ly[t-1]
   lk[t]=log((1-delta)*exp(lk[t-1])+exp(linv[t]))
end

@initialize model

end


text:
using BinDeps

@BinDeps.setup

xgboost = library_dependency("xgboost", aliases = ["libxgboost.so"])

if haskey(ENV, "XGBOOST_BUILD_VERSION") && ENV["XGBOOST_BUILD_VERSION"] == "master"
    libcheckout = `git checkout master`
    onload = "global const build_version = \"master\""
    info("Using the latest master version of the XGBoost library")
else
    libcheckout = `git checkout v0.60`
    onload = "global const build_version = \"0.60\""
    info("Using the latest stable version (0.60) of the XGBoost library")
end

provides(BuildProcess,
    FileRule(joinpath(BinDeps.libdir(xgboost), "libxgboost.so"),
        @build_steps begin
            CreateDirectory(BinDeps.srcdir(xgboost))
            @build_steps begin
                ChangeDirectory(BinDeps.srcdir(xgboost))
                `rm -rf xgboost`
                `git clone https://github.com/dmlc/xgboost.git --recursive`
            end

            @build_steps begin
                ChangeDirectory(joinpath(BinDeps.srcdir(xgboost), "xgboost"))
                libcheckout
                `bash build.sh`
                CreateDirectory(BinDeps.libdir(xgboost))
                `cp lib/libxgboost.so $(BinDeps.libdir(xgboost))`
            end
        end), xgboost, os = :Unix, onload = onload)

@BinDeps.install Dict(:xgboost => :_jl_libxgboost)


text:
@testset "_getdatesets" begin
    st = Date(2019, 1, 1)
    ed = Date(2019, 2, 1)
    date_range = st:Day(1):ed

    @testset "easy partition" begin

        holdout_dates = date_range[1:10]
        val_dates = date_range[11:end]

        result = DateSelectors._getdatesets(date_range, holdout_dates)

        @test result isa NamedTuple{(:validation, :holdout)}
        @test result.validation == first(result)
        @test result.holdout == last(result)

        @test Base.isiterable(typeof(result.validation))
        @test Base.isiterable(typeof(result.holdout))

        @test result.holdout == holdout_dates
        @test result.validation == val_dates
    end

    @testset "single date" begin
        date = rand(date_range)
        result = DateSelectors._getdatesets(date_range, date)

        @test result.holdout == [date]
        @test isempty(intersect(result...))
        @test sort(union(result...)) == date_range
    end

    @testset "random partition" begin
        dates = rand(date_range, 10)
        result = DateSelectors._getdatesets(date_range, dates)

        @test result.holdout != dates
        @test result.holdout == sort(dates)
        @test isempty(intersect(result...))
        @test sort(union(result...)) == date_range
    end

    @testset "parsing" begin
        # Just the type
        @test_throws ArgumentError parse(DateSelector, "NoneSelector")
        # Misspelling
        @test_throws ArgumentError parse(DateSelector, "NonSelector()")
        # Wrong expression
        @test_throws ArgumentError parse(DateSelector, "NoneSelector() isa NoneSelector")
        # Extra stuff in the string
        @test_throws Base.Meta.ParseError parse(DateSelector, "NoneSelector() 1")

        # Wrong args
        @test_throws MethodError parse(DateSelector, "PeriodicSelector()")
        @test_throws MethodError parse(DateSelector, "PeriodicSelector(14)")

        # Spacing
        @test parse(DateSelector, "   NoneSelector(  ) ") isa NoneSelector

        s = parse(DateSelector, "NoneSelector()")
        @test s == NoneSelector()
        # Check specific type works
        s = parse(NoneSelector, "NoneSelector()")
        @test s == NoneSelector()

        s = parse(DateSelector, "PeriodicSelector(Day(2))")
        @test s == PeriodicSelector(Day(2))
        s = parse(DateSelector, "PeriodicSelector(Day(6), Day(4), Day(1))")
        @test s == PeriodicSelector(Day(6), Day(4), Day(1))
        # Check specific type works
        s = parse(PeriodicSelector, "PeriodicSelector(Day(2))")
        @test s == PeriodicSelector(Day(2))

        s = parse(DateSelector, "RandomSelector(1)")
        @test s == RandomSelector(1, 1//2, Day(1), Day(0))
        s = parse(DateSelector, "RandomSelector(123, 1//4, Day(7),  Day(1))")
        @test s == RandomSelector(123, 1//4, Day(7), Day(1))
        # Check specific type works
        s = parse(RandomSelector, "RandomSelector(1)")
        @test s == RandomSelector(1, 1//2, Day(1), Day(0))

        s = PeriodicSelector(Day(2))
        @test s == parse(DateSelector, string(s))
    end

    @testset "show" begin
        if VERSION < v"1.5"
            # Check we have undone the period printing for DateSelectors alone
            @test repr(PeriodicSelector(Day(6), Day(4), Day(1))) != "PeriodicSelector(6 days, 4 days, 1 day)"
            @test repr((Day(6), Day(4), Day(1))) == "(6 days, 4 days, 1 day)"
        end

        @test repr(PeriodicSelector(Day(2))) == "PeriodicSelector(Day(2), Day(1), Day(0))"
        @test repr(RandomSelector(1)) == "RandomSelector(1, 1//2, Day(1), Day(0))"
        @test repr(NoneSelector()) == "NoneSelector()"
    end
end


text:
module TestExamples1DEuler

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 "Compressible Euler" begin
  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [2.2527950196212703e-8, 1.8187357193835156e-8, 7.705669939973104e-8],
      linf = [1.6205433861493646e-7, 1.465427772462391e-7, 5.372255111879554e-7])
  end

  @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_pure_fv.jl"),
      l2   = [0.019355699748523896, 0.022326984561234497, 0.02523665947241734],
      linf = [0.02895961127645519,  0.03293442484199227,  0.04246098278632804])
  end

  @trixi_testset "elixir_euler_density_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [0.0011482554820185795, 0.00011482554830363504, 5.741277417754598e-6],
      linf = [0.004090978306814375, 0.00040909783135059663, 2.045489209479001e-5])
  end

  @trixi_testset "elixir_euler_density_wave.jl with initial_condition_density_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [0.0037246420494099768, 0.0037246420494097534, 0.0018623210247056341],
      linf = [0.01853878721991542, 0.018538787219892106, 0.009269393609914633],
      initial_condition = initial_condition_density_pulse)
  end

  @trixi_testset "elixir_euler_density_wave.jl with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [7.71293052584723e-16, 1.9712947511091717e-14, 7.50672833504266e-15],
      linf = [3.774758283725532e-15, 6.733502644351574e-14, 2.4868995751603507e-14],
      initial_condition = initial_condition_constant)
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [3.8099996914101204e-6, 1.6745575717106341e-6, 7.732189531480852e-6],
      linf = [1.2971473393186272e-5, 9.270328934274374e-6, 3.092514399671842e-5])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.11915540925414216, 0.15489191247295198, 0.44543052524765375],
      linf = [0.2751485868543495, 0.2712764982000735, 0.9951407418216425])
  end

  @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.07905582221868049, 0.10180958900546237, 0.29596551476711125],
      linf = [0.23515297345769826, 0.2958208108392532, 0.8694224308790321],
      maxiters = 10,
      surface_flux = flux_kennedy_gruber,
      volume_flux = flux_kennedy_gruber)
  end

  @trixi_testset "elixir_euler_ec.jl with flux_shima_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.07909267609417114, 0.1018246500951966, 0.2959649187481973],
      linf = [0.23631829743146504, 0.2977756307879202, 0.8642794698697331],
      maxiters = 10,
      surface_flux = flux_shima_etal,
      volume_flux = flux_shima_etal)
  end

  @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.07905306555214126, 0.10181180378499956, 0.2959171937479504],
      linf = [0.24057642004451651, 0.29691454643616433, 0.886425723870524],
      maxiters = 10,
      surface_flux = flux_chandrashekar,
      volume_flux = flux_chandrashekar)
  end

  @trixi_testset "elixir_euler_ec.jl with flux_hll" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.07959780803600519, 0.10342491934977621, 0.2978851659149904],
      linf = [0.19228754121840885, 0.2524152253292552, 0.725604944702432],
      maxiters = 10,
      surface_flux = flux_hll,
      volume_flux = flux_ranocha)
  end

  @trixi_testset "elixir_euler_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.11665968950973675, 0.15105507394693413, 0.43503082674771115],
      linf = [0.1867400345208743, 0.24621854448555328, 0.703826406555577])
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [1.250005061244617, 0.06878411345533507, 0.9264328311018613],
      linf = [2.9766770877037168, 0.16838100902295852, 2.6655773445485798])
  end

  @trixi_testset "elixir_euler_sedov_blast_wave_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_pure_fv.jl"),
      l2   = [1.075075094036344, 0.06766902169711514, 0.9221426570128292],
      linf = [3.3941512671408542, 0.16862631133303882, 2.6572394126490315])
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl with pressure" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [1.297525985166995, 0.07964929522694145, 0.9269991156246368],
      linf = [3.1773015255764427, 0.21331831536493773, 2.6650170188241047],
      shock_indicator_variable = pressure,
      cfl = 0.2)
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl with density" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [1.2798798835860528, 0.07103461242058921, 0.9273792517187003],
      linf = [3.1087017048015824, 0.17734706962928956, 2.666689753470263],
      shock_indicator_variable = density,
      cfl = 0.2)
  end

  @trixi_testset "elixir_euler_positivity.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_positivity.jl"),
      l2   = [1.6493820253458906, 0.19793887460986834, 0.9783506076125921],
      linf = [4.71751203912051, 0.5272411022735763, 2.7426163947635844])
  end

  @trixi_testset "elixir_euler_blast_wave.jl" begin
  @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
      l2   = [0.21651329948737183, 0.28091709900008616, 0.5580778880050432],
      linf = [1.513525457073142, 1.5328754303137992, 2.0467706106669556],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"),
        l2   = [2.13605618e-01, 2.79953055e-01, 5.54424459e-01],
        linf = [1.55151701e+00, 1.55696782e+00, 2.05525953e+00],
        maxiters = 30)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"),
        l2   = [2.18148857e-01, 2.83182959e-01, 5.59096194e-01],
        linf = [1.62706876e+00, 1.61680275e+00, 2.05876517e+00],
        maxiters = 30)
  end
end

end # module


text:
using BinaryBuilder, Pkg

name = "Enzyme"
repo = "https://github.com/wsmoses/Enzyme.git"

version = v"0.0.1"

# Collection of sources required to build attr
sources = [GitSource(repo, "307f98aa9619b697120e78e2a4255e36d7b64993")]

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# Bash recipe for building across all platforms
script = raw"""
cd Enzyme
# install_license LICENSE.TXT

CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=RelWithDebInfo)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Explicitly use our cmake toolchain file and tell CMake we're cross-compiling
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING:BOOL=ON)

# Tell CMake where LLVM is
CMAKE_FLAGS+=(-DLLVM_DIR="${prefix}/lib/cmake/llvm")

# Force linking against shared lib
CMAKE_FLAGS+=(-DLLVM_LINK_LLVM_DYLIB=ON)

# Build the library
CMAKE_FLAGS+=(-DBUILD_SHARED_LIBS=ON)
cmake -B build -S enzyme -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["LLVMEnzyme-9", "LLVMEnzyme"], :libEnzyme),
]

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
function mh_sampler(x::Float64, g::Function, sd::Float64, gx::Float64)
  x1 = x + randn()*sd
  gx1 = g(x1)
  if gx1 - gx > log(rand())
    x = x1
    gx = gx1
  end
  return x,gx
end

function mh_sampler(x::Float64, g::Function, sd::Float64)
  gx = g(x)
  mh_sampler(x,g,sd,gx)
end

mh_sampler(x::Float64, g::Function) = mh_sampler(x,g,.5)

# Interface using Density type
mh_sampler(x::Float64, g::Density, sd::Float64, gx::Float64) = mh_sampler(x, Density.f, sd, gx)


text:
module BladeRF

export  devices,
        open,
        close,
        BladeRFDeviceInfo,
        BladeRFDevice,
        TX,
        RX,
        enable,
        disable,
        config,
        receive,
        samplerate

const libbladerf = "libbladerf"
typealias BladeRFStatus Cint

immutable BladeRFModule
    mod::Cint
end

const RX = BladeRFModule(0)
const TX = BladeRFModule(1)

# Thin wrapper around ccall to handle libbladerf's return/error codes
macro bladerfcall( fname, argtypes, args... )
    quote
        ret = ccall( ($fname, libbladerf), BladeRFStatus, $argtypes, $(args...) )
        ret >= 0 || error( "[$ret] "errorstring(ret) )
        ret
    end
end

immutable SerialNumber
    char1::Uint8
    char2::Uint8
    char3::Uint8
    char4::Uint8
    char5::Uint8
    char6::Uint8
    char7::Uint8
    char8::Uint8
    char9::Uint8
    char10::Uint8
    char11::Uint8
    char12::Uint8
    char13::Uint8
    char14::Uint8
    char15::Uint8
    char16::Uint8
    char17::Uint8
    char18::Uint8
    char19::Uint8
    char20::Uint8
    char21::Uint8
    char22::Uint8
    char23::Uint8
    char24::Uint8
    char25::Uint8
    char26::Uint8
    char27::Uint8
    char28::Uint8
    char29::Uint8
    char30::Uint8
    char31::Uint8
    char32::Uint8
    char33::Uint8
    SerialNumber() = new()
end

function Base.string( serial::SerialNumber )
    ascii( [serial.char1,serial.char2,serial.char3,serial.char4,serial.char5,serial.char6,serial.char7,serial.char8,serial.char9,serial.char10,serial.char11,serial.char12,serial.char13,serial.char14,serial.char15,serial.char16,serial.char17,serial.char18,serial.char19,serial.char20,serial.char21,serial.char22,serial.char23,serial.char24,serial.char25,serial.char26,serial.char27,serial.char28,serial.char29,serial.char30,serial.char31,serial.char32,serial.char33] )
end

immutable BladeRFDeviceInfo
    backend::Int32
    serial::SerialNumber
    usb_bus::Uint8
    usb_addr::Uint8
    instance::Uint32
end

type BladeRFDevice
    handle::Ptr{Void}
    function BladeRFDevice( handle )
        device = new( handle )
        finalizer( device, close )
        return device
    end
end

function Base.show( io::IO, devinfo::BladeRFDeviceInfo )
    println( "Backend:      $(devinfo.backend)" )
    println( "Serial #:     $(string(devinfo.serial))" )
    println( "USB Bust:     $(int(devinfo.usb_bus))" )
    println( "USB Address:  $(int(devinfo.usb_addr))" )
    println( "Instance:     $(int(devinfo.instance))\n" )
end

function errorstring( status::Integer )
    p = ccall( ( :bladerf_strerror, libbladerf ), Ptr{Uint8}, (BladeRFStatus,), status )
    return bytestring( p )
end

function devices()
    devlist  = Array(Ptr{BladeRFDeviceInfo},1)
    ndevices = @bladerfcall( :bladerf_get_device_list, ( Ptr{BladeRFDeviceInfo}, ), pointer(devlist) )
    return pointer_to_array( devlist[1], (ndevices,), true )
end

function Base.open( devinfo::BladeRFDeviceInfo )
    handle = [C_NULL]
    @bladerfcall( :bladerf_open_with_devinfo, ( Ptr{Ptr{Void}}, Ptr{BladeRFDeviceInfo} ), handle, [devinfo] )
    return BladeRFDevice(handle[1])
end

function Base.open()
    devinfo = devices()[1]
    device  = open( devinfo )
    return device
end

function Base.close( device::BladeRFDevice )
    if device.handle != C_NULL
        @bladerfcall( :bladerf_close, ( Ptr{Void}, ), device.handle )
    end
    device.handle = C_NULL
end

function enable( device::BladeRFDevice, rfmodule::BladeRFModule )
    @bladerfcall( :bladerf_enable_module, (Ptr{Void}, Cint, Bool), device.handle, rfmodule.mod, true )
    return nothing
end

function disable( device::BladeRFDevice, rfmodule::BladeRFModule )
    @bladerfcall( :bladerf_enable_module, (Ptr{Void}, Cint, Bool), device.handle, rfmodule.mod, false )
    return nothing
end

const BLADERF_FORMAT_SC16_Q11 = 0

function config( device::BladeRFDevice, rfmodule::BladeRFModule; buffersize=1024, transfers=4, buffers=transfers+1, timeout=1000 )
    @bladerfcall( :bladerf_sync_config, (Ptr{Void},Cint,Cint,Uint32,Uint32,Uint32,Uint32), device.handle, rfmodule.mod, BLADERF_FORMAT_SC16_Q11, buffers, buffersize, transfers, timeout )
    return nothing
end

function receive( device::BladeRFDevice, samples::Integer; timeout = 1000 )
    buffer = zeros( Complex{Int16}, samples )
    @bladerfcall( :bladerf_sync_rx, (Ptr{Void},Ptr{Complex{Int16}},Uint,Ptr{Void},Uint), device.handle, buffer, samples, [C_NULL], timeout )
    return buffer
end

# Get sample rate
function samplerate( device::BladeRFDevice, rfmodule::BladeRFModule )
    rate = Array(Cuint,1)
    @bladerfcall( :bladerf_get_sample_rate, (Ptr{Void},Cint,Ptr{Cuint}), device.handle, rfmodule.mod, rate )
    int(rate[1])
end

# Set sample rate
function samplerate( device::BladeRFDevice, rfmodule::BladeRFModule, rate::Integer )
    actualrate = Array(Cuint,1)
    @bladerfcall( :bladerf_set_sample_rate, (Ptr{Void},Cint,Cuint,Ptr{Cuint}), device.handle, rfmodule.mod, rate, actualrate )
    int(actualrate[1])
end

end # module


text:
include("solver_preamble.jl")
import SeDuMi
factory = with_optimizer(SeDuMi.Optimizer, fid=0)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config,
                  ["quadratic_infeasible_lyapunov_switched_system",
                   "quadratic_feasible_lyapunov_switched_system"])
end


text:
# Warning: long computations ahead!
# Run this script in parallel on a high-performance computing server
using SimpleBirthDeathProcess;
using DelimitedFiles;
using Statistics;

function run_simulations(N, t, S, parameter_set)
  tot_sim = 5 * N
  τ = t / S

  for (id, i, η) = parameter_set
    state = zeros(Int64, S + 1, tot_sim)

    n = 1
    while n <= tot_sim
      y = rand_discrete(i, S, τ, η)

      if y.state[2] > 0
        copyto!(state, 1 + (n - 1) * (S + 1), y.state, 1, S + 1)
        n += 1
      end
    end

    file_name = string("simulation_", lpad(id, 2, '0'), ".csv")
    fp = open(file_name, "w")
    write(fp, string("# i = ", i, ", t = ", t, ", λ = ", η[1], ", μ = ", η[2],
                     ", θ = ", η[1] - η[2], "\n"))
    writedlm(fp, state, ',')
    close(fp)
  end

  nothing
end

function find_mle(N, t, parameter_set)
  λ = zeros(Float64, 4, N)
  μ = zeros(Float64, 4, N)

  for u = [1; 3; 5]
    for (id, i, η) = parameter_set
      fill!(λ, zero(Float64))
      fill!(μ, zero(Float64))

      file_name = string("simulation_", lpad(id, 2, '0'), ".csv")
      state = readdlm(file_name, ',', Int64, '\n', comments=true)
      x = zeros(Int64, size(state, 1), u)

      for n = 1:N
        lower_bound = 1 + (n - 1) * u
        upper_bound = lower_bound + u - 1
        copyto!(x, state[:, lower_bound:upper_bound])

        for (iter, idx) = enumerate((1:8:9, 1:4:9, 1:2:9, 1:1:9))
          k = length(idx)
          S = k - 1
          τ = t / S
          z = x[idx, :]
          extinct = findall(z[2, :] .== 0)

          for q = extinct
            o = zeros(Int64, k)

            while o[2] == 0
              copyto!(o, rand_discrete(i, S, τ, η).state)
            end

            copyto!(z, 1 + (q - 1) * k, o, 1, k)
          end

          y = observation_discrete_time_equal(τ, z)
          ρ = mle(y)

          λ[iter, n] = ρ[1]
          μ[iter, n] = ρ[2]
        end
      end

      for (iter, S) = enumerate([1; 2; 4; 8])
        str = string("mle_", lpad(id, 2, '0'), "_i_", lpad(i, 4, '0'),
                     "_u_", u, "_S_", S, ".csv")
        fp = open(str, "w")
        write(fp, "lambda,mu\n")
        writedlm(fp, hcat(λ[iter, :], μ[iter, :]), ',')
        close(fp)
      end
    end
  end

  nothing
end

function evaluate_estimates(parameter_set)
  fp = open("simulation_results.csv", "w")
  write(fp, string("u,i,S,lambda,bias_lambda,rmse_lambda,mu,bias_mu,rmse_mu,",
                   "theta,bias_theta,rmse_theta\n"))

  for u = [1; 3; 5]
    for (id, i, η) = parameter_set
      θ = η[1] - η[2]

      for (iter, S) = enumerate([1; 2; 4; 8])
        file_name = string("mle_", lpad(id, 2, '0'), "_i_", lpad(i, 4, '0'),
                           "_u_", u, "_S_", S, ".csv")
        params = readdlm(file_name, ',', Float64, '\n', header=true)

        N = size(params[1], 1)
        λ = params[1][:, 1]
        μ = params[1][:, 2]

        bias = λ .- η[1]
        bias_lambda = mean(bias)
        mse_lambda = mean(bias.^2)

        bias = μ .- η[2]
        bias_mu = mean(bias)
        mse_mu = mean(bias.^2)

        bias = (λ .- μ) .- θ
        bias_theta = mean(bias)
        mse_theta = mean(bias.^2)

        str = string(u, ",", i, ",", S, ",", η[1], ",", bias_lambda, ",",
                     sqrt(mse_lambda), ",", η[2], ",", bias_mu, ",",
                     sqrt(mse_mu), ",", θ, ",", bias_theta, ",",
                     sqrt(mse_theta), "\n")
        write(fp, str)
      end
    end
  end

  close(fp)

  nothing
end

N = 1_000_000;
t = 10;
S = 8;

brg_i_0010_1 = (
# E[N(t)] = 2 * i, SD[N(t)] = 1.25 * i
(1, 10, [0.30541797643422590; 0.236103258378231370]),
# E[N(t)] = 2 * i, SD[N(t)] = 1.5 * i
(2, 10, [0.42455264809296650; 0.355237930036971960]),
# E[N(t)] = 2 * i, SD[N(t)] = 2.0 * i
(3, 10, [0.72780453958794260; 0.658489821531948000]));

brg_i_0010_2 = (
# E[N(t)] = 4 * i, SD[N(t)] = 1.25 * i
(4, 10, [0.15956825719140408; 0.020938821079415016]),
# E[N(t)] = 4 * i, SD[N(t)] = 1.5 * i
(5, 10, [0.19927981441098427; 0.060650378298995215]),
# E[N(t)] = 4 * i, SD[N(t)] = 2.0 * i
(6, 10, [0.30036377824264300; 0.161734342130653900]));

brg_i_0100_1 = (
# E[N(t)] = 2 * i, SD[N(t)] = 1.25 * i
(7, 100, [2.7422635330902840; 2.6729488150342893]),
# E[N(t)] = 2 * i, SD[N(t)] = 1.5 * i
(8, 100, [3.9336102496776895; 3.864295531621695]),
# E[N(t)] = 2 * i, SD[N(t)] = 2.0 * i
(9, 100, [6.9661291646274510; 6.896814446571456]));

brg_i_0100_2 = (
# E[N(t)] = 4 * i, SD[N(t)] = 1.25 * i
(10, 100, [0.97185010941009; 0.8332206732981009]),
# E[N(t)] = 4 * i, SD[N(t)] = 1.5 * i
(11, 100, [1.368965681605892; 1.2303362454939029]),
# E[N(t)] = 4 * i, SD[N(t)] = 2.0 * i
(12, 100, [2.379805319922479; 2.24117588381049]));

brg_i_1000_1 = (
# E[N(t)] = 2 * i, SD[N(t)] = 1.25 * i
(13, 1000, [27.110719099650860; 27.041404381594866]),
# E[N(t)] = 2 * i, SD[N(t)] = 1.5 * i
(14, 1000, [39.024186265524920; 38.954871547468926]),
# E[N(t)] = 2 * i, SD[N(t)] = 2.0 * i
(15, 1000, [69.349375415022520; 69.280060696966530]));

brg_i_1000_2 = (
# E[N(t)] = 4 * i, SD[N(t)] = 1.25 * i
(16, 1000, [9.0946686315969500; 8.9560391954849600]),
# E[N(t)] = 4 * i, SD[N(t)] = 1.5 * i
(17, 1000, [13.065824353554970; 12.927194917442980]),
# E[N(t)] = 4 * i, SD[N(t)] = 2.0 * i
(18, 1000, [23.174220736720837; 23.035591300608850]));

drg_i_0010_1 = (
# E[N(t)] = i / 2, SD[N(t)] = i / 4
(19, 10, [0.0519860385419959; 0.12130075659799043]),
# E[N(t)] = i / 2, SD[N(t)] = i / 2
(20, 10, [0.3119162312519754; 0.38123094930796990]),
# E[N(t)] = i / 2, SD[N(t)] = i
(21, 10, [1.3516370020918933; 1.42095172014788780]));

drg_i_0010_2 = (
# E[N(t)] = i / 4, SD[N(t)] = i / 4
(22, 10, [0.1617343421306539; 0.30036377824264300]),
# E[N(t)] = i / 4, SD[N(t)] = i / 2
(23, 10, [0.8548815226905992; 0.99351095880258830]),
# E[N(t)] = i / 4, SD[N(t)] = i
(24, 10, [3.6274702449303806; 3.76609968104236970]));

drg_i_0100_1 = (
# E[N(t)] = i / 2, SD[N(t)] = i / 4
(25, 100, [0.8317766166719344; 0.9010913347279289]),
# E[N(t)] = i / 2, SD[N(t)] = i / 2
(26, 100, [3.4310785437717293; 3.5003932618277240]),
# E[N(t)] = i / 2, SD[N(t)] = i
(27, 100, [13.828286252170908; 13.897600970226904]));

drg_i_0100_2 = (
# E[N(t)] = i / 4, SD[N(t)] = i / 4
(28, 100, [2.2411758838104900; 2.3798053199224790]),
# E[N(t)] = i / 4, SD[N(t)] = i / 2
(29, 100, [9.1726476894099440; 9.3112771255219330]),
# E[N(t)] = i / 4, SD[N(t)] = i
(30, 100, [36.898534911807760; 37.037164347919740]));

drg_i_1000_1 = (
# E[N(t)] = i / 2, SD[N(t)] = i / 4
(31, 1000, [8.6296823979713190; 8.6989971160273140]),
# E[N(t)] = i / 2, SD[N(t)] = i / 2
(32, 1000, [34.622701668969270; 34.692016387025260]),
# E[N(t)] = i / 2, SD[N(t)] = i
(33, 1000, [138.59477875296108; 138.66409347101705]));

drg_i_1000_2 = (
# E[N(t)] = i / 4, SD[N(t)] = i / 4
(34, 1000, [23.035591300608850; 23.1742207367208370]),
# E[N(t)] = i / 4, SD[N(t)] = i / 2
(35, 1000, [92.350309356603380; 92.488938792715370]),
# E[N(t)] = i / 4, SD[N(t)] = i
(36, 1000, [369.60918158058150; 369.74781101669350]));


text:
mutable struct BagNode{T}
    @atomic next::Union{BagNode{T},Nothing}
    value::T

    BagNode{T}() where {T} = new{T}(nothing)
    BagNode{T}(next::Union{BagNode{T},Nothing}, value::T) where {T} = new{T}(next, value)
end

# TODO: find a proper/better bag algorithm?
mutable struct Bag{T}
    @atomic next::Union{BagNode{T},Nothing}
end

Bag() = Bag{Any}()
Bag{T}() where {T} = Bag{T}(nothing)

Base.eltype(::Type{Bag{T}}) where {T} = T
Base.IteratorSize(::Type{<:Bag}) = Base.SizeUnknown()

function isdeleted end

function Base.iterate(bag::Bag, (prev, curr) = (bag, @atomic bag.next))
    curr === nothing && return nothing
    next = @atomic curr.next
    while isdeleted(curr.value)
        if curr.value isa Message
            @trace(
                label = :cleaning_offer,
                offerid = offerid(curr.value.offer),
                taskid = objectid(current_task()),
                offer = curr.value.offer,
                msgs = bag,
            )
        end
        # `prev` may be phisically removed while `curr` is phisically removed.
        # But this is OK since `curr` is already logically removed.
        # TODO: check this
        curr, ok = @atomicreplace(prev.next, curr => next)
        if ok
            curr = next
        end
        curr === nothing && return nothing
        next = @atomic curr.next
    end
    return (curr.value, (curr, next))
end

function Base.push!(bag::Bag{T}, v) where {T}
    v = convert(T, v)
    next = @atomic bag.next
    node = BagNode{T}(next, v)
    while true
        next, success = @atomicreplace(bag.next, next => node)
        success && return bag
        node = BagNode{T}(next, v)
    end
end

function Base.in(v::T, bag::Bag{T}) where {T}
    for x in bag
        x === v && return true
    end
    return false
end

function nitems(bag::Bag)
    n = 0
    for _ in bag
        n += 1
    end
    return n
end

function tryremove!(bag::Bag{T}, item::T) where {T}
    prev = bag
    y = iterate(bag)
    while true
        y === nothing && return false
        value, (curr, next) = y
        if value === item
            @atomicreplace(prev.next, curr => next)
            return true
        end
        prev = curr
        y = iterate(bag, (curr, next))
    end
end


text:
# a < b < c, a+b+c=1000
for a = 1:333 # 3a < a+b+c < 1000, a < 1000/3
	for c = 334:1000 # 3c > a+b+c = 1000, c > 1000/3
		b = 1000 -a -c
		if b <= 0 continue end
		if a^2 + b^2 == c^2
			println(a, '\t', b, '\t', c)
		end
	end
end

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
#######
# Universal Kalman filtering methods

function predict!(kf::KalmanFilter)
    kf.x = ap(kf.f,kf.x)
    kf
end

function predict(kf::KalmanFilter)
    predict!(copy(kf))
end

function update(kf::KalmanFilter,y::Observation)
    update!(copy(kf),y)
end

function update!(kf::KalmanFilter,y::Observation)
    (res,ph,s) = covs(kf,y)
    xn = kf.x.x + ph * (s\res)
    pn = kf.x.p - ph * (s'\ph')

    # This is an ugly hack which works for now
    if typeof(kf.x) <: AbstractUnscentedState
        kf.x = UnscentedState(xn,pn,kf.x.α,kf.x.β,kf.x.κ)
    else
        kf.x = State(xn,pn)
    end
    kf
end

function predictupdate(kf::KalmanFilter,y::Observation)
    update(predict(kf),y)
end

function predictupdate!(kf::KalmanFilter,y::Observation)
    update!(predict!(kf),y)
end