text:
import Base:writemime

function html_body(app)
    """
        <title>WebRTC.jl</title>

        <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
        <meta charset="utf-8">
        <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no">
        <link rel="author" type="text/html" href="https://plus.google.com/+MuazKhan">
        <meta name="author" content="Muaz Khan">
        <meta http-equiv="X-UA-Compatible" content="IE=edge,chrome=1">

        <meta name="title" content="Video Conferencing" />
        <meta name="description" content="Video Conferencing using WebRTC." />
        <meta name="keywords" content="WebRTC, Video Conferencing, Demo, Example, Experiment" />

        <link rel="stylesheet" href="https://cdn.webrtc-experiment.com/style.css">

                        WebRTC.jl

                        <input type = "text" id = "conference-name">
                        <input type = "button" id = "setup-new-room" class = "setup" value="Setup Conference" style="float: right;">

    """
end

function Base.writemime(io::IO, ::MIME"text/html", app::WebRTCApplication)
    print(io, html_body(app))
end


text:
##
#
# Generate all non-zero roots of polynomials with coefficients 1 or -1 of
# degree at most 25.
#
# Warning: this may take a while (a few hours in my case) and generate almost
#          50 GB file (3 221 224 897 roots).
#

using ChristensenSet

polynomials = PolynomialIterator([-1.0, 1.0], 25)

@time find_roots!(polynomials, "ones_Float64.dat")

ext: jl


text:
module SlurpTests
using  ExpressionPatterns.Destructuring
using  ExpressionPatterns.Analyzer
using  ExpressionPatterns.PatternStructure.SlurpTypes
using  Base.Test

# generic greedy

@test isa(analyze(:(*{x+y,z},1+1,1,*{b})).child.children[1].head, GenericGreedySlurp)

@macrods t1(*{x+y,z},1+1,1,*{b}) (x,y,z)

@test @t1(10+2, a, 30+54, b, 53+1, c, 1+1, 1, 1+1, 1) == ([10,30,53,1], [2,54,1,1], [:a,:b,:c,1])

# generic lazy

@test isa(analyze(:(:?{x+y,z},1+1,1,*{b})).child.children[1].head, GenericLazySlurp)

@macrods t2(:?{x+y,z},1+1,1,*{b}) (x,y,z)

@test @t2(10+2, a, 30+54, b, 53+1, c, 1+1, 1, 1+1, 1) == ([10,30,53], [2,54,1], [:a,:b,:c])

# last, simple

@test isa(analyze(:(a,b,*{c},d,e)).child.children[3].head, SimpleLastSlurp)

@macrods t3(a,b,*{c},d,e) c

@test @t3(1,2,3,4,5,6,7) == [3,4,5]

# simple, until, greedy

@test isa(analyze(:(*{a},3,*{b})).child.children[1].head, SimpleGreedySlurpUntil)

@macrods t4(*{a},3,*{b}) a

@test @t4(0,1,2,3,4,5,6) == [0,1,2]

# simple, until, lazy

@test isa(analyze(:(:?{a},3,*{b})).child.children[1].head, SimpleLazySlurpUntil)

@macrods t5(:?{a},3,*{b}) a

@test @t5(0,1,2,3,4,5,6) == [0,1,2]

end


text:
using Random

@testset "Profiles" begin

σ = rand()
g = Gaussian(σ,1)
@test g(σ) ≈ exp(-1)/sqrt(π*σ^2)

@test radius(g) == σ
@test strength(g) == 1

σ = 0.2
A = 1
dg = DGaussian(σ,A)
x = 0.1
@test dg(x) == -2*A*x/sqrt(π)/σ^3*exp(-x^2/σ^2)

end

@testset "SpatialFields" begin

σ = 0.2
A = 1

g = EmptySpatialField()
@test g(rand(2)...) == 0.0

g = SpatialGaussian(σ,σ,0.0,0.5,1)
@test g(0,0.5+σ) ≈ g(σ,0.5) ≈ g(-σ,0.5) ≈ g(0,0.5-σ) ≈ exp(-1)/(π*σ^2)

u = 1
v = 0
gc = SpatialGaussian(σ,σ,0.0,0.5,1,u,v)
t = 1
@test gc(u*t,0.5+σ+v*t,t) ≈ gc(σ+u*t,0.5+v*t,t) ≈ gc(-σ+u*t,0.5+v*t,t) ≈ gc(0+u*t,0.5-σ+v*t,t) ≈ exp(-1)/(π*σ^2)

g = EmptySpatialField()
for x in [-0.5,0,0.5], y in [-0.5,0,0.5]
  g += SpatialGaussian(0.2,0.5,x,y,1)
end
@test g(0.5,0) ≈ 5.5357580598289235

g2 = -g
@test g2(0.5,0) ≈ -5.5357580598289235

gr = PhysicalGrid((-3.0,3.0),(-2.0,2.0),0.03)
w = Nodes(Dual,size(gr))
xg, yg = coordinates(w,gr)

gfield = GeneratedField(w,g,gr)

@test gfield()[104,70] == g(xg[104],yg[70])
@test datatype(gfield) == typeof(w)

wfield = PulseField(gfield,0.5,0.1)

@test maximum(wfield(0.5)) ≈ maximum(gfield()) ≈ 5.532960088678624
@test maximum(wfield(2)) ≈ 0.0
@test datatype(wfield) == datatype(gfield) <: ScalarGridData

q = Edges(Primal,size(gr))
gauss = SpatialGaussian(0.5,0,0,1)
gaussfield = GeneratedField(q,gauss,EmptySpatialField(),gr)
q .= gaussfield()

ffield = PulseField(gaussfield,0.5,0.1)

@test maximum(ffield(0.5).u) ≈ maximum(q.u) ≈ 1.272094144652253
@test maximum(ffield(2).u) ≈ 0.0
@test maximum(ffield(0.5).v) ≈ maximum(q.v) ≈ 0.0
@test maximum(ffield(2).v) ≈ maximum(q.v) ≈ 0.0

@test datatype(ffield) == datatype(gaussfield) <: VectorGridData

gq = EdgeGradient(Dual,size(gr))
tfield = GeneratedField(gq,[g,gauss,EmptySpatialField(),EmptySpatialField()],gr)
gq .= tfield()
@test maximum(gq.dudx) ≈ 5.532960088678624
@test maximum(gq.dudy) ≈ 1.2732395447351625
@test maximum(gq.dvdx) == maximum(gq.dvdy) == 0.0

ffield = PulseField(tfield,0.5,0.1)
@test datatype(ffield) <: TensorGridData

end

@testset "Derivatives of gaussians" begin

  σx = 0.5
  σy = 0.5
  x0 = 0
  y0 = 0
  A = 1
  dgaussx = SpatialGaussian(σx,σy,x0,y0,A,deriv=1)
  dgaussy = SpatialGaussian(σx,σy,x0,y0,A,deriv=2)
  x = 0.1
  y = 0.2
  @test dgaussx(0.1,0.2) == dgaussy(0.2,0.1) ≈ -2*A*x/π/σx^3/σy*exp(-x^2/σx^2)*exp(-y^2/σy^2)

end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
#! format: off
var_map(s::String) = var_map(Val(Symbol(s)))
var_map(::Val{T}) where {T} = nothing

var_map(::Val{Symbol("ρ")}) = ("rho", ())
var_map(::Val{Symbol("ρu[1]")}) = ("u_mean", (:ρ,))
var_map(::Val{Symbol("ρu[2]")}) = ("v_mean", (:ρ,))
var_map(::Val{Symbol("moisture.ρq_tot")}) = ("qt_mean", (:ρ,))
var_map(::Val{Symbol("turbconv.updraft[1].ρa")}) = ("updraft_fraction", (:ρ,))
var_map(::Val{Symbol("turbconv.updraft[1].ρaw")}) = ("updraft_w", (:ρ, :a))
var_map(::Val{Symbol("turbconv.updraft[1].ρaq_tot")}) = ("updraft_qt", (:ρ, :a))
var_map(::Val{Symbol("turbconv.updraft[1].ρaθ_liq")}) = ("updraft_thetali", (:ρ, :a))
var_map(::Val{Symbol("turbconv.environment.ρatke")}) = ("tke_mean", (:ρ, :a))
var_map(::Val{Symbol("turbconv.environment.ρaθ_liq_cv")}) = ("env_thetali2", (:ρ, :a))
var_map(::Val{Symbol("turbconv.environment.ρaq_tot_cv")}) = ("env_qt2", (:ρ, :a))
#! format: on

max_forks_repo_name: jm-c/CLIMA


text:
using Quaternions
using LinearAlgebra
using Random
using Test

@testset "Octonion" begin
    @testset "type aliases" begin
        @test OctonionF16 === Octonion{Float16}
        @test OctonionF32 === Octonion{Float32}
        @test OctonionF64 === Octonion{Float64}
    end

    @testset "Constructors" begin
        @testset "from coefficients" begin
            cs = [(1, 2.0, 3.0f0, 4//1, 5, 6, 7, 8), (1//1, 2.0f0, 3.0f0, 4, 5, 6, 7, 8)]
            @testset for coef in cs, T in (Float32, Float64, Int), norm in (true, false)
                q = @inferred Octonion{T}(coef..., norm)
                @test q isa Octonion{T}
                @test q.norm === norm
                @test q === Octonion{T}(convert.(T, coef)..., norm)
                q2 = @inferred Octonion(convert.(T, coef)..., norm)
                @test Octonion(convert.(T, coef)..., norm) === q
                if !norm
                    @test Octonion(convert.(T, coef)...) === q
                end
            end
        end
        @testset "from real" begin
            @testset for x in (-1//1, 1.0, 2.0), T in (Float32, Float64, Int, Rational{Int})
                coef = T.((x, zeros(7)...))
                @test @inferred(Octonion{T}(x)) === Octonion{T}(coef..., isone(abs(x)))
                @test @inferred(Octonion(T(x))) === Octonion{T}(coef..., isone(abs(x)))
            end
        end
        @testset "from complex" begin
            @testset for z in (1 + 0im, -im, 1 + 2im),
                T in (Float32, Float64, Int, Rational{Int})

                coef = T.((reim(z)..., zeros(6)...))
                z2 = Complex{T}(z)
                norm = isone(abs(z))
                @test Octonion{T}(z) === Octonion{T}(coef..., norm)
                @test @inferred(Octonion(z2)) === Octonion{T}(coef..., norm)
            end
        end
        @testset "from quaternion" begin
            qs = (Quaternion(1, 2, 3, 4), QuaternionF64(0, 1, 0, 0, true))
            @testset for q in qs, T in (Float32, Float64)
                coef = T.((q.s, q.v1, q.v2, q.v3, zeros(4)...))
                q2 = Quaternion{T}(q)
                @test @inferred(Octonion{T}(q)) === Octonion{T}(coef..., q.norm)
                @test @inferred(Octonion(q2)) === Octonion{T}(coef..., q.norm)
            end
        end
        @testset "from octonion" begin
            os = (
                Octonion(1, 2, 3, 4, 5, 6, 7, 8), OctonionF64(0, 1, 0, 0, 0, 0, 0, 0, true)
            )
            @testset for o in os, T in (Float32, Float64)
                coef = T.((o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7))
                @test @inferred(Octonion{T}(o)) === Octonion{T}(coef..., o.norm)
                @test @inferred(Octonion(o)) === o
            end
        end
        @testset "from vector" begin
            s = randn()
            v = randn(7)
            @test @inferred(Octonion(s, v)) === Octonion(s, v...)
            @test @inferred(Octonion(v)) === Octonion(0, v)
        end
    end

    @testset "==" begin
        @test Octonion(1.0, 2, 3, 4, 5, 6, 7, 8) == Octonion(1, 2, 3, 4, 5, 6, 7, 8)
        @test Octonion(1.0, 2, 3, 4, 5, 6, 7, 8) != Octonion(1, 2, 3, 4, 1, 2, 3, 4)
        @test Octonion(1, 0, 0, 0, 0, 0, 0, 0, false) ==
            Octonion(1, 0, 0, 0, 0, 0, 0, 0, true) # test that .norm field does not affect equality
    end

    @testset "convert" begin
        @test convert(Octonion{Float64}, 1) === Octonion(1.0)
        @test convert(Octonion{Float64}, Complex(1, 2)) ===
            Octonion(1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
        @test convert(Octonion{Float64}, Quaternion(1, 2, 3, 4)) ===
            Octonion(1.0, 2.0, 3.0, 4.0, 0.0, 0.0, 0.0, 0.0)
        @test convert(Octonion{Float64}, Quaternion(0, 1, 0, 0, true)) ===
            Octonion(0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, true)
        @test convert(Octonion{Float64}, Octonion(1, 2, 3, 4, 5, 6, 7, 8)) ===
            Octonion(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
        @test convert(Octonion{Float64}, Octonion(0, 0, 0, 0, 1, 0, 0, 0, true)) ===
            Octonion(0.0, 0.0, 0.0, 0.0, 1, 0, 0, 0, true)
    end

    @testset "promote" begin
        @test promote(Octonion(1.0, 2:8...), 1.0) === (Octonion(1.0, 2:8...), Octonion(1.0))
        @test promote(Octonion(1.0f0, 2:8...), 2.0) ===
            (Octonion(1.0, 2:8...), Octonion(2.0))
        @test promote(Octonion(1.0f0), 2 + 3im) ===
            (Octonion(1.0f0), Octonion(2.0f0 + 3.0f0im))
        @test promote(Octonion(1.0f0), Quaternion(1, 2, 3, 4)) ===
            (Octonion(1.0f0), Octonion(1.0f0:4.0f0..., fill(0, 4)...))
        @test promote(Octonion(1.0f0), Octonion(2.0)) === (Octonion(1.0), Octonion(2.0))

        @test Octonion(1) == 1.0
        @test Octonion(1, 2, fill(0, 6)...) == Complex(1.0, 2.0)
        @test Octonion(1) == 1.0
        @test Octonion(1:4..., fill(0, 4)...) == Quaternion(1.0:4.0...)
    end

    @testset "shorthands" begin
        @test octo(1) === Octonion(1) # checking the .norm field in particular
        @test octo(1, 0, 0, 0, 0, 0, 0, 0) === Octonion(1, 0, 0, 0, 0, 0, 0, 0) # checking the .norm field in particular
        @test octo(1, 2, 3, 4, 5, 6, 7, 8) === Octonion(1, 2, 3, 4, 5, 6, 7, 8)
        @test octo(Octonion(1, 0, 0, 0, 0, 0, 0, 0)) === Octonion(1, 0, 0, 0, 0, 0, 0, 0) # checking the .norm field in particular
        @test octo(Octonion(1, 2, 3, 4, 5, 6, 7, 8)) === Octonion(1, 2, 3, 4, 5, 6, 7, 8)
        @test octo(1, 0, 0, 0, 0, 0, 0, 0, false).norm == false # respect the .norm input (even if wrong)
        @test octo(1, 2, 3, 4, 5, 6, 7, 8, true).norm == true # respect the .norm input (even if wrong)
        @test octo(1, collect(2:8)) === Octonion(1:8...)
        @test octo(collect(2:8)) === Octonion(0, 2:8...)
    end

    @testset "random generation" begin
        @testset "octorand" begin
            o = octorand()
            @test o isa Octonion
            @test !o.norm
        end

        @testset "rand($H)" for H in (OctonionF32, OctonionF64)
            rng = Random.MersenneTwister(42)
            o = rand(rng, H)
            @test o isa H
            @test !o.norm

            os = rand(rng, H, 1000)
            @test eltype(os) === H
            @test length(os) == 1000
            xs = map(os) do o
                return [real(o); imag_part(o)...]
            end
            xs_mean = sum(xs) / length(xs)
            xs_var = sum(x -> abs2.(x .- xs_mean), xs) / (length(xs) - 1)
            @test all(isapprox.(xs_mean, 0.5; atol=0.1))
            @test all(isapprox.(xs_var, 1 / 12; atol=0.01))
        end

        @testset "randn($H)" for H in (OctonionF32, OctonionF64)
            rng = Random.MersenneTwister(42)
            o = randn(rng, H)
            @test o isa H
            @test !o.norm

            os = randn(rng, H, 10000)
            @test eltype(os) === H
            @test length(os) == 10000
            xs = map(os) do o
                return [real(o); imag_part(o)...]
            end
            xs_mean = sum(xs) / length(xs)
            xs_var = sum(x -> abs2.(x .- xs_mean), xs) / (length(xs) - 1)
            @test all(isapprox.(xs_mean, 0; atol=0.1))
            @test all(isapprox.(xs_var, 1 / 8; atol=0.1))
        end
    end

    @testset "basic" begin
        q = randn(OctonionF64)
        qnorm = normalize(q)
        @test real(q) === q.s
        @test_throws MethodError imag(q)
        @test @test_deprecated(Quaternions.imag(q)) == [q.v1, q.v2, q.v3, q.v4, q.v5, q.v6, q.v7]
        @test imag_part(q) === (q.v1, q.v2, q.v3, q.v4, q.v5, q.v6, q.v7)
        @test conj(q) ===
            Octonion(q.s, -q.v1, -q.v2, -q.v3, -q.v4, -q.v5, -q.v6, -q.v7, q.norm)
        @test conj(qnorm) === Octonion(
            qnorm.s,
            -qnorm.v1,
            -qnorm.v2,
            -qnorm.v3,
            -qnorm.v4,
            -qnorm.v5,
            -qnorm.v6,
            -qnorm.v7,
            qnorm.norm,
        )
        @test conj(conj(q)) === q
        @test conj(conj(qnorm)) === qnorm
        @test float(Octonion(1:8...)) === Octonion(1.0:8.0...)
        @test Quaternions.abs_imag(q) ==
            abs(Octonion(0, q.v1, q.v2, q.v3, q.v4, q.v5, q.v6, q.v7))
    end

    @testset "algebraic properties" begin
        for _ in 1:100, T in (Float32, Float64, Int32, Int64)
            if T <: Integer
                q, q1, q2, q3 = [octo(rand((-T(100)):T(100), 8)...) for _ in 1:4]
                c1, c2 = [complex(rand((-T(100)):T(100), 2)...) for _ in 1:2]
            else
                q, q1, q2, q3 = randn(Octonion{T}, 4)
                c1, c2 = randn(Complex{T}, 2)
            end

            # skewfield
            test_group(q1, q2, q3, +, zero(q), -)
            # test all group properties but associativity
            test_neutral(q1, one(q1), *)
            test_inverse(q1, one(q1), *, inv)
            test_multiplicative(q1, q2, *, norm)

            # complex embedding
            test_multiplicative(c1, c2, *, Octonion)
            test_multiplicative(c1, c2, +, Octonion)
        end
    end

    @testset "isreal" begin
        @test isreal(octo(1))
        @test !isreal(octo(1, 1, 0, 0, 0, 0, 0, 0))
        @test !isreal(octo(1, 0, 1, 0, 0, 0, 0, 0))
        @test !isreal(octo(1, 0, 0, 1, 0, 0, 0, 0))
        @test !isreal(octo(1, 0, 0, 0, 1, 0, 0, 0))
        @test !isreal(octo(1, 0, 0, 0, 0, 1, 0, 0))
        @test !isreal(octo(1, 0, 0, 0, 0, 0, 1, 0))
        @test !isreal(octo(1, 0, 0, 0, 0, 0, 0, 1))
    end

    @testset "iszero" begin
        @test iszero(octo(0))
        @test !iszero(octo(1))
        @test !iszero(octo(0, 1, 0, 0, 0, 0, 0, 0))
        @test !iszero(octo(0, 0, 1, 0, 0, 0, 0, 0))
        @test !iszero(octo(0, 0, 0, 1, 0, 0, 0, 0))
        @test !iszero(octo(0, 0, 0, 0, 1, 0, 0, 0))
        @test !iszero(octo(0, 0, 0, 0, 0, 1, 0, 0))
        @test !iszero(octo(0, 0, 0, 0, 0, 0, 1, 0))
        @test !iszero(octo(0, 0, 0, 0, 0, 0, 0, 1))
        @test !iszero(octo(0, 0, 0, 0, 0, 0, 0, 0, true))
    end

    @testset "isone" begin
        @test isone(octo(1))
        @test !isone(octo(-1))
        @test !isone(octo(0, 1, 0, 0, 0, 0, 0, 0))
        @test !isone(octo(1, 1, 0, 0, 0, 0, 0, 0))
        @test !isone(octo(1, 0, 1, 0, 0, 0, 0, 0))
        @test !isone(octo(1, 0, 0, 1, 0, 0, 0, 0))
        @test !isone(octo(1, 0, 0, 0, 1, 0, 0, 0))
        @test !isone(octo(1, 0, 0, 0, 0, 1, 0, 0))
        @test !isone(octo(1, 0, 0, 0, 0, 0, 1, 0))
        @test !isone(octo(1, 0, 0, 0, 0, 0, 0, 1))
    end

    @testset "isfinite" begin
        @test isfinite(octo(1:8...))
        for val in (Inf, -Inf, NaN)
            @test !isfinite(octo(val, 0, 0, 0, 0, 0, 0, 0))
            @test !isfinite(octo(0, val, 0, 0, 0, 0, 0, 0))
            @test !isfinite(octo(0, 0, val, 0, 0, 0, 0, 0))
            @test !isfinite(octo(0, 0, 0, val, 0, 0, 0, 0))
            @test !isfinite(octo(0, 0, 0, 0, val, 0, 0, 0))
            @test !isfinite(octo(0, 0, 0, 0, 0, val, 0, 0))
            @test !isfinite(octo(0, 0, 0, 0, 0, 0, val, 0))
            @test !isfinite(octo(0, 0, 0, 0, 0, 0, 0, val))
            @test isfinite(octo(fill(val, 8)..., true))
        end
    end

    @testset "isinf" begin
        @test !isinf(octo(1, 2, 3, 4, 5, 6, 7, 8))
        @test !isinf(octo(1, 2, 3, 4, 5, 6, 7, NaN))
        for inf in (Inf, -Inf)
            @test isinf(octo(inf, 0, 0, 0, 0, 0, 0, 0))
            @test isinf(octo(0, inf, 0, 0, 0, 0, 0, 0))
            @test isinf(octo(0, 0, inf, 0, 0, 0, 0, 0))
            @test isinf(octo(0, 0, 0, inf, 0, 0, 0, 0))
            @test isinf(octo(0, 0, 0, 0, inf, 0, 0, 0))
            @test isinf(octo(0, 0, 0, 0, 0, inf, 0, 0))
            @test isinf(octo(0, 0, 0, 0, 0, 0, inf, 0))
            @test isinf(octo(0, 0, 0, 0, 0, 0, 0, inf))
            @test !isinf(octo(fill(inf, 8)..., true))
        end
    end

    @testset "isnan" begin
        @test !isnan(octo(1, 2, 3, 4, 5, 6, 7, 8))
        @test !isnan(octo(1, 2, 3, 4, 5, 6, 7, Inf))
        @test !isnan(octo(1, 2, 3, 4, 5, 6, 7, -Inf))
        @test isnan(octo(NaN, 2, 3, 4, 5, 6, 7, 8))
        @test isnan(octo(1, NaN, 3, 4, 5, 6, 7, 8))
        @test isnan(octo(1, 2, NaN, 4, 5, 6, 7, 8))
        @test isnan(octo(1, 2, 3, NaN, 5, 6, 7, 8))
        @test isnan(octo(1, 2, 3, 4, NaN, 6, 7, 8))
        @test isnan(octo(1, 2, 3, 4, 5, NaN, 7, 8))
        @test isnan(octo(1, 2, 3, 4, 5, 6, NaN, 8))
        @test isnan(octo(1, 2, 3, 4, 5, 6, 7, NaN))
    end

    @testset "*" begin
        # verify basic correctness
        q0 = Octonion(1,0,0,0,0,0,0,0)
        q1 = Octonion(0,1,0,0,0,0,0,0)
        q2 = Octonion(0,0,1,0,0,0,0,0)
        q3 = Octonion(0,0,0,1,0,0,0,0)
        q4 = Octonion(0,0,0,0,1,0,0,0)
        q5 = Octonion(0,0,0,0,0,1,0,0)
        q6 = Octonion(0,0,0,0,0,0,1,0)
        q7 = Octonion(0,0,0,0,0,0,0,1)
        @test q0 * q0 == q0
        @test q0 * q1 == q1
        @test q0 * q2 == q2
        @test q0 * q3 == q3
        @test q0 * q4 == q4
        @test q0 * q5 == q5
        @test q0 * q6 == q6
        @test q0 * q7 == q7
        @test q1 * q0 == q1
        @test q1 * q1 == -q0
        @test q1 * q2 == q3
        @test q1 * q3 == -q2
        @test q1 * q4 == -q7
        @test q1 * q5 == -q6
        @test q1 * q6 == q5
        @test q1 * q7 == q4
        @test q2 * q0 == q2
        @test q2 * q1 == -q3
        @test q2 * q2 == -q0
        @test q2 * q3 == q1
        @test q2 * q4 == q6
        @test q2 * q5 == -q7
        @test q2 * q6 == -q4
        @test q2 * q7 == q5
        @test q3 * q0 == q3
        @test q3 * q1 == q2
        @test q3 * q2 == -q1
        @test q3 * q3 == -q0
        @test q3 * q4 == -q5
        @test q3 * q5 == q4
        @test q3 * q6 == -q7
        @test q3 * q7 == q6
        @test q4 * q0 == q4
        @test q4 * q1 == q7
        @test q4 * q2 == -q6
        @test q4 * q3 == q5
        @test q4 * q4 == -q0
        @test q4 * q5 == -q3
        @test q4 * q6 == q2
        @test q4 * q7 == -q1
        @test q5 * q0 == q5
        @test q5 * q1 == q6
        @test q5 * q2 == q7
        @test q5 * q3 == -q4
        @test q5 * q4 == q3
        @test q5 * q5 == -q0
        @test q5 * q6 == -q1
        @test q5 * q7 == -q2
        @test q6 * q0 == q6
        @test q6 * q1 == -q5
        @test q6 * q2 == q4
        @test q6 * q3 == q7
        @test q6 * q4 == -q2
        @test q6 * q5 == q1
        @test q6 * q6 == -q0
        @test q6 * q7 == -q3
        @test q7 * q0 == q7
        @test q7 * q1 == -q4
        @test q7 * q2 == -q5
        @test q7 * q3 == -q6
        @test q7 * q4 == q1
        @test q7 * q5 == q2
        @test q7 * q6 == q3
        @test q7 * q7 == -q0

        @testset "* same between Quaternions and Octonion" begin
            # make sure this tracks the `*` tests for Quaternions
            q1 = Octonion(1,0,0,0,0,0,0,0)
            qi = Octonion(0,1,0,0,0,0,0,0)
            qj = Octonion(0,0,1,0,0,0,0,0)
            qk = Octonion(0,0,0,1,0,0,0,0)
            @test q1 * q1 == q1
            @test q1 * qi == qi
            @test q1 * qj == qj
            @test q1 * qk == qk
            @test qi * q1 == qi
            @test qi * qi == -q1
            @test qi * qj == qk
            @test qi * qk == -qj
            @test qj * q1 == qj
            @test qj * qi == -qk
            @test qj * qj == -q1
            @test qj * qk == qi
            @test qk * q1 == qk
            @test qk * qi == qj
            @test qk * qj == -qi
            @test qk * qk == -q1
        end
    end

    @testset "abs2" for _ in 1:100, T in (Float16, Float32, Float64)
        o = rand(Octonion{T})
        @test abs2(o) == o'*o
    end

    @testset "/" begin
        for _ in 1:100
            o, o2 = randn(OctonionF64, 2)
            x = randn()
            @test o / o ≈ o \ o ≈ one(o)
            @test o / o2 ≈ o * inv(o2)
            @test o2 \ o ≈ inv(o2) * o
            @test o / x ≈ x \ o ≈ inv(x) * o
        end
    end

    @testset "^" begin
        @testset "^(::Octonion, ::Real)" begin
            for _ in 1:100
                o = randn(OctonionF64)
                @test @inferred(o^2.0) ≈ o * o
                @test o^1.0 ≈ o
                @test o^-1.0 ≈ inv(o)
                @test o^1.3 ≈ exp(1.3 * log(o))
                @test o^7.8 ≈ exp(7.8 * log(o))
                @test o^1.3f0 ≈ exp(1.3f0 * log(o))
                @test o^7.8f0 ≈ exp(7.8f0 * log(o))
            end
        end
        @testset "^(::Octonion, ::Octonion)" begin
            @test octo(Float64(ℯ))^octo(0, 0, 0, 0, 0, 0, π / 2, 0) ≈
                octo(0, 0, 0, 0, 0, 0, 1, 0)
            z = (3.5 + 2.3im)^(0.2 + 1.7im)
            @test octo(3.5, 0, 0, 0, 0, 0, 2.3, 0)^octo(0.2, 0, 0, 0, 0, 0, 1.7, 0) ≈
                octo(real(z), 0, 0, 0, 0, 0, imag(z), 0)
            for _ in 1:100
                q, p = randn(OctonionF64, 2)
                @test @inferred(q^p) ≈ exp(p * log(q))
            end
        end
    end

    @testset "non-analytic functions" begin
        unary_funs = [conj, abs, abs2, norm, sign]
        # since every octonion is conjugate to a quaternion,
        # one can establish correctness as follows:
        @testset for fun in unary_funs
            for _ in 1:100
                o1, o2 = randn(OctonionF64, 2)
                q = randn(QuaternionF64)
                o = octo(q)
                @test @inferred(fun(o)) ≈ fun(q)
                @test o2 * fun(o1) * inv(o2) ≈ fun(o2 * o1 * inv(o2))
            end
        end
    end

    @testset "analytic functions" begin
        unary_funs = [sqrt, inv, exp, log]
        # since every octonion is conjugate to a quaternion,
        # one can establish correctness as follows:
        @testset for fun in unary_funs
            for _ in 1:100
                o1, o2 = randn(OctonionF64, 2)
                q = randn(QuaternionF64)
                o = octo(q)
                @test @inferred(fun(o)) ≈ fun(q)
                @test o2 * fun(o1) * inv(o2) ≈ fun(o2 * o1 * inv(o2))
            end
        end

        @testset "identities" begin
            for _ in 1:100
                o = randn(OctonionF64)
                @test inv(o) * o ≈ o * inv(o) ≈ one(o)
                @test sqrt(o) * sqrt(o) ≈ o
                @test exp(log(o)) ≈ o
                @test exp(zero(o)) === one(o)
                @test log(one(o)) === zero(o)
            end
            @test log(zero(OctonionF64)) === octo(-Inf)
        end
    end

    @testset "normalize" begin
        for _ in 1:100
            q = randn(OctonionF64)
            qnorm = @inferred normalize(q)
            @test abs(qnorm) ≈ 1
            @test qnorm.norm
            @test q ≈ abs(q) * qnorm
            @test normalize(qnorm) === qnorm
        end
        @test_broken @inferred(normalize(octo(1:8...)))
    end

    @testset "normalizea" begin
        for _ in 1:100
            q = randn(OctonionF64)
            qnorm, a = @inferred normalizea(q)
            @test abs(qnorm) ≈ 1
            @test qnorm.norm
            @test a isa Real
            @test a ≈ abs(q)
            @test q ≈ a * qnorm
            @test normalizea(qnorm) === (qnorm, one(real(q)))
        end
        @test_broken @inferred(normalizea(octo(1:8...)))
    end
end


text:
using AdventOfCode, Test, Mocking, HTTP, Dates
Mocking.activate()
@testset "All tests" begin
    AOC = AdventOfCode
    @testset "_is_unlocked" begin
        @test AOC._is_unlocked(2018, 1)
        @test (@test_logs (:warn, "Advent of Code for year 3000 and day 1 hasn't unlocked yet.") !AOC._is_unlocked(3000, 1))
    end
    @testset "_get_cookies" begin
        cleanup = false
        if "AOC_SESSION" ∈ keys(ENV)
            temp = pop!(ENV, "AOC_SESSION")
            cleanup = true
        end
        @test_throws ErrorException("Session cookie in ENV[\"AOC_SESSION\"] needed to download data.") AOC._get_cookies()
        ENV["AOC_SESSION"] = "COOKIE"
        @test AOC._get_cookies() == Dict("session" => "COOKIE")

        if cleanup
            ENV["AOC_SESSION"] = temp
        end
    end
    @testset "_base_url" begin
        @test AOC._base_url(2019, 1) == "https://adventofcode.com/2019/day/1"
    end
    @testset "_template" begin
        @test AOC._template(2019, 1) == """
        # https://adventofcode.com/2019/day/1
        using AdventOfCode

        input = readlines("data/2019/day_1.txt")

        function part_1(input)
            nothing
        end
        @info part_1(input)

        function part_2(input)
            nothing
        end
        @info part_2(input)
        """
    end
    get_patch = @patch HTTP.get(a...; kw...) = (status = 200, body = "TESTDATA")
    get_error_patch = @patch HTTP.get(a...; kw...) = (status = 404, body = "ERROR")
    @testset "_download_data" begin
        apply(get_patch) do
            @test AdventOfCode._download_data(2019, 1) == "TESTDATA"
        end
        apply(get_error_patch) do
            @test_throws ErrorException("Unable to download data") AdventOfCode._download_data(2019, 1)
        end
    end
    @testset "_setup_data_file" begin
        apply(get_patch) do
            try
                AdventOfCode._setup_data_file(2019, 1)
                @test isdir("data/2019")
                @test isfile("data/2019/day_1.txt")
                @test readlines("data/2019/day_1.txt") == ["TESTDATA"]
                @test_logs (:warn, r"will not redownload") AdventOfCode._setup_data_file(2019, 1)
            finally
                rm("data", recursive = true, force = true)
            end
        end
    end
    @testset "setup_files" begin
        @testset "works for full method" begin
            apply(get_patch) do
                try
                    path = AdventOfCode.setup_files(2019, 1, force = true)
                    dir_path = @__DIR__
                    joinpath(dir_path, "src", "2019", "day_1.jl")
                    @test isdir("data/2019")
                    @test isfile("data/2019/day_1.txt")
                    @test readlines("data/2019/day_1.txt") == ["TESTDATA"]
                    @test isdir("src/2019")
                    @test isfile("src/2019/day_1.jl")
                    @test String(read("src/2019/day_1.jl")) == AdventOfCode._template(2019, 1)
                    @test_logs (:warn, r"data/2019/day_1.txt already exists.") (:warn, r"src/2019/day_1.jl already exists.") AdventOfCode.setup_files(2019, 1, force = false)
                finally
                    rm("data", recursive = true, force = true)
                    rm("src", recursive = true, force = true)
                end
            end
        end
    end
end


text:
module Intelligence

using CSV
using DataFrames

export foxtrot2scores, golf2scores

function str2int(letter::AbstractString)::Int
    letter = string(letter)
    @assert(occursin(r"[A-Z]+", letter), "r\"[A-Z]+\" does not occur in $letter")
    ch = letter[1]
    i = Int(ch) - 64
end

n_equal(a::Array, b::Array)::Number = count(map(t -> t[1] == t[2], zip(a, b)))
n_equal(a::Array, b::AbstractString)::Number = n_equal(a, map(i -> b, 1:length(a)))
n_weetniet(preds::Array{Int64,1})::Int64 = length(filter(v -> v == "Weet ik niet", preds))
get_answers(i::Number, df::DataFrame)::Vector = Vector(df[i, :][3:end])
get_numbers(i::Number, df::DataFrame)::Array{Int64,1} = map(str2int, get_answers(i, df))

function n_correct(df::DataFrame, solutions::Array)::Array{Int64,1}
    map(i -> n_equal(get_numbers(i, df), solutions), 1:nrow(df))
end

function n_incorrect(df::DataFrame, solutions::Array)::Array{Int64,1}
    map(n -> length(solutions) - n, n_correct(df, solutions))
end

function n_weetniet(df::DataFrame)::Array{Int64,1}
    map(i -> n_equal(get_answers(i, df), "Weet ik niet"), 1:nrow(df))
end

add_column(df::DataFrame, name::Symbol, data::Array)::DataFrame = hcat(df, DataFrame(name => data))

function foxtrot2scores(foxtrot::DataFrame)::DataFrame
    # Source: https://github.com/compsy/hoe_gek_is/blob/master/db/questionnaires/intelligentie.rb
    solutions = [4, 3, 5, 2, 2, 4, 5, 3, 1, 4, 5]

    scores = select(foxtrot, [:id, :completed_at])
    scores = add_column(scores, :correct, n_correct(foxtrot, solutions))
    scores = add_column(scores, :incorrect, n_incorrect(foxtrot, solutions))
    scores = add_column(scores, :dontknow, n_weetniet(foxtrot))

    return scores
end

function golf2scores(golf::DataFrame)::DataFrame
    solutions = [6, 5, 3, 2, 2, 6, 1, 7, 3, 7, 1, 5, 7, 1, 6, 5, 2, 3, 5, 6, 1, 2, 7, 3]

    scores = select(golf, [:id, :completed_at])
    scores = add_column(scores, :correct, n_correct(golf, solutions))
    scores = add_column(scores, :incorrect, n_incorrect(golf, solutions))
    scores = add_column(scores, :dontknow, n_weetniet(golf))

    return scores
end

end # module


text:
# Infastructure to support generating overloads from rules.
function __init__()
    # Need to refresh rules when a package is loaded
    push!(Base.package_callbacks, pkgid -> refresh_rules())
end

# Holds all the hook functions that are invokes when a new rule is defined
const RRULE_DEFINITION_HOOKS = Function[]
const FRULE_DEFINITION_HOOKS = Function[]
_hook_list(::typeof(rrule)) = RRULE_DEFINITION_HOOKS
_hook_list(::typeof(frule)) = FRULE_DEFINITION_HOOKS

"""
    on_new_rule(hook, frule | rrule)

Register a `hook` function to run when new rules are defined.
The hook receives a signature type-type as input, and generally will use `eval` to define
an overload of an AD system's overloaded type
For example, using the signature type `Tuple{typeof(+), Real, Real}` to make
`+(::DualNumber, ::DualNumber)` call the `frule` for `+`.
A signature type tuple always has the form:
`Tuple{typeof(operation), typeof{pos_arg1}, typeof{pos_arg2}...}`, where `pos_arg1` is the
first positional argument.

The hooks are automatically run on new rules whenever a package is loaded.
They can be manually triggered by [`refresh_rules`](@ref).
When a hook is first registered with `on_new_rule` it is run on all existing rules.
"""
function on_new_rule(hook_fun, rule_kind)
    # apply the hook to the existing rules
    ret = map(_rule_list(rule_kind)) do method
        sig = _primal_sig(rule_kind, method)
        _safe_hook_fun(hook_fun, sig)
    end

    # register hook for new rules -- so all new rules get this function applied
    push!(_hook_list(rule_kind), hook_fun)
    return ret
end

"""
    clear_new_rule_hooks!(frule|rrule)

Clears all hooks that were registered with corresponding [`on_new_rule`](@ref).
This is useful for while working interactively to define your rule generating hooks.
If you previously wrong an incorrect hook, you can use this to get rid of the old one.

!!! warning
    This absolutely should not be used in a package, as it will break any other AD system
    using the rule hooks that might happen to be loaded.
"""
clear_new_rule_hooks!(rule_kind) = empty!(_hook_list(rule_kind))

"""
    _rule_list(frule | rrule)

Returns a list of all the methods of the currently defined rules of the given kind.
Excluding the fallback rule that returns `nothing` for every input.
"""
function _rule_list end
# The fallback rules are the only rules defined in ChainRulesCore & that is how we skip them
_rule_list(rule_kind) = (m for m in methods(rule_kind) if m.module != @__MODULE__)

const LAST_REFRESH_RRULE = Ref(0)
const LAST_REFRESH_FRULE = Ref(0)
last_refresh(::typeof(frule)) = LAST_REFRESH_FRULE
last_refresh(::typeof(rrule)) = LAST_REFRESH_RRULE

"""
    refresh_rules()
    refresh_rules(frule | rrule)

This triggers all [`on_new_rule`](@ref) hooks to run on any newly defined rules.
It is *automatically* run when ever a package is loaded.
It can also be manually called to run it directly, for example if a rule was defined
in the REPL or within the same file as the AD function.
"""
function refresh_rules()
    refresh_rules(frule);
    refresh_rules(rrule)
end

function refresh_rules(rule_kind)
    isempty(_rule_list(rule_kind)) && return  # if no hooks, exit early, nothing to run
    already_done_world_age = last_refresh(rule_kind)[]
    for method in _rule_list(rule_kind)
        _defined_world(method) < already_done_world_age && continue
        sig = _primal_sig(rule_kind, method)
        _trigger_new_rule_hooks(rule_kind, sig)
    end

    last_refresh(rule_kind)[] = _current_world()
    return nothing
end

@static if VERSION >= v"1.2"
    _current_world() = Base.get_world_counter()
    _defined_world(method) = method.primary_world
else
    _current_world() = ccall(:jl_get_world_counter, UInt, ())
    _defined_world(method) = method.min_world
end

"""
    _primal_sig(frule|rule, rule_method | rule_sig)

Returns the signature as a `Tuple{function_type, arg1_type, arg2_type,...}`.
"""
_primal_sig(rule_kind, method::Method) = _primal_sig(rule_kind, method.sig)
function _primal_sig(::typeof(frule), rule_sig::DataType)
    @assert rule_sig.parameters[1] == typeof(frule)
    # need to skip frule and the deriviative info, so starting from the 3rd
    return Tuple{rule_sig.parameters[3:end]...}
end
function _primal_sig(::typeof(rrule), rule_sig::DataType)
    @assert rule_sig.parameters[1] == typeof(rrule)
    # need to skip rrule so starting from the 2rd
    return Tuple{rule_sig.parameters[2:end]...}
end
function _primal_sig(rule_kind, rule_sig::UnionAll)
    # This looks a lot like Base.unwrap_unionall and Base.rewrap_unionall, but using those
    # seems not to work
    p_sig = _primal_sig(rule_kind, rule_sig.body)
    return UnionAll(rule_sig.var, p_sig)
end

function _trigger_new_rule_hooks(rule_kind, sig)
    for hook_fun in _hook_list(rule_kind)
        _safe_hook_fun(hook_fun, sig)
    end
end

function _safe_hook_fun(hook_fun, sig)
    try
        hook_fun(sig)
    catch err
        @error "Error triggering hook" hook_fun sig exception=(err, catch_backtrace())
    end
end


text:
"""Gas diffusivity in air under STP condition [m^2 s^-1]."""
const diffus_air_stp = Dict{String, Float64}(
    "h2o" => 2.178e-5,
    "co2" => 1.381e-5,
    "ch4" => 1.952e-5,
    "co" => 1.807e-5,
    "so2" => 1.089e-5,
    "o3" => 1.444e-5,
    "nh3" => 1.978e-5,
    "n2o" => 1.436e-5,
    "no" => 1.802e-5,
    "no2" => 1.361e-5,
    "n2" => 1.788e-5,
    "o2" => 1.820e-5,
    "cos" => 1.381e-5 / 1.21)

"""
Pre-exponential factor and activation energy for the calculation of gas
diffusivity in water [m^2 s^-1] under STP condition.
"""
const diffus_water_stp = Dict{String, Pair{Float64, Float64}}(
    "he" => Pair(818e-9, 11.70e3),
    "ne" => Pair(1608e-9, 14.84e3),
    "kr" => Pair(6393e-9, 20.20e3),
    "xe" => Pair(9007e-9, 21.61e3),
    "rn" => Pair(15877e-9, 23.26e3),
    "h2" => Pair(3338e-9, 16.06e3),
    "ch4" => Pair(3047e-9, 18.36e3),
    "co2" => Pair(5019e-9, 19.51e3),
    "cos" => Pair(4.735872481253359e-6, 19336.20405260121),
    "co" => Pair(0.407e-4, 24518.24),
    "no" => Pair(39.8e-4, 34978.24))

"""Clapp-Hornberger shape parameters for tortuosity effect."""
const soil_shape_params = Dict{String, Float64}(
    "sand" => 4.05,
    "loamy sand" => 4.38,
    "sandy loam" => 4.9,
    "silt loam" => 5.30,
    "loam" => 5.39,
    "sandy clay loam" => 7.12,
    "silty clay loam" => 7.75,
    "clay loam" => 8.52,
    "sandy clay" => 10.4,
    "silty clay" => 10.4,
    "clay" => 11.4)

"""
    diffus_air(species::String, temp[, pressure=Constants.atm])

Calculate the diffusivity of a gas in air.

Supported gases include H2O, CO2, CH4, CO, SO2, O3, NH3, N2O, NO, NO2, N2, O2
and COS.

# Parameters

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `temp`: Temperature in Kelvin.
* `pressure` (*optional*): Ambient pressure in Pascal. Default is the standard
  atmospheric pressure.

# Returns

Diffusivity of the gas in air [m^2 s^-1].

# See also

[`diffus_soil_air`](@ref): Calculate the diffusivity of a gas in soil air.

# References

[M98] Massman, W. J. (1998). A review of the molecular diffusivities of
    H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2 in air, O2 and N2
    near STP. *Atmos. Environ.*, 32(6), 1111--1127.
    <https://doi.org/10.1016/S1352-2310(97)00391-9>.

[S10] Seibt, U. et al. (2010). A kinetic analysis of leaf uptake of COS
    and its relation to transpiration, photosynthesis and carbon isotope
    fractionation. *Biogeosci.*, 7, 333--341.
    <https://doi.org/10.5194/bg-7-333-2010>.
"""
function diffus_air(species::String, temp, pressure=Constants.atm)
    if !haskey(diffus_air_stp, species)
        throw(ArgumentError("species '" * species * "' not supported"))
    end

    return @. (diffus_air_stp[species] * (Constants.atm / pressure) *
               (temp / Constants.T_0) ^ 1.81)
end

"""
    diffus_water(species::String, temp)

Calculate the diffusivity of a gas in water.

Supported gases include He, Ne, Kr, Xe, Rn, H2, CH4, CO2, CO, NO, and COS.

# Parameters

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `temp`: Temperature in Kelvin.

# Returns

Diffusivity of the gas in water [m^2 s^-1].

# See also

[`diffus_soil_water`](@ref): Calculate the diffusivity of a gas in soil water.

# References

[J87] Jähne, B. et al. (1987). Measurement of the diffusion
    coefficients of sparingly soluble gases in water. *J. Geophys. Res.*,
    92(C10), 10767--10776. <https://doi.org/10.1029/JC092iC10p10767>.

[WH68] Wise, D. L. and Houghton, G. (1968). Diffusion coefficients of
    neon, krypton, xenon, carbon monoxide and nitric oxide in water
    at 10-6O C. *Chem. Eng. Sci.*, 23, 1211--1216.
    <https://doi.org/10.1016/0009-2509(68)89029-3>.

[UFUA96] Ulshöfer, V. S., Flöck, O. R., Uher, G., and Andreae, M. O.
    (1996). Photochemical production and air-sea exchange of sulfide in
    the eastern Mediterranean Sea. *Mar. Chem.*, 53, 25--39.
    <https://doi.org/10.1016/0304-4203(96)00010-2>.
"""
function diffus_water(species::String, temp)
    if !haskey(diffus_water_stp, species)
        throw(ArgumentError("species \"" * species * "\" not supported"))
    end

    # diffus_water_stp[species] is a pair of (pre-exp factor, E_act)
    return @. (diffus_water_stp[species].first *
               exp(- diffus_water_stp[species].second / (Constants.R * temp)))
end

"""
    diffus_soil_air(species::String, texture::String, temp, theta_sat,
                    theta_w[, pressure=Constants.atm])

Calculate the diffusivity of a gas in soil air.

Supported gases include H2O, CO2, CH4, CO, SO2, O3, NH3, N2O, NO, NO2, N2, O2
and COS.

# Parameters

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `texture`: Soil texture name in lower case.
* `temp`: Temperature in Kelvin.
* `theta_sat`: Total porosity of the soil [m^3 m^-3].
* `theta_w`: Water-filled porosity of the soil [m^3 m^-3].
* `pressure` (*optional*): Ambient pressure in Pascal. Default is the standard
  atmospheric pressure.

# Returns

Diffusivity of the gas in soil air [m^2 s^-1].

# See also

* [`diffus_air`](@ref): Calculate the diffusivity of a gas in air.
* [`diffus_soil_water`](@ref): Calculate the diffusivity of a gas in soil
  water.
* [`diffus_soil`](@ref): Calculate the total diffusivity of a gas in soil
  gaseous and aqueous phases.

# References

[CH78] Clapp, R. B. and Hornberger, G. M. (1978). Empirical equations
    for some soil hydraulic properties. *Water Resources Res.*, 14(4),
    601--604. <https://doi.org/10.1029/WR014i004p00601>.

[M03] Moldrup, P., Olesen, T., Komatsu, T, Yoshikawa, S, Schjønning, P, and
    Rolston, D. E. (2003). Modeling diffusion and reaction in soils: X. A
    unifying model for solute and gas diffusivity in unsaturated soil. *Soil
    Sci.*, 168(5), 321--337.
    <https://doi.org/10.1097/01.ss.0000070907.55992.3c>.
"""
function diffus_soil_air(species::String, texture::String, temp, theta_sat,
                         theta_w, pressure=Constants.atm)
    if !haskey(soil_shape_params, texture)
        throw(ArgumentError("soil texture \"" * texture * "\" not supported"))
    end

    # tortuosity in soil air
    theta_a = theta_sat - theta_w
    tau_a = @. (theta_a * (theta_a / theta_sat) ^
                (3. / soil_shape_params[texture]))
    return @. (diffus_air(species, temp, pressure) * theta_a * tau_a)
end

"""
    diffus_soil_water(species::String, texture::String, temp, theta_sat,
                      theta_w)

Calculate the diffusivity of a gas in soil water.

Supported gases include He, Ne, Kr, Xe, Rn, H2, CH4, CO2, CO, NO, and COS.

# Parameters

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `texture`: Soil texture name in lower case.
* `temp`: Temperature in Kelvin.
* `theta_sat`: Total porosity of the soil [m^3 m^-3].
* `theta_w`: Water-filled porosity of the soil [m^3 m^-3].

# Returns

Diffusivity of the gas in soil water [m^2 s^-1].

# See also

* [`diffus_water`](@ref): Calculate the diffusivity of a gas in water.
* [`diffus_soil_air`](@ref): Calculate the diffusivity of a gas in soil air.
* [`diffus_soil`](@ref): Calculate the total diffusivity of a gas in soil
  gaseous and aqueous phases.

# References

[CH78] Clapp, R. B. and Hornberger, G. M. (1978). Empirical equations
    for some soil hydraulic properties. *Water Resources Res.*, 14(4),
    601--604. <https://doi.org/10.1029/WR014i004p00601>.

[M03] Moldrup, P., Olesen, T., Komatsu, T, Yoshikawa, S, Schjønning, P, and
    Rolston, D. E. (2003). Modeling diffusion and reaction in soils: X. A
    unifying model for solute and gas diffusivity in unsaturated soil. *Soil
    Sci.*, 168(5), 321--337.
    <https://doi.org/10.1097/01.ss.0000070907.55992.3c>.
"""
function diffus_soil_water(species::String, texture::String, temp, theta_sat,
                           theta_w)
    if !haskey(soil_shape_params, texture)
        throw(ArgumentError("soil texture \"" * texture * "\" not supported"))
    end

    # tortuosity in soil water
    tau_w = @. (theta_w * (theta_w / theta_sat) ^
                (soil_shape_params[texture] / 3. - 1.))
    return @. (diffus_water(species, temp) * theta_w * tau_w)
end

"""
    diffus_soil(species::String, texture::String, temp, theta_sat, theta_w[,
                pressure=Constants.atm])

Calculate the total diffusivity of a gas in soil gaseous and aqueous phases.

Supported gases include CO2, CO, NO, CH4, and COS.

# Parameters

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `texture`: Soil texture name in lower case.
* `temp`: Temperature in Kelvin.
* `theta_sat`: Total porosity of the soil [m^3 m^-3].
* `theta_w`: Water-filled porosity of the soil [m^3 m^-3].
* `pressure` (*optional*): Ambient pressure in Pascal. Default is the standard
  atmospheric pressure.

# Returns

Dual-phase diffusivity of the gas in soil [m^2 s^-1].

# See also

* [`diffus_soil_air`](@ref): Calculate the diffusivity of a gas in soil air.
* [`diffus_soil_water`](@ref): Calculate the diffusivity of a gas in soil
  water.

# References

[CH78] Clapp, R. B. and Hornberger, G. M. (1978). Empirical equations
    for some soil hydraulic properties. *Water Resources Res.*, 14(4),
    601--604. <https://doi.org/10.1029/WR014i004p00601>.

[M03] Moldrup, P., Olesen, T., Komatsu, T, Yoshikawa, S, Schjønning, P, and
    Rolston, D. E. (2003). Modeling diffusion and reaction in soils: X. A
    unifying model for solute and gas diffusivity in unsaturated soil. *Soil
    Sci.*, 168(5), 321--337.
    <https://doi.org/10.1097/01.ss.0000070907.55992.3c>.
"""
function diffus_soil(species::String, texture::String, temp, theta_sat,
                     theta_w, pressure=Constants.atm)
    if !haskey(soil_shape_params, texture)
        throw(ArgumentError("soil texture \"" * texture * "\" not supported"))
    end

    if !(haskey(diffus_air_stp, species) && haskey(diffus_water_stp, species))
        throw(ArgumentError("species \"" * species * "\" not supported"))
    end

    return diffus_soil_water(species, texture, temp, theta_sat, theta_w) +
        diffus_soil_air(species, texture, temp, theta_sat, theta_w, pressure)
end


text:
# Source for Cartesian package. License is MIT https://github.com/ja72/Cartesian.jl/blob/master/LICENSE

"""

Definitions for cartesian vectors and matrices used in 3D mechanics.
The repository is located at https://github.com/ja72/Cartesian.jl

    $(EXPORTS)

    $(SIGNATURES)

    $(METHODLIST)

"""
module Cartesian
using LinearAlgebra
using StaticArrays
using Rotations
using DocStringExtensions

using Base: @propagate_inbounds

import Base.convert, Base.*
import LinearAlgebra.cross, LinearAlgebra.dot, LinearAlgebra.normalize, LinearAlgebra.norm

export
    Vector3,
    Matrix3,
    dot,
    cross,
    cross2,
    slerp,
    outer,
    inv,
    solve

include("Vector3.jl")
include("Matrix3.jl")

end


text:
" Implements a CoupledConvection operator "
type CoupledConvection <: Kernel
    u::Variable
    other_var::Variable
end

import MOOSE.computeQpResidual

@inline function computeQpResidual(kernel::CoupledConvection, qp::Int64, i::Int64)
    u = kernel.u

    other_var = kernel.other_var

    return (other_var.grad[qp] ⋅ u.grad[qp]) * u.phi[qp][i]
end

import MOOSE.computeQpJacobian

@inline function computeQpJacobian(kernel::CoupledConvection, v::Variable, qp::Int64, i::Int64, j::Int64)::Float64
    u = kernel.u

    other_var = kernel.other_var

    if u.id == v.id
        return (other_var.grad[qp] ⋅ u.grad_phi[qp][j]) * u.phi[qp][i]
    elseif other_var.id == v.id
        return (other_var.grad_phi[qp][j] ⋅ u.grad[qp]) * u.phi[qp][i]
    end

    return 0
end

import MOOSE.coupledVars

@inline function coupledVars(kernel::CoupledConvection)::Array{Variable}
    return [kernel.u, kernel.other_var]
end


text:
using Multilane
using StatsBase
using MCTS
using JLD

behaviors=Dict{String,Any}(
    "correlated" => standard_uniform(1.0, correlated=true)
)

dpws = DPWSolver(depth=20,
                 n_iterations=500,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

solvers = Dict{String, Any}(
    "dpw"=>dpws,
    "assume_normal"=>SingleBehaviorSolver(dpws, Multilane.NORMAL)
)

curve = TestSet(lambda=[1.0, 2.0, 4.0, 8.0, 16.0, 32.0], N=500, p_appear=1.0)

tests = []
for b in [2.5, 3.0, 3.5, 4.0]
    push!(tests, TestSet(curve, solver_key="dpw",
                         behaviors="correlated",
                         key=@sprintf("upper_bound_%02d", b*10),
                         brake_threshold=b))
    push!(tests, TestSet(curve, solver_key="assume_normal",
                         behaviors="correlated",
                         key=@sprintf("assume_normal_%02d", b*10),
                         brake_threshold=b))
end

objects = gen_initials(tests, behaviors=behaviors, generate_physical=true)

@show objects["param_table"]
objects["solvers"] = solvers

objects["note"] = "Tests over a range of braking values with correlated behaviors."

sbatch_spawn(tests, objects,
             batch_size=50,
             time_per_batch="1:00:00",
             submit_command="sbatch",
             template_name="sherlock.sh")

# sbatch_spawn(tests, objects,
#              batch_size=12,
#              time_per_batch="1:00:00",
#              submit_command="bash",
#              template_name="theresa.sh")
#
# results = evaluate(tests, objects, parallel=true)
#
# filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
# results["histories"] = nothing
# save(filename, results)
# println("results saved to $filename")


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to ECEF to ECI transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Get the current EOP Data.
#
# TODO: The EOP data obtained from IERS website does not match the values in the
# examples in [1]. However, this should be enough, because 1) the individual
# functions at the low level are tested using the same values of [1], and 2) the
# difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("./eop_IAU1980.txt",  :IAU1980)
eop_iau2000a = read_iers_eop("./eop_IAU2000A.txt", :IAU2000A)

# File: ./src/transformations/ecef_to_eci.jl
# ==========================================

# Functions: rECEFtoECI
# ---------------------

################################################################################
#                                  IAU-76/FK5
################################################################################

## ITRF to GCRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_gcrf = 5102.50895790  i + 6123.01140070   j + 6378.13692820   k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => GCRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_GCRF_ITRF = rECEFtoECI(ITRF(), GCRF(), JD_UTC, eop_iau1980)
    r_gcrf = D_GCRF_ITRF*r_itrf

    @test r_gcrf[1] ≈ +5102.50895790 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=3e-4

    q_GCRF_ITRF = rECEFtoECI(Quaternion, ITRF(), GCRF(), JD_UTC, eop_iau1980)
    r_gcrf = vect(conj(q_GCRF_ITRF)*r_itrf*q_GCRF_ITRF)

    @test r_gcrf[1] ≈ +5102.50895790 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=3e-4
end

## ITRF to J2000
## =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => J2000" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_J2000_ITRF = rECEFtoECI(ITRF(), J2000(), JD_UTC, eop_iau1980)
    r_j2000 = D_J2000_ITRF*r_itrf

    @test r_j2000[1] ≈ +5102.50960000 atol=3e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=3e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=3e-4

    q_J2000_ITRF = rECEFtoECI(Quaternion, ITRF(), J2000(), JD_UTC, eop_iau1980)
    r_j2000 = vect(conj(q_J2000_ITRF)*r_itrf*q_J2000_ITRF)

    @test r_j2000[1] ≈ +5102.50960000 atol=3e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=3e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=3e-4
end

## ITRF to TOD
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => TOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_TOD_ITRF = rECEFtoECI(ITRF(), TOD(), JD_UTC, eop_iau1980)
    r_tod = D_TOD_ITRF*r_itrf

    @test r_tod[1] ≈ +5094.51620300 atol=3e-4
    @test r_tod[2] ≈ +6127.36527840 atol=3e-4
    @test r_tod[3] ≈ +6380.34453270 atol=3e-4

    q_TOD_ITRF = rECEFtoECI(Quaternion, ITRF(), TOD(), JD_UTC, eop_iau1980)
    r_tod = vect(conj(q_TOD_ITRF)*r_itrf*q_TOD_ITRF)

    @test r_tod[1] ≈ +5094.51620300 atol=3e-4
    @test r_tod[2] ≈ +6127.36527840 atol=3e-4
    @test r_tod[3] ≈ +6380.34453270 atol=3e-4
end

## ITRF to MOD
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => MOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_MOD_ITRF = rECEFtoECI(ITRF(), MOD(), JD_UTC, eop_iau1980)
    r_mod = D_MOD_ITRF*r_itrf

    @test r_mod[1] ≈ +5094.02837450 atol=3e-4
    @test r_mod[2] ≈ +6127.87081640 atol=3e-4
    @test r_mod[3] ≈ +6380.24851640 atol=3e-4

    q_MOD_ITRF = rECEFtoECI(Quaternion, ITRF(), MOD(), JD_UTC, eop_iau1980)
    r_mod = vect(conj(q_MOD_ITRF)*r_itrf*q_MOD_ITRF)

    @test r_mod[1] ≈ +5094.02837450 atol=3e-4
    @test r_mod[2] ≈ +6127.87081640 atol=3e-4
    @test r_mod[3] ≈ +6380.24851640 atol=3e-4
end

## ITRF to TEME
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830   i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => TEME" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_TEME_ITRF = rECEFtoECI(ITRF(), TEME(), JD_UTC, eop_iau1980)
    r_teme = D_TEME_ITRF*r_itrf

    @test r_teme[1] ≈ +5094.18016210 atol=3e-4
    @test r_teme[2] ≈ +6127.64465950 atol=3e-4
    @test r_teme[3] ≈ +6380.34453270 atol=3e-4

    q_TEME_ITRF = rECEFtoECI(Quaternion, ITRF(), TEME(), JD_UTC, eop_iau1980)
    r_teme = vect(conj(q_TEME_ITRF)*r_itrf*q_TEME_ITRF)

    @test r_teme[1] ≈ +5094.18016210 atol=3e-4
    @test r_teme[2] ≈ +6127.64465950 atol=3e-4
    @test r_teme[3] ≈ +6380.34453270 atol=3e-4
end

## PEF to GCRF
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_gcrf = 5102.50895790  i + 6123.01140070   j + 6378.13692820   k [km]
#
################################################################################

@testset "Function rECEFtoECI PEF => GCRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]

    D_GCRF_PEF = rECEFtoECI(PEF(), GCRF(), JD_UTC, eop_iau1980)
    r_gcrf = D_GCRF_PEF*r_pef

    @test r_gcrf[1] ≈ +5102.50895790 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=3e-4

    q_GCRF_PEF = rECEFtoECI(Quaternion, PEF(), GCRF(), JD_UTC, eop_iau1980)
    r_gcrf = vect(conj(q_GCRF_PEF)*r_pef*q_GCRF_PEF)

    @test r_gcrf[1] ≈ +5102.50895790 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=3e-4
end

## PEF to J2000
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
################################################################################

@testset "Function rECEFtoECI PEF => J2000" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]

    D_J2000_PEF = rECEFtoECI(PEF(), J2000(), JD_UTC, eop_iau1980)
    r_j2000 = D_J2000_PEF*r_pef

    @test r_j2000[1] ≈ +5102.50960000 atol=3e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=3e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=3e-4

    q_J2000_PEF = rECEFtoECI(Quaternion, PEF(), J2000(), JD_UTC, eop_iau1980)
    r_j2000 = vect(conj(q_J2000_PEF)*r_pef*q_J2000_PEF)

    @test r_j2000[1] ≈ +5102.50960000 atol=3e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=3e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=3e-4
end

## PEF to TOD
## ==========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
# If not EOP correction is used, then we get:
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#
# Notice, however, that this result is obtained by using the correct UT1.
#
################################################################################

@testset "Function rECEFtoECI PEF => TOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]

    D_TOD_PEF = rECEFtoECI(PEF(), TOD(), JD_UTC, eop_iau1980)
    r_tod = D_TOD_PEF*r_pef

    @test r_tod[1] ≈ +5094.51620300 atol=3e-4
    @test r_tod[2] ≈ +6127.36527840 atol=3e-4
    @test r_tod[3] ≈ +6380.34453270 atol=3e-4

    q_TOD_PEF = rECEFtoECI(Quaternion, PEF(), TOD(), JD_UTC, eop_iau1980)
    r_tod = vect(conj(q_TOD_PEF)*r_pef*q_TOD_PEF)

    @test r_tod[1] ≈ +5094.51620300 atol=3e-4
    @test r_tod[2] ≈ +6127.36527840 atol=3e-4
    @test r_tod[3] ≈ +6380.34453270 atol=3e-4

    # No EOP corrections
    # ==================

    # Notice that if we use JD_UT1, then the correction to TT will be wrong.
    # However, this will lead to a much smaller error than assuming that UTC =
    # UT1.

    JD_UT1 = DatetoJD(2004,4,6,7,51,28.386009) - 0.4399619/86400

    D_TOD_PEF = rECEFtoECI(PEF(), TOD(), JD_UT1)
    r_tod = D_TOD_PEF*r_pef

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    q_TOD_PEF = rECEFtoECI(Quaternion, PEF(), TOD(), JD_UT1)
    r_tod = vect(conj(q_TOD_PEF)*r_pef*q_TOD_PEF)

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7
end

## PEF to MOD
## ==========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
# If not EOP correction is used, then we get:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#
# Notice, however, that this result is obtained by using the correct UT1.
#
################################################################################

@testset "Function rECEFtoECI PEF => MOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]

    D_MOD_PEF = rECEFtoECI(PEF(), MOD(), JD_UTC, eop_iau1980)
    r_mod = D_MOD_PEF*r_pef

    @test r_mod[1] ≈ +5094.02837450 atol=3e-4
    @test r_mod[2] ≈ +6127.87081640 atol=3e-4
    @test r_mod[3] ≈ +6380.24851640 atol=3e-4

    q_MOD_PEF = rECEFtoECI(Quaternion, PEF(), MOD(), JD_UTC, eop_iau1980)
    r_mod = vect(conj(q_MOD_PEF)*r_pef*q_MOD_PEF)

    @test r_mod[1] ≈ +5094.02837450 atol=3e-4
    @test r_mod[2] ≈ +6127.87081640 atol=3e-4
    @test r_mod[3] ≈ +6380.24851640 atol=3e-4

    # No EOP corrections
    # ==================

    # Notice that if we use JD_UT1, then the correction to TT will be wrong.
    # However, this will lead to a much smaller error than assuming that UTC =
    # UT1.

    JD_UT1 = DatetoJD(2004,4,6,7,51,28.386009) - 0.4399619/86400

    D_MOD_PEF = rECEFtoECI(PEF(), MOD(), JD_UT1)
    r_mod = D_MOD_PEF*r_pef

    @test r_mod[1] ≈ +5094.02901670 atol=8e-6
    @test r_mod[2] ≈ +6127.87093630 atol=8e-6
    @test r_mod[3] ≈ +6380.24788850 atol=8e-6

    q_MOD_PEF = rECEFtoECI(Quaternion, PEF(), MOD(), JD_UT1)
    r_mod = vect(conj(q_MOD_PEF)*r_pef*q_MOD_PEF)

    @test r_mod[1] ≈ +5094.02901670 atol=8e-6
    @test r_mod[2] ≈ +6127.87093630 atol=8e-6
    @test r_mod[3] ≈ +6380.24788850 atol=8e-6
end

## PEF to TEME
## ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECEFtoECI PEF => TEME" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]

    D_TEME_PEF = rECEFtoECI(PEF(), TEME(), JD_UTC, eop_iau1980)
    r_teme = D_TEME_PEF*r_pef

    @test r_teme[1] ≈ +5094.18016210 atol=3e-4
    @test r_teme[2] ≈ +6127.64465950 atol=3e-4
    @test r_teme[3] ≈ +6380.34453270 atol=3e-4

    q_TEME_PEF = rECEFtoECI(Quaternion, PEF(), TEME(), JD_UTC, eop_iau1980)
    r_teme = vect(conj(q_TEME_PEF)*r_pef*q_TEME_PEF)

    @test r_teme[1] ≈ +5094.18016210 atol=3e-4
    @test r_teme[2] ≈ +6127.64465950 atol=3e-4
    @test r_teme[3] ≈ +6380.34453270 atol=3e-4
end

################################################################################
#                                IAU-2006/2010
################################################################################

## ITRF => CIRS
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => CIRS" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_CIRS_ITRF = rECEFtoECI(ITRF(), CIRS(), JD_UTC, eop_iau2000a)
    r_cirs = D_CIRS_ITRF*r_itrf

    @test r_cirs[1] ≈ +5100.01840470 atol=3e-4
    @test r_cirs[2] ≈ +6122.78636480 atol=3e-4
    @test r_cirs[3] ≈ +6380.34453270 atol=3e-4

    q_CIRS_ITRF = rECEFtoECI(Quaternion, ITRF(), CIRS(), JD_UTC, eop_iau2000a)
    r_cirs = vect(q_CIRS_ITRF\r_itrf*q_CIRS_ITRF)

    @test r_cirs[1] ≈ +5100.01840470 atol=3e-4
    @test r_cirs[2] ≈ +6122.78636480 atol=3e-4
    @test r_cirs[3] ≈ +6380.34453270 atol=3e-4
end

## ITRF => GCRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#
# one gets:
#
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#
################################################################################

@testset "Function rECEFtoECI ITRF => GCRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]

    D_GCRF_ITRF = rECEFtoECI(ITRF(), GCRF(), JD_UTC, eop_iau2000a)
    r_gcrf = D_GCRF_ITRF*r_itrf

    @test r_gcrf[1] ≈ +5102.50895290 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol=3e-4

    q_GCRF_ITRF = rECEFtoECI(Quaternion, ITRF(), GCRF(), JD_UTC, eop_iau2000a)
    r_gcrf = vect(q_GCRF_ITRF\r_itrf*q_GCRF_ITRF)

    @test r_gcrf[1] ≈ +5102.50895290 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol=3e-4
end

## TIRS => CIRS
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_tirs = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#
################################################################################

@testset "Function rECEFtoECI TIRS => CIRS" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_tirs = [-1033.47503120; 7901.30558560; 6380.34453270]

    D_CIRS_TIRS = rECEFtoECI(TIRS(), CIRS(), JD_UTC, eop_iau2000a)
    r_cirs = D_CIRS_TIRS*r_tirs

    @test r_cirs[1] ≈ +5100.01840470 atol=3e-4
    @test r_cirs[2] ≈ +6122.78636480 atol=3e-4
    @test r_cirs[3] ≈ +6380.34453270 atol=3e-4

    q_CIRS_TIRS = rECEFtoECI(Quaternion, TIRS(), CIRS(), JD_UTC, eop_iau2000a)
    r_cirs = vect(q_CIRS_TIRS\r_tirs*q_CIRS_TIRS)

    @test r_cirs[1] ≈ +5100.01840470 atol=3e-4
    @test r_cirs[2] ≈ +6122.78636480 atol=3e-4
    @test r_cirs[3] ≈ +6380.34453270 atol=3e-4
end

## TIRS => GCRF
## ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_tirs = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#
################################################################################

@testset "Function rECEFtoECI TIRS => GCRF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)
    r_tirs = [-1033.47503120; 7901.30558560; 6380.34453270]

    D_GCRF_TIRS = rECEFtoECI(TIRS(), GCRF(), JD_UTC, eop_iau2000a)
    r_gcrf = D_GCRF_TIRS*r_tirs

    @test r_gcrf[1] ≈ +5102.50895290 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol=3e-4

    q_GCRF_TIRS = rECEFtoECI(Quaternion, TIRS(), GCRF(), JD_UTC, eop_iau2000a)
    r_gcrf = vect(q_GCRF_TIRS\r_tirs*q_GCRF_TIRS)

    @test r_gcrf[1] ≈ +5102.50895290 atol=3e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol=3e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol=3e-4
end


text:
#### randomized SVD (from Jiahao Chen, based on http://arxiv.org/pdf/0909.4061.pdf)

import Base.LinAlg.SVD

#The simplest possible randomized svd
#Inputs
#   A: input matrix
#   n: Number of singular value/vector pairs to find
#   p: Number of extra vectors to include in computation
function rsvd(A, n, p=0)
    Q = rrange(A, n, p=p)
    rsvd_direct(A, Q)
end

#Algorithm 4.4: randomized subspace iteration
#A must support size(A), multiply and transpose multiply
#p is the oversampling parameter
#q controls the accuracy of the subspace found; it is the "number of power iterations"
#A good heuristic is that when the original scheme produces a basis whose
#approximation error is within a factor C of the optimum, the power scheme produces
#an approximation error within C^(1/(2q+1)) of the optimum.
function rrange(A, l::Integer; p::Integer=5, q::Integer=3)
    p≥0 || error()
    m, n = size(A)
    l <= m || error("Cannot find $l linearly independent vectors of $m x $n matrix")
    Ω = randn(n, l+p)
    Q = q_from_qr(A*Ω)
    for t=1:q
        Q = q_from_qr(A'*Q)
        Q = q_from_qr(A*Q)
    end
    Q = p==0 ? Q : Q[:,1:l]
end

function q_from_qr(Y, l::Integer=-1)
    Q = full(qrfact!(Y)[:Q])
    Q = l<0 ? Q : Q[:,1:l]
end

#Algorithm 5.1: direct SVD
#More accurate
function rsvd_direct(A, Q)
    B=Q'A
    S=svdfact!(B)
    SVD(Q*S[:U], S[:S], S[:Vt])
end

function onepass_svd(A::AbstractArray, r::Int)
  m, n = size(A)
  k = 2r + 1
  l = 4r + 3
  Omega = randn(n,k)
  Psi   = randn(m,l)
  Y     = A*Omega
  W     = A'*Psi
  Q,_ = qr(view(Y,:,1:k))
  B = view(W,:,1:l) / (Q'*view(Psi,:,1:l)) # Q's.Psi is k x l, its pinv is l x k, so B is n x k
  mysvd,_ = svds(B, nsv=r) # U is n x r
  return SVD(Q*mysvd.Vt, mysvd.S, mysvd.U)
end


text:
module Allsvenskan

using DataFrames
using CSV
using Missings
using StringDistances
using Dates
using Printf

export
	addgame!,
	table,
	mergetables,
    allgamesbyseason,
    allgames,
    alltables,
    allteams,
    marathon,
    showmarathon,
    showhead2head,
    showgames

include("table.jl")
include("utils.jl")
include("show.jl")
include("city.jl")

end


text:
module SemanticReleaseDemo

# Write your package code here.

f() = "hello world"

end


text:
using RiskAdjustedLinearizations, Test, JLD2
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "wachter_disaster_risk", "wachter.jl"))

# Solve model
m_wachter = WachterDisasterRisk()
m = inplace_wachter_disaster_risk(m_wachter)
try
    solve!(m, m.z, m.y; verbose = :none)
catch e
    local sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "iterative_sss_output.jld2"), "r")
    update!(m, sssout["z"], sssout["y"], sssout["Psi"])
end

# Simulate with no shocks
horizon = 100
zero_shocks = zeros(3, horizon) # 100 fake draws
@testset "Simulate Wachter (2013) with no shocks" begin
    state1, jump1 = simulate(m, horizon)
    state2, jump2 = simulate(m, horizon, m.z)
    state3, jump3 = simulate(m, horizon, zero_shocks)
    state4, jump4 = simulate(m, horizon, zero_shocks, m.z)
    state5, jump5 = simulate(m, horizon, 1.01 * m.z) # perturb away from steady state
    state6, jump6 = simulate(m, horizon, zero_shocks, 1.01 * m.z)
    state7, jump7 = simulate(m, vec(zero_shocks[:, 1]), 1.01 * m.z)

    @test state1 ≈ state2
    @test state1 ≈ state3
    @test state1 ≈ state4
    @test jump1 ≈ jump2
    @test jump1 ≈ jump3
    @test jump1 ≈ jump4
    @test state5 ≈ state6
    @test jump5 ≈ jump6
    @test state1 ≈ repeat(m.z, 1, horizon) # check state1 remains at steady state
    @test jump1 ≈ repeat(m.y, 1, horizon)
    @test state7 ≈ state6[:, 1]
    @test jump7 ≈ jump6[:, 1]
end

# 100 draws from a standard multivariate normal, technically not the right distribution but it's fine
shocks = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "wachter_shocks.jld2"), "r")["shocks"]

@testset "Simulate Wachter (2013) with shocks" begin
    state1, jump1 = simulate(m, horizon)
    state3, jump3 = simulate(m, horizon, shocks)
    state4, jump4 = simulate(m, horizon, shocks, m.z)
    state5, jump5 = simulate(m, shocks[:, 1], m.z)

    @test state3 ≈ state4
    @test jump3 ≈ jump4
    @test !(state1 ≈ state4)
    @test !(jump1 ≈ jump4)
    @test state5 ≈ state4[:, 1]
    @test jump5 ≈ jump4[:, 1]
end


text:
struct KMarkovDecUpdater <: Updater
    up::KMarkovUpdater
    problem::UrbanPOMDP
end

const KMarkovDecBelief = Vector{Tuple{Symbol, Array{Float64, 2}}}

function initialize_dec_belief(up::KMarkovDecUpdater, pomdp::UrbanPOMDP, rng::AbstractRNG)
    s = initialstate(pomdp, rng)
    o = initialobs(pomdp, s, rng)
    inputs = decompose_input(pomdp, o)
    init_b = initialize_dec_belief(up, inputs)
    return init_b
end

function initialize_dec_belief(up::KMarkovDecUpdater, pomdp::UrbanPOMDP, s::UrbanState, rng::AbstractRNG)
    o = initialobs(pomdp, s, rng)
    inputs = decompose_input(pomdp, o)
    init_b = initialize_dec_belief(up, inputs)
    return init_b
end

function initialize_dec_belief(up::KMarkovDecUpdater, pomdp::UrbanPOMDP, o::UrbanObs)
    inputs = decompose_input(pomdp, o)
    init_b = initialize_dec_belief(up, inputs)
    return init_b
end

function initialize_dec_belief(up::KMarkovDecUpdater, dec_o::Vector{Tuple{Symbol, Vector{Float64}}})
    init_b = Tuple{Symbol, Array{Float64,2}}[]
    for (prob_key, obs) in dec_o
        b_stacked = hcat(initialize_belief(up.up, fill(obs, up.up.k))...)
        push!(init_b, (prob_key, b_stacked))
    end
    return init_b
end

function POMDPs.update(up::KMarkovDecUpdater, b_old::KMarkovDecBelief, a::UrbanAction, o::Array{Float64})
    b_new = Tuple{Symbol, Array{Float64, 2}}[]
    inputs = decompose_input(up.problem, o)
    for (i,(prob_key, b)) in enumerate(b_old)
        b_vec = [b[:,i] for i=1:up.up.k]
        b_stacked = hcat(update(up.up, b_vec, a, inputs[i][2][:])...)
        push!(b_new, (prob_key, b_stacked))
    end
    return b_new
end

struct PreviousObsDecUpdater <: Updater
    up::PreviousObservationUpdater
    problem::UrbanPOMDP
end

const PreviousObsDecBelief = Vector{Tuple{Symbol, Array{Float64, 1}}}

function initialize_dec_belief(up::PreviousObsDecUpdater, pomdp::UrbanPOMDP, rng::AbstractRNG)
    s = initialstate(pomdp, rng)
    o = initialobs(pomdp, s, rng)
    inputs = decompose_input(pomdp, o)
    return inputs
end
function initialize_dec_belief(up::PreviousObsDecUpdater, pomdp::UrbanPOMDP, s::UrbanState, rng::AbstractRNG)
    o = initialobs(pomdp, s, rng)
    inputs = decompose_input(pomdp, o)
    return inputs
end

function POMDPs.update(up::PreviousObsDecUpdater, b_old::PreviousObsDecBelief, a::UrbanAction, o::Array{Float64})
    b_new = Tuple{Symbol, Array{Float64, 1}}[]
    inputs = decompose_input(up.problem, o)
    for (i,(prob_key, b)) in enumerate(b_old)
        push!(b_new, (prob_key, update(up.up, b, a, inputs[i][2][:])))
    end
    return b_new
end

struct DecomposedPolicy{A} <: Policy
    problem::POMDP
    problem_map::Dict
    policy_map::Dict
    fusion
    action_map::Vector{A}
end

# some default fusion function
const SUM = x -> sum(x)
const MIN = x -> minimum(hcat(l), 2)

function POMDPs.action(policy::DecomposedPolicy{A}, beliefs::Union{KMarkovDecBelief, PreviousObsDecBelief}) where A
    val = value(policy, beliefs)
    @assert length(val) == n_actions(policy.problem)
    return policy.action_map[argmax(val)]
end

function POMDPs.value(policy::DecomposedPolicy{A}, beliefs::Union{KMarkovDecBelief, PreviousObsDecBelief}) where A
    n = length(beliefs)
    val = fill(zeros(n_actions(policy.problem),1), n)
    for (i, (prob_key, b)) in enumerate(beliefs)
        val[i] = value(policy.policy_map[prob_key], b)
    end
    return policy.fusion(val)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# Autogenerated wrapper script for CUDA_full_jll for aarch64-linux-gnu
JLLWrappers.@generate_wrapper_header("CUDA_full")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("../matplotlib_settings.jl")

using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.SymmetricLoop

using ChargedParticleDynamics.GuidingCenter4d: plot_integral_error, plot_loop, plot_trajectories

include("guiding_center_4d_settings_vprk.jl")
include("guiding_center_4d_poincare_invariant_1st.jl")