text:
"""
  transformer_encoder(planes, depth, nheads; mlp_ratio = 4.0, dropout = 0.)

Transformer as used in the base ViT architecture.
([reference](https://arxiv.org/abs/2010.11929)).

# Arguments
- `planes`: number of input channels
- `depth`: number of attention blocks
- `nheads`: number of attention heads
- `mlp_ratio`: ratio of MLP layers to the number of input channels
- `dropout`: dropout rate
"""
function transformer_encoder(planes, depth, nheads; mlp_ratio = 4.0, dropout = 0.)
  layers = [Chain(SkipConnection(prenorm(planes, MHAttention(planes, nheads; attn_drop = dropout,
                                                             proj_drop = dropout)), +),
                  SkipConnection(prenorm(planes, mlp_block(planes, floor(Int, mlp_ratio * planes);
                                                           dropout)), +))
            for _ in 1:depth]
  Chain(layers)
end

"""
    vit(imsize::Dims{2} = (256, 256); inchannels = 3, patch_size::Dims{2} = (16, 16),
        embedplanes = 768, depth = 6, nheads = 16, mlp_ratio = 4.0, dropout = 0.1,
        emb_dropout = 0.1, pool = :class, nclasses = 1000)

Creates a Vision Transformer (ViT) model.
([reference](https://arxiv.org/abs/2010.11929)).

# Arguments
- `imsize`: image size
- `inchannels`: number of input channels
- `patch_size`: size of the patches
- `embedplanes`: the number of channels after the patch embedding
- `depth`: number of blocks in the transformer
- `nheads`: number of attention heads in the transformer
- `mlpplanes`: number of hidden channels in the MLP block in the transformer
- `dropout`: dropout rate
- `emb_dropout`: dropout rate for the positional embedding layer
- `pool`: pooling type, either :class or :mean
- `nclasses`: number of classes in the output
"""
function vit(imsize::Dims{2} = (256, 256); inchannels = 3, patch_size::Dims{2} = (16, 16),
             embedplanes = 768, depth = 6, nheads = 16, mlp_ratio = 4.0, dropout = 0.1,
             emb_dropout = 0.1, pool = :class, nclasses = 1000)

  @assert pool in [:class, :mean]
  "Pool type must be either :class (class token) or :mean (mean pooling)"
  npatches = prod(imsize .÷ patch_size)
  return Chain(Chain(PatchEmbedding(imsize; inchannels, patch_size, embedplanes),
                     ClassTokens(embedplanes),
                     ViPosEmbedding(embedplanes, npatches + 1),
                     Dropout(emb_dropout),
                     transformer_encoder(embedplanes, depth, nheads; mlp_ratio, dropout),
                     (pool == :class) ? x -> x[:, 1, :] : seconddimmean),
               Chain(LayerNorm(embedplanes), Dense(embedplanes, nclasses, tanh_fast)))
end

vit_configs = Dict(:tiny => (depth = 12, embedplanes = 192, nheads = 3),
                   :small => (depth = 12, embedplanes = 384, nheads = 6),
                   :base => (depth = 12, embedplanes = 768, nheads = 12),
                   :large => (depth = 24, embedplanes = 1024, nheads = 16),
                   :huge => (depth = 32, embedplanes = 1280, nheads = 16),
                   :giant => (depth = 40, embedplanes = 1408, nheads = 16, mlp_ratio = 48/11),
                   :gigantic => (depth = 48, embedplanes = 1664, nheads = 16, mlp_ratio = 64/13))

"""
    ViT(mode::Symbol = base; imsize::Dims{2} = (256, 256), inchannels = 3,
        patch_size::Dims{2} = (16, 16), pool = :class, nclasses = 1000)

Creates a Vision Transformer (ViT) model.
([reference](https://arxiv.org/abs/2010.11929)).

# Arguments
- `mode`: the model configuration, one of [:tiny, :small, :base, :large, :huge, :giant, :gigantic]
- `imsize`: image size
- `inchannels`: number of input channels
- `patch_size`: size of the patches
- `pool`: pooling type, either :class or :mean
- `nclasses`: number of classes in the output

See also [`Metalhead.vit`](#).
"""
struct ViT
  layers
end

function ViT(mode::Symbol = :base; imsize::Dims{2} = (256, 256), inchannels = 3,
             patch_size::Dims{2} = (16, 16), pool = :class, nclasses = 1000)
  @assert mode in keys(vit_configs) "`mode` must be one of $(keys(vit_configs))"
  kwargs = vit_configs[mode]
  layers = vit(imsize; inchannels, patch_size, nclasses, pool, kwargs...)

  ViT(layers)
end

(m::ViT)(x) = m.layers(x)

backbone(m::ViT) = m.layers[1]
classifier(m::ViT) = m.layers[2]

@functor ViT


text:
using HTTP, Logging
import Sockets

const Token = UInt64

struct Connection
    ws
    api_addr
end

abstract type WSAPIService <: Plugin end
mutable struct WSAPIServiceImpl<: WSAPIService
    manager_factory
    connections::Dict{Token, Connection}
    helper
    socket
    WSAPIServiceImpl(; wsapi_managerfactory, options...) = new(wsapi_managerfactory, Dict())
end
Plugins.symbol(plugin::WSAPIServiceImpl) = :wsapi
__init__() = Plugins.register(WSAPIServiceImpl)

mutable struct WSAPIHelper <: Actor{Any}
    manager_addr::Addr
    wsapi::WSAPIServiceImpl
    core
    WSAPIHelper(manager_addr, wsapi) = new(manager_addr, wsapi)
end

function spawn_helper(wsapi::WSAPIServiceImpl, scheduler)
    manager = wsapi.manager_factory(typeof(emptycore(scheduler.service)))
    manager_addr = spawn(scheduler.service, manager)
    wsapi.helper = WSAPIHelper(manager_addr, wsapi)
    spawn(scheduler.service, wsapi.helper)
end

Circo.schedule_start(wsapi::WSAPIServiceImpl, scheduler::AbstractScheduler{TMsg}) where {TMsg} = begin
    spawn_helper(wsapi, scheduler)
    listenport = 7274 + port(postcode(scheduler)) - CircoCore.PORT_RANGE[1] # PAPI, Programmable API
    ipaddr = Sockets.IPv4(0) # TODO config
    try
        wsapi.socket = Sockets.listen(Sockets.InetAddr(ipaddr, listenport))
        @info "Programmable API websocket listening on $(ipaddr):$(listenport)"
    catch e
        @warn "Programmable API unable to listen on $(ipaddr):$(listenport)", e
    end
    @async HTTP.listen(ipaddr, listenport; server=wsapi.socket) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http; binary=true) do ws
                try
                    @debug "Programmable API got WS connection", ws
                    handle_connection(wsapi, ws, scheduler)
                catch e
                    @info "wsapi: $e"
                end
            end
        end
    end
end

function Circo.schedule_stop(wsapi::WSAPIServiceImpl, scheduler)
    for conn in values(wsapi.connections)
        try
            close(conn.ws)
        catch e end
    end
    isdefined(wsapi, :socket) && close(wsapi.socket)
end

function create_api(wsapi::WSAPIServiceImpl, ws, scheduler)
    token = rand(ProgrammableAPI.PAPIReqId)
    conn = Connection(ws, nothing)
    wsapi.connections[token] = conn
    send(scheduler.service, wsapi.helper, wsapi.helper.manager_addr,
        CreateEngine(token, addr(wsapi.helper)))
    return token
end

function handle_connection(wsapi::WSAPIServiceImpl, ws, scheduler)
    @debug "Programmable API ws handle_connection on thread $(Threads.threadid())"
    token = create_api(wsapi, ws, scheduler)
    readloop(wsapi, token, scheduler)
end

Circo.onmessage(me::WSAPIHelper, msg::EngineCreated, service) = begin
    connection = get(me.wsapi.connections, msg.token, nothing)
    if isnothing(connection)
        @info "Got EngineCreated with unknown token $(msg.token)"
        return nothing
    end
    ws = connection.ws
    conn = Connection(ws, msg.addr)
    me.wsapi.connections[msg.token] = conn
    return nothing
end

function sendws(rawmsg, ws)
    try
        write(ws, rawmsg)
    catch e
        @error "Unable to write to websocket." exception=(e, catch_backtrace())
    end
end

function handleinput(wsapi::WSAPIServiceImpl, buf, connection, scheduler)
    send(scheduler.service, wsapi.helper, connection.api_addr,
        PAPIInput(addr(wsapi.helper), buf))
    return nothing
end

function _close_connection(wsapi, token, service)
    try
        connection = get(wsapi.connections, token, nothing)
        isnothing(connection) && return nothing
        # TODO send(service, wsapi.helper, connection.api_addr, DestroyEngine(connection.token))
        delete!(wsapi.connections, token)
        isnothing(connection.ws) && return nothing
        close(connection.ws)
    catch e
        @debug "Error closing websocket", e
    end
    return nothing
end

function readloop(wsapi::WSAPIServiceImpl, token, scheduler)
    connection = wsapi.connections[token]
    wait_count = 0
    while isnothing(connection.api_addr) && wait_count < 1000
        sleep(0.001)
        wait_count += 1
        connection = wsapi.connections[token]
    end
    if isnothing(connection.api_addr)
        @error "API endpoint creation timeouted."
        return _close_connection(wsapi, token, scheduler.service)
    end
    @debug "Created API endpoint at $(connection.api_addr)"
    ws = connection.ws
    buf = nothing
    try
        while !eof(ws)
            try
                buf = readavailable(ws)
            catch e
                @debug "Websocket closed: $e"
                return
            end
            handleinput(wsapi, buf, connection, scheduler)
        end
    catch e
        if e isa MethodError && e.f == convert
            @show e
            @info "Erroneous websocket frame: ", buf
        else
            # TODO this causes segfault on 1.5.0 with multithreading
            if Threads.nthreads() == 1
                @error "Exception while handling websocket frame" exception=(e, catch_backtrace())
            else
                @error "Exception while handling websocket frame: $e"
                @error "Cannot print stack trace due to an unknown issue in Base or HTTP.jl. Rerun with JULIA_NUM_THREADS=1 to get more info"
            end
        end
    end
    @debug "PAPI Websocket closed", ws
end

Circo.onmessage(me::WSAPIHelper, msg::PAPIOutput, service) = begin
    try
        ws = me.wsapi.connections[msg.token].ws
        write(ws, msg.body)
        return nothing
    catch e
        @debug "Error writing to websocket", e
    end
    _close_connection(me.wsapi, msg.token, service)
end

ext: jl


text:
# REQUIRE
#   Distributions
#   PyCall
#   SkLearn

export GaussianMixtureRegressionDriver

# Based on Comparison of Parametric and Non-Parametric Approaches for Vehicle Speed Prediction
type GaussianMixtureRegressionDriver{A,F<:AbstractFeatureExtractor} <: DriverModel{A, IntegratedContinuous}
    rec::SceneRecord
    context::IntegratedContinuous

    gmr::GMR
    extractor::F
    features::Vector{Float64}
    action::Vector{Float64}
end

get_name(::GaussianMixtureRegressionDriver) = "GMR"
action_context(driver::GaussianMixtureRegressionDriver) = driver.context
function reset_hidden_state!(driver::GaussianMixtureRegressionDriver)
    empty!(driver.rec)
    driver
end
function observe!{A,F}(driver::GaussianMixtureRegressionDriver{A,F}, scene::Scene, roadway::Roadway, egoid::Int)
    update!(driver.rec, scene)
    vehicle_index = get_index_of_first_vehicle_with_id(scene, egoid)
    pull_features!(driver.extractor, driver.features, driver.rec, roadway, vehicle_index)
    driver.gmr(driver.features)
    driver
end
Base.rand{A,F}(driver::GaussianMixtureRegressionDriver{A,F}) = convert(A, rand(driver.gmr.mixture_Act_given_Obs))
Distributions.pdf{A,F}(driver::GaussianMixtureRegressionDriver{A,F}, a::A) = pdf(driver.gmr.mixture_Act_given_Obs, convert(Vector{Float64}, a))
Distributions.logpdf{A,F}(driver::GaussianMixtureRegressionDriver{A,F}, a::A) = logpdf(driver.gmr.mixture_Act_given_Obs, convert(Vector{Float64}, a))


text:
using KernelFunctions
using MLDataUtils
using Zygote
using Flux
using Distributions, LinearAlgebra
using Plots

Flux.@functor SqExponentialKernel
Flux.@functor ScaleTransform
Flux.@functor KernelSum
Flux.@functor Matern32Kernel

xmin = -3; xmax = 3
x = range(xmin,xmax,length=100)
x_test = range(xmin,xmax,length=300)
x,y = noisy_function(sinc,x;noise=0.1)
X = reshape(x,:,1)
X_test = reshape(x_test,:,1)
k = SqExponentialKernel(1.0)#+Matern32Kernel(2.0)
λ = [-1.0]
f(x,k,λ) = kernelmatrix(k,x,X,obsdim=1)*inv(kernelmatrix(k,X,obsdim=1)+exp(λ[1])*I)*y
f(X,k,1.0)
loss(k,λ) = f(X,k,λ) |> ŷ ->sum(y-ŷ)/length(y)+exp(λ[1])*norm(ŷ)
loss(k,λ)
ps = Flux.params(k)
push!(ps,λ)
opt = Flux.Momentum(0.1)
##
for i in 1:10
    grads = Zygote.gradient(()->loss(k,λ),ps)
    Flux.Optimise.update!(opt,ps,grads)
    p = Plots.scatter(x,y,lab="data",title="Loss = $(loss(k,λ))")
    Plots.plot!(x_test,f(X_test,k,λ),lab="Prediction",lw=3.0)
    display(p)
end


text:
using HyperGraphs

x = HyperGraph([HyperEdge([1, 2, 3]), HyperEdge([2, 3, 4])])
srcs = [[1], [4]]; tgts = [[1, 2, 3], [2, 3, 4]]
chx = ChemicalHyperGraph([ChemicalHyperEdge(src, tgt) for (src, tgt) in zip(srcs, tgts)])

# num_has_vertex
e = HyperEdge([1, 2, 3, 1])
@test HyperGraphs.num_has_vertex(e, 1) == 2
@test HyperGraphs.num_has_vertex(e, 4) == 0
che = ChemicalHyperEdge([1], [1, 2, 3])
@test HyperGraphs.num_has_vertex(che, 1) == 2
@test HyperGraphs.num_has_vertex(che, 4) == 0
@test HyperGraphs.num_has_vertex(x, 2) == 2
@test HyperGraphs.num_has_vertex(chx, 2) == 2

# vertices_to_indices
@test all([HyperGraphs.vertices_to_indices(x)[i] == i for i in 1:nv(x)])
@test all([HyperGraphs.vertices_to_indices(chx)[i] == i for i in 1:nv(chx)])

# stoichiometries
@test HyperGraphs.stoichiometries(ChemicalHyperEdge([:X], [:X]), :X) == ([1], [1])

# get_hyperedge_type
@test HyperGraphs.get_hyperedge_type(HyperGraph) == HyperEdge
@test HyperGraphs.get_hyperedge_type(HyperGraph{String}) == HyperEdge
@test HyperGraphs.get_hyperedge_type(ChemicalHyperGraph) == ChemicalHyperEdge
@test HyperGraphs.get_hyperedge_type(ChemicalHyperGraph{Symbol}) == ChemicalHyperEdge


text:
using OneHot, Test, Random
using OneHot: argmax_

@testset "argmax_" begin
    @inferred argmax_(randn(3,2,4); dims=2)
    @test size(argmax_(randn(3,2,4); dims=2)) == (3,4)
    A = [3 4; 1 2]
    i = argmax_(A; dims=2)
    @test A[i] == [4, 2]
end

@testset "decode" begin
    A = randn(5,5,4)
    cl = OneHot.decode(A)
    for I in CartesianIndices(cl)
        @test maximum(A[:,I]) == A[cl[I],I]
    end
    @test size(cl) == size(A)[2:end]

    X = OneHotArray(rand(1:4,5,5), 4)
    cl = OneHot.decode(X)
    for I in CartesianIndices(cl)
        @test maximum(X[:,I]) == X[cl[I],I]
    end
end

q = 10
A = rand(1:q, 4,5,6)
X = @inferred OneHot.encode(A, 1:q)
@test size(X) == (q, size(A)...)
@test vec(X) == [X[i] for i=1:length(X)]
for i in CartesianIndices(A), a = 1:q
    @test X[a,i] == (A[i] == a)
end
@inferred OneHot.decode(X)
@test OneHot.decode(X) == A
@test OneHot.decode(Array(X)) == A
@inferred OneHot.decode(Array(X))


text:
module DoubleDouble

export DoubleF64, ComplexDF64

import Printf
import Base: +, -, *, /, ^, <, ==, <=

# Basic building blocks
"""
    quick_two_sum(a, b)

Computes `s = fl(a+b)` and `e = err(a+b)`. Assumes `|a| ≥ |b|`.
"""
@inline function quick_two_sum(a::Float64, b::Float64)
    s = a + b
    e = b - (s - a)
    s, e
end

"""
    two_sum(a, b)

Computes `s = fl(a+b)` and `e = err(a+b)`.
"""
@inline function two_sum(a::Float64, b::Float64)
    s = a + b
    v = s - a
    e = (a - (s - v)) + (b - v)

    s, e
end

"""
    two_diff(a, b)

Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff(a::Float64, b::Float64)
    s = a - b
    v = s - a
    e = (a - (s - v)) - (b + v)

    s, e
end

"""
    two_prod(a, b)

Computes `s = fl(a*b)` and `e = err(a*b)`.
"""
@inline function two_prod(a::Float64, b::Float64)
    p = a * b
    e = fma(a, b, -p)
    p, e
end

"""
    two_square(a)

Computes `s = fl(a*a)` and `e = err(a*a)`. Faster than [`two_prod(a, a)`](@ref).
"""
@inline function two_square(a::Float64)
    p = a * a
    e = fma(a, a, -p)
    p, e
end

"""
    DoubleF64(x [, mode::ComputeMode]) <: AbstractFloat
"""
struct DoubleF64 <: AbstractFloat
    hi::Float64
    lo::Float64
end

DoubleF64(x::DoubleF64) = DoubleF64(x.hi, x.lo)
function DoubleF64(x::Float64)
    DoubleF64(x, isinf(x) ? Inf : 0.0)
end

function DoubleF64(x::Float32)
    DoubleF64(convert(Float64, x), isinf(x) ? Inf : 0.0)
end
function DoubleF64(x::Float16)
    DoubleF64(convert(Float64, x), isinf(x) ? Inf : 0.0)
end
function DoubleF64(x::Integer)
    DoubleF64(convert(Float64, x), isinf(x) ? Inf : 0.0)
end
function DoubleF64(x::BigFloat)
    z = convert(Float64, x)
    DoubleF64(z, convert(Float64, x - z))
end
function DoubleF64(x::Irrational)
    DoubleF64(big(x))
end
function DoubleF64(x::Rational)
    DoubleF64(numerator(x)) / DoubleF64(denominator(x))
end

const ComplexDF64 = Complex{DoubleF64}

hi(x::DoubleF64) = x.hi
lo(x::DoubleF64) = x.lo

Base.isbits(::DoubleF64) = true
Base.isbitstype(::Type{DoubleF64}) = true

Base.zero(::DoubleF64) = DoubleF64(0.0, 0.0)
Base.zero(::Type{DoubleF64}) = DoubleF64(0.0, 0.0)

Base.one(::DoubleF64) = DoubleF64(1.0, 0.0)
Base.one(::Type{DoubleF64}) = DoubleF64(1.0, 0.0)

Base.convert(::Type{T}, a::DoubleF64) where {T<:AbstractFloat} = convert(T, a.hi)
Base.convert(::Type{BigFloat}, a::DoubleF64) = big(a.hi) + big(a.lo)
Base.BigFloat(a::DoubleF64) = convert(BigFloat, a)
Base.convert(::Type{T}, a::DoubleF64) where {T<:Integer} = convert(T, a.hi)
Base.convert(::Type{Integer}, a::DoubleF64) = convert(Int64, a.hi)
Base.convert(::Type{BigInt}, a::DoubleF64) = convert(BigInt, big(a.hi) + big(a.lo))

Base.convert(::Type{DoubleF64}, x::DoubleF64) = x
Base.convert(::Type{DoubleF64}, x::AbstractFloat) = DoubleF64(x)
Base.convert(::Type{DoubleF64}, x::Irrational) = DoubleF64(x)
Base.convert(::Type{DoubleF64}, x::Integer) = DoubleF64(x)

Base.promote_rule(::Type{DoubleF64}, ::Type{<:Integer}) = DoubleF64
Base.promote_rule(::Type{DoubleF64}, ::Type{BigInt}) = BigFloat
Base.promote_rule(::Type{DoubleF64}, ::Type{BigFloat}) = BigFloat
Base.promote_rule(::Type{DoubleF64}, ::Type{Float64}) = DoubleF64
Base.promote_rule(::Type{DoubleF64}, ::Type{Float32}) = DoubleF64
Base.promote_rule(::Type{DoubleF64}, ::Type{Float16}) = DoubleF64

Base.big(x::DoubleF64) = big(x.hi) + big(x.lo)

const double_π = DoubleF64(pi)
const double_2pi = DoubleF64(2 * BigFloat(Base.pi))
const double_pi = DoubleF64(BigFloat(Base.pi))
const double_pi2 = DoubleF64(BigFloat(Base.pi) * 0.5)
const double_pi4 = DoubleF64(BigFloat(Base.pi) * 0.25)
const double_pi16 = DoubleF64(BigFloat(Base.pi) * (1 / 16))
const double_3pi4 = DoubleF64(BigFloat(Base.pi) * 0.75)
const double_nan = DoubleF64(NaN, NaN)
const double_inf = DoubleF64(Inf)
const double_e = convert(DoubleF64, Base.MathConstants.e)

const double_log2 = convert(DoubleF64, log(BigFloat(2.0)))
const double_log10 = convert(DoubleF64, log(BigFloat(10.0)))

const double_eps = 4.93038065763132e-32 # 2^-104

#
# ADDITION
#
"""
    wide_add(a::Float64, b::Float64)

Add two `Float64`s with `DoubleF64` precision.
"""
@inline function wide_add(a::Float64, b::Float64)
    hi, lo = two_sum(a, b)

    DoubleF64(hi, lo)
end

@inline function +(a::DoubleF64, b::Float64)
    hi, lo = two_sum(a.hi, b)
    lo += a.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleF64(hi, lo)
end
+(a::Float64, b::DoubleF64) = b + a
+(a::Integer, b::DoubleF64) = b + float(a)
+(a::DoubleF64, b::Integer) = a + float(b)

@inline function +(a::DoubleF64, b::DoubleF64)
    hi, lo = two_sum(a.hi, b.hi)
    lo += (a.lo + b.lo)
    hi, lo = quick_two_sum(hi, lo)

    DoubleF64(hi, lo)
end

#
# SUBSTRACTION
#
"""
    wide_sub(a::Float64, b::Float64)

Subtract two `Float64`s with `DoubleF64` precision.
"""
@inline function wide_sub(a::Float64, b::Float64)
    hi, lo = two_diff(a, b)

    DoubleF64(hi, lo)
end

@inline function -(a::DoubleF64, b::Float64)
    hi, lo = two_diff(a.hi, b)
    lo += a.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleF64(hi, lo)
end

@inline function -(a::Float64, b::DoubleF64)
    hi, lo = two_diff(a, b.hi)
    lo -= b.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleF64(hi, lo)
end

@inline function -(a::DoubleF64, b::DoubleF64)
    hi, lo = two_diff(a.hi, b.hi)
    lo += a.lo
    lo -= b.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleF64(hi, lo)
end
-(a::DoubleF64) = DoubleF64(-a.hi, -a.lo)

#
# MULTIPLICATION
#
"""
    wide_mutiply(a::Float64, b::Float64)

Multiply two `Float64`s with `DoubleF64` precision.
"""
@inline function wide_mul(a::Float64, b::Float64)
    hi, lo = two_prod(a, b)
    DoubleF64(hi, lo)
end

@inline function *(a::DoubleF64, b::Float64)
    p1, p2 = two_prod(a.hi, b)
    p2 += a.lo * b
    p1, p2 = quick_two_sum(p1, p2)

    DoubleF64(p1, p2)
end
*(a::Float64, b::DoubleF64) = b * a
*(a::Integer, b::DoubleF64) = b * float(a)
*(a::Bool, b::DoubleF64) = b * convert(Float64, a)
*(a::DoubleF64, b::Integer) = a * float(b)

@inline function *(a::DoubleF64, b::DoubleF64)
    p1, p2 = two_prod(a.hi, b.hi)
    p2 += a.hi * b.lo + a.lo * b.hi
    p1, p2 = quick_two_sum(p1, p2)

    DoubleF64(p1, p2)
end

#
# DIVISION
#
"""
    wide_div(a::Float64, b::Float64)

Divide two `Float64`s with `DoubleF64` precision.
"""
@inline function wide_div(a::Float64, b::Float64)
    q1 = a / b
    # Compute a - q1 * b
    p1, p2 = two_prod(q1, b)
    s, e = two_diff(a, p1)
    e -= p2

    # get next approximation
    q2 = (s + e) / b

    s, e = quick_two_sum(q1, q2)

    DoubleF64(s, e)
end

@inline function /(a::DoubleF64, b::Float64)
    q1 = a.hi / b

    # Compute  this - q1 * d
    p1, p2 = two_prod(q1, b)
    s, e = two_diff(a.hi, p1)
    e += a.lo
    e -= p2

    # get next approximation.
    q2 = (s + e) / b

    # renormalize
    hi, lo = quick_two_sum(q1, q2)

    DoubleF64(hi, lo)
end

@inline function /(a::DoubleF64, b::DoubleF64)
    q1 = a.hi / b.hi  # approximate quotient

    # compute  this - q1 * dd
    r = b * q1
    s1, s2 = two_diff(a.hi, r.hi)
    s2 -= r.lo
    s2 += a.lo

    # get next approximation
    q2 = (s1 + s2) / b.hi

    # renormalize
    hi, lo = quick_two_sum(q1, q2)

    DoubleF64(hi, lo)
end

/(a::Float64, b::DoubleF64) = DoubleF64(a) / b
/(a::Integer, b::DoubleF64) = float(a) / b
/(a::DoubleF64, b::Integer) = a / float(b)

Base.inv(a::DoubleF64) = 1.0 / a

Base.rem(a::DoubleF64, b::DoubleF64) = a - round(a / b) * b
@inline function Base.divrem(a::DoubleF64, b::DoubleF64)
    n = round(a / b)
    n, a - n * b
end

function Base.mod(x::DoubleF64, y::DoubleF64)
    n = round(a / b)
    return (a - b * n)
end

#
# POWERS
#

"""
    square(x::DoubleF64)

Compute `x * x` in a more efficient way.
"""
@inline function square(a::DoubleF64)
    p1, p2 = two_square(a.hi)
    p2 += 2.0 * a.hi * a.lo
    p2 += a.lo * a.lo
    hi, lo = quick_two_sum(p1, p2)
    return DoubleF64(hi, lo)
end
square(a::Float64) = wide_square(a)
"""
    wide_square(x::Float64)

Convert `x` to a DoubleF64 and then compute `x*x`.
"""
@inline function wide_square(a::Float64)
    hi, lo = two_square(a)
    return DoubleF64(hi, lo)
end
# Implementation adapted from Base
@inline function power_by_squaring(x::DoubleF64, p::Integer)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return square(x)
    end
    P = abs(p)
    t = trailing_zeros(P) + 1
    P >>= t
    while (t -= 1) > 0
        x = square(x)
    end
    y = x
    while P > 0
        t = trailing_zeros(P) + 1
        P >>= t
        while (t -= 1) >= 0
            x = square(x)
        end
        y *= x
    end
    if (p < 0)
        return 1.0 / y
    end
    y
end

^(a::DoubleF64, p::Integer) = power_by_squaring(a, p)
^(a::DoubleF64, b::Real) = iszero(a) ? one(a) : exp(b * log(a))

@inline function Base.sqrt(a::DoubleF64)
    # Strategy:  Use Karp's trick:  if x is an approximation
    # to sqrt(a), then
    #
    #  sqrt(a) = a*x + [a - (a*x)^2] * x / 2   (approx)
    #
    # The approximation is accurate to twice the accuracy of x.
    # Also, the multiplication (a*x) and [-]*x can be done with
    # only half the precision.
    if a.hi < 0
        throw(DomainError("sqrt will only return a complex result if called with a complex argument."))
    end
    if iszero(a)
        return zero(a)
    end

    x = inv(sqrt(a.hi))
    ax = a.hi * x

    wide_add(ax, (a - square(ax)).hi * (x * 0.5))
end

"""
    wide_sqrt(x::Float64)

Convert `x` to a DoubleF64 and then compute `sqrt(x)`.
"""
wide_sqrt(a::Float64) = sqrt(DoubleF64(a))

# COMPARISON & EQUALITY
<(a::DoubleF64, b::DoubleF64) = a.hi + a.lo < b.hi + b.lo
<(a::DoubleF64, b::Float64) = a.hi < b || (a.hi == b) && a.lo < 0.0
<(a::Float64, b::DoubleF64) = a < b.hi || (a == b.hi) && b.lo > 0.0

Base.isless(a::DoubleF64, b::DoubleF64) = isless(a.hi + a.lo, b.hi + b.lo)

<=(a::DoubleF64, b::DoubleF64) = !(b < a)
<=(a::DoubleF64, b::Float64) = !(b < a)
<=(a::Float64, b::DoubleF64) = !(b < a)

==(a::DoubleF64, b::Float64) = a.hi == b && a.lo == 0.0
==(a::Float64, b::DoubleF64) = b == a
Base.iszero(a::DoubleF64) = a.hi == 0.0

Base.isone(a::DoubleF64) = a.hi == 1.0 && a.lo == 0.0
Base.abs(a::DoubleF64) = a.hi < 0.0 ? -a : a

Base.eps(::DoubleF64) = 4.93038065763132e-32 # 2^-104
Base.eps(::Type{DoubleF64}) = 4.93038065763132e-32 # 2^-104

# Base.realmin(::Type{DoubleF64}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
# Base.realmin(::Type{DoubleF64}) = 2.0041683600089728e-292 # = 2^(-1022 + 53)
# Base.realmax(::Type{DoubleF64}) = DoubleF64(1.79769313486231570815e+308, 9.97920154767359795037e+291);
# Base.realmax(::Type{DoubleF64}) = DoubleF64(1.79769313486231570815e+308, 9.97920154767359795037e+291);

Base.isnan(a::DoubleF64) = isnan(a.hi) || isnan(a.lo)
Base.isinf(a::DoubleF64) = isinf(a.hi)
Base.isfinite(a::DoubleF64) = isfinite(a.hi)

#
# ROUNDING
#
@inline function Base.round(a::DoubleF64, r::RoundingMode = RoundNearest)
    hi = round(a.hi, r)
    lo = 0.0

    if hi == a.hi
        # High word is an integer already.  Round the low word.
        lo = round(a.lo, r)

        # Renormalize. This is needed if hi = some integer, lo = 1/2.
        hi, lo = quick_two_sum(hi, lo)
    else
        # High word is not an integer.
        if abs(hi - a.hi) == 0.5 && a.lo < 0.0
            # There is a tie in the high word, consult the low word to break the tie.
            hi -= 1.0
        end
    end

    DoubleF64(hi, lo)
end

@inline function Base.floor(a::DoubleF64)
    hi = floor(a.hi)
    lo = 0.0

    if hi == a.hi
        lo = floor(a.lo)
        hi, lo = quick_two_sum(hi, lo)
    end

    DoubleF64(hi, lo)
end

@inline function Base.floor(::Type{I}, a::DoubleF64) where {I<:Integer}
    hi = floor(I, a.hi)
    lo = zero(I)

    if hi == a.hi
        lo = floor(I, a.lo)
    end

    hi + lo
end

@inline function Base.ceil(a::DoubleF64)
    hi = ceil(a.hi)
    lo = 0.0

    if hi == a.hi
        lo = ceil(a.lo)
        hi, lo = quick_two_sum(hi, lo)
    end

    DoubleF64(hi, lo)
end

@inline function Base.ceil(::Type{I}, a::DoubleF64) where {I<:Integer}
    hi = ceil(I, a.hi)
    lo = zero(I)

    if hi == a.hi
        lo = ceil(I, a.lo)
    end

    hi + lo
end

Base.trunc(a::DoubleF64) = a.hi ≥ 0.0 ? floor(a) : ceil(a)
Base.trunc(::Type{I}, a::DoubleF64) where {I<:Integer} =
    a.hi ≥ 0.0 ? floor(I, a) : ceil(I, a)
Base.isinteger(x::DoubleF64) = iszero(x - trunc(x))

# import Random: AbstractRNG, GLOBAL_RNG
# function Base.rand(rng::AbstractRNG, S::Type{DoubleF64})
#     u = rand(rng, UInt64)
#     f = Float64(u)
#     uf = UInt64(f)
#     ur = uf > u ? uf - u : u - uf
#     DoubleF64(5.421010862427522e-20 * f, 5.421010862427522e-20 * Float64(ur))
# end
# Base.rand(::Type{DoubleF64}) = rand(GLOBAL_RNG, DoubleF64)
#
# function Base.rand(rng::AbstractRNG, T::Type{<:DoubleF64}, dims::Vararg{Int, N}) where N
#     rands = Array{T}(dims)
#     for l in eachindex(rands)
#         rands[l] = rand(rng, T)
#     end
#     rands
# end
# Base.rand(T::Type{<:DoubleF64}, dims::Vararg{Int, N}) where N = rand(GLOBAL_RNG, T, dims)
#
# Base.rand(::Type{Complex{DoubleF64}}) = rand(Complex{DoubleF64})
# Base.rand(::Type{Complex{DoubleF64}}, dims::Vararg{Int, N}) where N = rand(Complex{DoubleF64}, dims)
# Base.rand(rng::AbstractRNG, ::Type{Complex{DoubleF64}}) = rand(rng, Complex{DoubleF64})
# Base.rand(rng::AbstractRNG, ::Type{Complex{DoubleF64}}, dims::Vararg{Int, N}) where N = rand(rng, Complex{DoubleF64}, dims)

function Base.decompose(a::DoubleF64)::Tuple{Int128,Int,Int}
    hi, lo = a.hi, a.lo
    num1, pow1, den1 = Base.decompose(hi)
    num2, pow2, den2 = Base.decompose(lo)

    num = Int128(num1)

    pdiff = pow1 - pow2
    shift = min(pdiff, 52)
    signed_num = den1 * (Int128(num) << shift) # den1 is +1/-1
    signed_num += den2 * (num2 >> (pdiff - shift)) # den2 is +1/-1

    num = abs(signed_num)
    den = signed_num ≥ 0 ? 1 : -1
    pow = pow1 - shift

    num, pow, den
end

# Precomputed valuee
const inv_fact = [DoubleF64(1.0 / BigFloat(factorial(Int64(k)))) for k = 3:17]
const ninv_fact = length(inv_fact)

Base.ldexp(a::DoubleF64, exp::Int) = DoubleF64(ldexp(a.hi, exp), ldexp(a.lo, exp))
"""
	mul_pwr2(a::DoubleF64, b::Float64)

`a` * `b`,  where `b` is a power of 2.
"""
mul_pwr2(a::DoubleF64, b::Float64) = DoubleF64(a.hi * b, a.lo * b)

function Base.exp(a::DoubleF64)
    # Strategy:  We first reduce the size of x by noting that
    #
    #       exp(kr + m * log(2)) = 2^m * exp(r)^k
    #
    # where m and k are integers.  By choosing m appropriately
    # we can make |kr| <= log(2) / 2 = 0.347.  Then exp(r) is
    # evaluated using the familiar Taylor series.  Reducing the
    # argument substantially speeds up the convergence.

    k = 512.0
    inv_k = 1.0 / k

    if a.hi ≤ -709.0
        return zero(a)
    end

    if a.hi ≥ 709.0
        convert(DoubleF64, Inf)
    end

    if iszero(a)
        return one(a)
    end

    if isone(a)
        return double_e
    end

    m = floor(a.hi / double_log2.hi + 0.5)
    r = mul_pwr2(a - double_log2 * m, inv_k)

    #
    p = square(r)
    s = r + mul_pwr2(p, 0.5)
    p *= r
    t = p * inv_fact[1]

    i = 1
    while true
        s += t
        p *= r
        i += 1
        t = p * inv_fact[i]

        if abs(convert(Float64, t)) < inv_k * double_eps || i > 6
            break
        end
    end

    s += t

    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s = mul_pwr2(s, 2.0) + square(s)
    s += 1.0

    return ldexp(s, convert(Int, m))
end

Base.log10(a::DoubleF64) = log(a) / double_log10
Base.log2(a::DoubleF64) = log(a) / double_log2
function Base.log(a::DoubleF64)
    # Strategy.  The Taylor series for log converges much more
    # slowly than that of exp, due to the lack of the factorial
    # term in the denominator.  Hence this routine instead tries
    # to determine the root of the function
    #
    #     f(x) = exp(x) - a
    #
    # using Newton iteration.  The iteration is given by
    #
    #     x' = x - f(x)/f'(x)
    #        = x - (1 - a * exp(-x))
    #        = x + a * exp(-x) - 1.
    #
    # Only one iteration is needed, since Newton's iteration
    # approximately doubles the number of digits per iteration.

    if isone(a)
        return zero(DoubleF64)
    end

    if a.hi ≤ 0.0
        throw(DomainError("Passed non-positive argument to log."))
    end

    x = DoubleF64(log(a.hi)) # Initial approximation

    x + a * exp(-x) - 1.0
end

const sin_table = [DoubleF64(sin(k * big(π) * 0.0625)) for k = 1:4]
const cos_table = [DoubleF64(cos(k * big(π) * 0.0625)) for k = 1:4]

"""
	sin_taylor(a::DoubleF64)

Computes sin(a) using Taylor series. Assumes |a| <= π/32.
"""
function sin_taylor(a::DoubleF64)
    thresh = 0.5 * abs(convert(Float64, a)) * double_eps

    if iszero(a)
        return zero(a)
    end

    i = 1
    x = -square(a)
    s = a
    r = a
    while true
        r *= x
        t = r * inv_fact[i]
        s += t
        i += 2
        if i > ninv_fact || abs(convert(Float64, t)) < thresh
            break
        end
    end

    s
end

"""
	cos_taylor(a::DoubleF64)

Computes cos(a) using Taylor series. Assumes |a| <= π/32.
"""
function cos_taylor(a::DoubleF64)
    thresh = 0.5 * double_eps

    if iszero(a)
        return one(a)
    end

    x = -square(a)
    r = x
    s = 1.0 + mul_pwr2(r, 0.5)
    i = 2
    while true
        r *= x
        t = r * inv_fact[i]
        s += t
        i += 2
        if i > ninv_fact || abs(convert(Float64, t)) < thresh
            break
        end
    end

    s
end

function sincos_taylor(a::DoubleF64)
    if iszero(a)
        zero(a), one(a)
    end

    sin_a = sin_taylor(a)
    cos_a = sqrt(1.0 - square(sin_a))

    sin_a, cos_a
end

function Base.sin(a::DoubleF64)
    # Strategy.  To compute sin(x), we choose integers a, b so that
    #
    #    x = s + a * (pi/2) + b * (pi/16)
    #
    #  and |s| <= pi/32.  Using the fact that
    #
    #    sin(pi/16) = 0.5 * sqrt(2 - sqrt(2 + sqrt(2)))
    #
    #  we can compute sin(x) from sin(s), cos(s).  This greatly
    #  increases the convergence of the sine Taylor series.

    if iszero(a)
        return zero(a)
    end

    # approximately reduce modulo 2*pi
    z = round(a / double_2pi)
    r = a - double_2pi * z

    # approximately reduce modulo pi/2 and then modulo pi/16.

    q = floor(r.hi / double_pi2.hi + 0.5)
    t = r - double_pi2 * q
    j = convert(Int, q)
    q = floor(t.hi / double_pi16.hi + 0.5)
    t -= double_pi16 * q
    k = convert(Int, q)
    abs_k = abs(k)

    if j < -2 || j > 2
        # Cannot reduce modulo pi/2.
        return double_nan
    end

    if (abs_k > 4)
        # Cannot reduce modulo pi/16.
        return double_nan
    end

    if k == 0
        if j == 0
            return sin_taylor(t)
        elseif j == 1
            return cos_taylor(t)
        elseif j == -1
            return -cos_taylor(t)
        else
            return -sin_taylor(t)
        end
    end

    u = cos_table[abs_k]
    v = sin_table[abs_k]
    sin_t, cos_t = sincos_taylor(t)

    if j == 0
        r = k > 0 ? u * sin_t + v * cos_t : u * sin_t - v * cos_t
    elseif j == 1
        r = k > 0 ? u * cos_t - v * sin_t : u * cos_t + v * sin_t
    elseif j == -1
        r = k > 0 ? v * sin_t - u * cos_t : -u * cos_t - v * sin_t
    else
        r = k > 0 ? -u * sin_t - v * cos_t : v * cos_t - u * sin_t
    end

    r
end

function Base.cos(a::DoubleF64)
    if iszero(a)
        return one(a)
    end

    # approximately reduce modulo 2*pi
    z = round(a / double_2pi)
    r = a - double_2pi * z

    # approximately reduce modulo pi/2 and then modulo pi/16.

    q = floor(r.hi / double_pi2.hi + 0.5)
    t = r - double_pi2 * q
    j = convert(Int, q)
    q = floor(t.hi / double_pi16.hi + 0.5)
    t -= double_pi16 * q
    k = convert(Int, q)
    abs_k = abs(k)

    if j < -2 || j > 2
        # Cannot reduce modulo pi/2.
        return double_nan
    end

    if (abs_k > 4)
        # Cannot reduce modulo pi/16.
        return double_nan
    end

    if k == 0
        if j == 0
            return cos_taylor(t)
        elseif j == 1
            return -sin_taylor(t)
        elseif j == -1
            return sin_taylor(t)
        else
            return -cos_taylor(t)
        end
    end

    u = cos_table[abs_k]
    v = sin_table[abs_k]
    sin_t, cos_t = sincos_taylor(t)

    if j == 0
        r = k > 0 ? u * cos_t - v * sin_t : u * cos_t + v * sin_t
    elseif j == 1
        r = k > 0 ? -u * sin_t - v * cos_t : v * cos_t - u * sin_t
    elseif j == -1
        r = k > 0 ? u * sin_t + v * cos_t : u * sin_t - v * cos_t
    else
        r = k > 0 ? v * sin_t - u * cos_t : -u * cos_t - v * sin_t
    end

    r
end

function Base.sincos(a::DoubleF64)
    if iszero(a)
        return zero(a), one(a)
    end

    # approximately reduce modulo 2*pi
    z = round(a / double_2pi)
    r = a - double_2pi * z

    # approximately reduce modulo pi/2 and then modulo pi/16.
    q = floor(r.hi / double_pi2.hi + 0.5)
    t = r - double_pi2 * q
    j = convert(Int, q)
    abs_j = abs(j)
    q = floor(t.hi / double_pi16.hi + 0.5)
    t -= double_pi16 * q
    k = convert(Int, q)
    abs_k = abs(k)

    if abs_j > 2
        # Cannot reduce modulo pi/2.
        return double_nan, double_nan
    end

    if abs_k > 4
        # Cannot reduce modulo pi/16.
        return double_nan, double_nan
    end

    sin_t, cos_t = sincos_taylor(t)

    if abs_k == 0
        s = sin_t
        c = cos_t
    else
        u = cos_table[abs_k]
        v = sin_table[abs_k]

        s = k > 0 ? u * sin_t + v * cos_t : u * sin_t - v * cos_t
        c = k > 0 ? u * cos_t - v * sin_t : u * cos_t + v * sin_t
    end

    if j == 0
        s, c
    elseif j == 1
        c, -s
    elseif j == -1
        -c, s
    else
        -s, -c
    end
end

Base.atan(a::DoubleF64) = atan(a, one(a))

function Base.atan(y::DoubleF64, x::DoubleF64)
    # Strategy: Instead of using Taylor series to compute
    # arctan, we instead use Newton's iteration to solve
    # the equation
    #
    #    sin(z) = y/r    or    cos(z) = x/r
    #
    # where r = sqrt(x^2 + y^2).
    # The iteration is given by
    #
    #    z' = z + (y - sin(z)) / cos(z)          (for equation 1)
    #    z' = z - (x - cos(z)) / sin(z)          (for equation 2)
    #
    # Here, x and y are normalized so that x^2 + y^2 = 1.
    # If |x| > |y|, then first iteration is used since the
    # denominator is larger.  Otherwise, the second is used.
    if iszero(x)
        if iszero(y)
            double_nan
        end

        return y.hi > 0.0 ? double_pi2 : -double_pi2
    elseif iszero(y)
        return x.hi > 0.0 ? zero(x) : double_pi
    end

    if x == y
        return y.hi > 0.0 ? double_pi4 : -double_3pi4
    end

    if x == -y
        return y.hi > 0.0 ? double_3pi4 : -double_pi4
    end

    r = sqrt(square(x) + square(y))
    xx = x / r
    yy = y / r

    # Compute double precision approximation to atan.
    z = DoubleF64(atan(convert(Float64, y), convert(Float64, x)))

    if abs(xx.hi) > abs(yy.hi)
        # Use Newton iteration 1.  z' = z + (y - sin(z)) / cos(z)
        sin_z, cos_z = sincos(z)
        z += (yy - sin_z) / cos_z
    else
        # Use Newton iteration 2.  z' = z - (x - cos(z)) / sin(z)
        sin_z, cos_z = sincos(z)
        z -= (xx - cos_z) / sin_z
    end

    z
end

function Base.tan(a::DoubleF64)
    s, c = sincos(a)
    s / c
end

function Base.asin(a::DoubleF64)
    abs_a = abs(a)

    if abs_a > 1.0
        throw(DomainError())
    end

    if isone(abs_a)
        return a.hi > 0.0 ? double_pi2 : -double_pi2
    end

    atan(a, sqrt(1.0 - square(a)))
end

function Base.acos(a::DoubleF64)
    abs_a = abs(a)

    if abs_a > 1.0
        throw(DomainError())
    end

    if isone(abs_a)
        return a.hi > 0.0 ? zero(a) : -double_pi
    end

    atan(sqrt(1.0 - square(a)), a)
end

function Base.show(io::IO, x::DoubleF64)
    Printf.@printf io "%.32g" big(x)  # crude approximation to valid number of digits
end

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Declaration of constants.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# This regex match all the ANSI escape sequences that defines decorations.
const _REGEX_ANSI_SEQUENCES = r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])"


text:
exp(π) - π

# Solve x^2 + 5x + 6 = 0 using the quadratic formula (real arithmetics)

# Coefficients a,b,c in ax^2 + bx + c = 0
a = 1
b = 5
c = 6

# The quadratic formula
d = sqrt(b^2 - 4*a*c)
r1 = (-b - d) / 2a
r2 = (-b + d) / 2a
println("The roots are ", r1, " and ", r2)

function myfunc(x,y)
    x + y
end

myfunc(1,2)

myfunc2(x,y) = x + 2y
myfunc2(3,5)

myfunc3 = (x,y) -> x + 3y
myfunc3(3,5)

# Same thing without giving the function a name:
((x,y) -> x + 3y)(3,5)

function mynewfunc(x,y)
    out1 = x + y
    out2 = out1 * (2x + y)
    out1, out2
end

y1, y2 = mynewfunc(2,1)

function real_roots_of_quadratic(a,b,c)
    # Compute the real roots of the quadratic ax^2 + bx + c = 0
    d = sqrt(b^2 - 4*a*c)
    r1 = (-b - d) / 2a
    r2 = (-b + d) / 2a
    r1, r2
end

real_roots_of_quadratic(1, 5, 6)

real_roots_of_quadratic(-1, 5, 6)


text:
#=

  Euler #47 in Julia.

  Problem 47:
  """
  The first two consecutive numbers to have two distinct prime factors are:

  14 = 2 x 7
  15 = 3 x 5

  The first three consecutive numbers to have three distinct
  prime factors are:

  644 = 2^2 x 7 x 23
  645 = 3 x 5 x 43
  646 = 2 x 17 x 19.

  Find the first four consecutive integers to have four distinct primes
  factors. What is the first of these numbers?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.01653001s
function euler47a()
    maxn = 1_000_000
    f = zeros(maxn)
    for i in 2:maxn-1
        if f[i] == 0
            for j in 2*i:i:maxn-1
                f[j] += 1
            end
        end
    end

    goal = [4,4,4,4]
    found = 0
    for i in 1:maxn-3
        if [f[i],f[i+1],f[i+2],f[i+3]] == goal
            found = i
            break
        end
    end
    return found
end

run_euler(euler47a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using NeXLCore
using NeXLMatrixCorrection
using Random
using Test

@testset "XPPU" begin
    rgen = MersenneTwister(0xBADF00D)

    unknown, standard, cxr, e0, θtoa = "K240", "SiO2", n"O K-L3", 20.0e3, deg2rad(40.0)
    unk_mat = material(
        unknown, #
        n"O" => 0.340023,
        n"Mg" => 0.030154,
        n"Si" => 0.186986,
        n"Ti" => 0.059950, #
        n"Zn" => 0.040168,
        n"Zr" => 0.074030,
        n"Ba" => 0.268689,
    )
    std_mat = parse(Material, "SiO2")

    # For comparison with direct calculation...
    zaf = zafcorrection(
        XPP,
        NullFluorescence,
        Coating,
        unk_mat,
        std_mat,
        inner(cxr),
        e0,
        unkCoating = carboncoating(10.0),
        stdCoating = carboncoating(11.0),
    )

    m = NeXLMatrixCorrection.m(XPP, inner(cxr))

    input_uvs = uvs(
        (MassFractionLabel(unknown, elm) => uv(unk_mat[elm], 0.01 * unk_mat[elm]) for elm in keys(unk_mat))...,
        (AtomicWeightLabel(unknown, elm) => uv(a(elm, unk_mat), 0.001 * a(elm, unk_mat)) for elm in keys(unk_mat))...,
        (JzLabel(elm) => Ju(elm) for elm in keys(unk_mat))...,
        NeXLMatrixCorrection.E0Label(unknown) => uv(e0, 0.1e3),
        NeXLMatrixCorrection.mLabel(inner(cxr)) => uv(m, 0.01 * m),
    )

    maintain = MaintainInputs([NeXLMatrixCorrection.mLabel], input_uvs)
    mjz = StepMJZbarb(unknown, [n"O", n"Mg", n"Si", n"Ti", n"Zn", n"Zr", n"Ba"]) | maintain

    mjz_res = mjz(input_uvs)
    mjz_mcres = mcpropagate(mjz, input_uvs, 10000, parallel = false, rng = rgen)

    @test isapprox(value(mjz_res[NeXLMatrixCorrection.BigMLabel(unknown)]), 0.466, atol = 0.001)
    @test isapprox(value(mjz_res[NeXLMatrixCorrection.ZbarbLabel(unknown)]), 22.5, atol = 0.1)
    @test isapprox(value(mjz_res[NeXLMatrixCorrection.E0keVLabel(unknown)]), 20.0, atol = 0.1)
    @test isapprox(value(mjz_res[NeXLMatrixCorrection.JLabel(unknown)]), 0.217, atol = 0.001)
    @test isapprox(σ(mjz_res[NeXLMatrixCorrection.BigMLabel(unknown)]), 0.0023, atol = 0.0001)
    @test isapprox(σ(mjz_res[NeXLMatrixCorrection.ZbarbLabel(unknown)]), 0.25, atol = 0.03)
    @test isapprox(σ(mjz_res[NeXLMatrixCorrection.E0keVLabel(unknown)]), 0.1, atol = 0.01)
    @test isapprox(σ(mjz_res[NeXLMatrixCorrection.JLabel(unknown)]), 0.00198, atol = 0.0001)

    # println("Analytical Result")
    # print(mjz_res)
    # println("MC Result")
    # print(mjz_mcres)

    dpt = StepDPT(unknown, inner(cxr)) | AllInputs()

    dpt_res = dpt(mjz_res)
    dpt_mcres = mcpropagate(dpt ∘ mjz, input_uvs, 1000, parallel = false, rng = rgen)

    @test isapprox(value(dpt_res[NeXLMatrixCorrection.DkLabel(unknown, inner(cxr), 1)]), 6.6e-6, atol = 0.01e-6)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.DkLabel(unknown, inner(cxr), 2)]), 1.49e-05, atol = 0.01e-5)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.DkLabel(unknown, inner(cxr), 3)]), 1.02e-5, atol = 0.01e-5)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.PkLabel(unknown, inner(cxr), 1)]), 0.780, atol = 0.001)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.PkLabel(unknown, inner(cxr), 2)]), 0.100, atol = 0.001)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.PkLabel(unknown, inner(cxr), 3)]), -0.446, atol = 0.001)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.TkLabel(unknown, inner(cxr), 1)]), 0.911, atol = 0.001)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.TkLabel(unknown, inner(cxr), 2)]), 0.231, atol = 0.001)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.TkLabel(unknown, inner(cxr), 3)]), -0.315, atol = 0.001)

    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.DkLabel(unknown, inner(cxr), 1)]), 0.0, atol = 0.001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.DkLabel(unknown, inner(cxr), 2)]), 4.3e-09, atol = 0.1e-9)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.DkLabel(unknown, inner(cxr), 3)]), 9.3e-08, atol = 0.1e-8)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.PkLabel(unknown, inner(cxr), 1)]), 0.0, atol = 0.001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.PkLabel(unknown, inner(cxr), 2)]), 0.0, atol = 0.001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.PkLabel(unknown, inner(cxr), 3)]), 0.00049, atol = 0.00001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.TkLabel(unknown, inner(cxr), 1)]), 0.009, atol = 0.001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.TkLabel(unknown, inner(cxr), 2)]), 0.009, atol = 0.001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.TkLabel(unknown, inner(cxr), 3)]), 0.0088, atol = 0.0001)

    @test isapprox(value(dpt_res[NeXLMatrixCorrection.BigMLabel(unknown)]), 0.466, atol = 0.001)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.ZbarbLabel(unknown)]), 22.5, atol = 0.1)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.E0keVLabel(unknown)]), 20.0, atol = 0.1)
    @test isapprox(value(dpt_res[NeXLMatrixCorrection.JLabel(unknown)]), 0.217, atol = 0.001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.BigMLabel(unknown)]), 0.0023, atol = 0.0001)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.ZbarbLabel(unknown)]), 0.25, atol = 0.03)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.E0keVLabel(unknown)]), 0.1, atol = 0.01)
    @test isapprox(σ(dpt_res[NeXLMatrixCorrection.JLabel(unknown)]), 0.00198, atol = 0.0001)

    # println("Analytical Result")
    # print(dpt_res)
    # println("MC Result")
    # print(dpt_mcres)

    retain = MaintainInputs([NeXLMatrixCorrection.E0keVLabel, NeXLMatrixCorrection.ZbarbLabel], dpt_res)

    qla = StepQlaOoS(unknown, inner(cxr)) | retain
    qla_res = qla(dpt_res)
    qla_model = qla ∘ dpt ∘ mjz
    qla_mcres = mcpropagate(qla_model, input_uvs, 10000, parallel = false, rng = rgen)

    #println("Analytical Result")
    #print(qla_res)
    #println("MC Result")
    #print(qla_mcres)

    @test isapprox(value(qla_res[NeXLMatrixCorrection.QlaLabel(unknown, inner(cxr))]), 0.548, atol = 0.001)
    @test isapprox(value(qla_res[NeXLMatrixCorrection.OoSLabel(unknown, inner(cxr))]), 0.00108, atol = 0.00001)
    @test isapprox(value(qla_res[NeXLMatrixCorrection.ηLabel(unknown)]), 0.253, atol = 0.001)
    @test isapprox(value(qla_res[NeXLMatrixCorrection.JU0Label(unknown)]), 99.8, atol = 0.1)
    @test isapprox(value(qla_res[NeXLMatrixCorrection.WbarLabel(unknown)]), 0.665, atol = 0.001)
    @test isapprox(value(qla_res[NeXLMatrixCorrection.qLabel(unknown)]), 0.987, atol = 0.001)

    @test isapprox(σ(qla_res[NeXLMatrixCorrection.QlaLabel(unknown, inner(cxr))]), 0.0173, atol = 0.001)
    @test isapprox(σ(qla_res[NeXLMatrixCorrection.OoSLabel(unknown, inner(cxr))]), 2.85e-5, atol = 0.1e-5)
    @test isapprox(σ(qla_res[NeXLMatrixCorrection.ηLabel(unknown)]), 0.00218, atol = 0.0001)
    @test isapprox(σ(qla_res[NeXLMatrixCorrection.JU0Label(unknown)]), 0.682, atol = 0.001)
    @test isapprox(σ(qla_res[NeXLMatrixCorrection.WbarLabel(unknown)]), 0.000667, atol = 0.00001)
    @test isapprox(σ(qla_res[NeXLMatrixCorrection.qLabel(unknown)]), 0.00595, atol = 0.0001)

    maintain = MaintainInputs(
        [
            NeXLMatrixCorrection.E0keVLabel,
            NeXLMatrixCorrection.ZbarbLabel,
            NeXLMatrixCorrection.OoSLabel,
            NeXLMatrixCorrection.QlaLabel,
        ],
        qla_res,
    )

    rp = NeXLMatrixCorrection.StepRPhi0(unknown, inner(cxr)) | maintain
    rp_res = rp(qla_res)
    rp_model = rp ∘ qla_model
    rp_mcres = mcpropagate(rp_model, input_uvs, 1000, parallel = false, rng = rgen)

    #println("Analytical Result")
    #print(rp_res)
    #println("MC Result")
    #print(rp_mcres)

    @test isapprox(value(rp_res[NeXLMatrixCorrection.RLabel(unknown, inner(cxr))]), 0.852, atol = 0.001)
    @test isapprox(value(rp_res[NeXLMatrixCorrection.ϕ0Label(unknown, inner(cxr))]), zaf[1].za.ϕ0, atol = 1.0e-7)

    @test isapprox(σ(rp_res[NeXLMatrixCorrection.RLabel(unknown, inner(cxr))]), 0.00146, atol = 0.001)
    @test isapprox(σ(rp_res[NeXLMatrixCorrection.ϕ0Label(unknown, inner(cxr))]), 0.00643, atol = 0.0001)

    maintain = MaintainInputs(
        [NeXLMatrixCorrection.ZbarbLabel, NeXLMatrixCorrection.ϕ0Label, NeXLMatrixCorrection.E0keVLabel],
        rp_res,
    )

    frbar = NeXLMatrixCorrection.StepFRBar(unknown, inner(cxr)) | maintain
    frbar_res = frbar(rp_res)
    frbar_model = frbar ∘ rp_model
    frbar_mcres = mcpropagate(frbar_model, input_uvs, 1000, parallel = false, rng = rgen)

    #println("Analytical Result")
    #print(frbar_res)
    #println("MC Result")
    #print(frbar_mcres)

    @test isapprox(value(frbar_res[NeXLMatrixCorrection.FLabel(unknown, inner(cxr))]), 0.00168, atol = 0.00001)
    @test isapprox(value(frbar_res[NeXLMatrixCorrection.RbarLabel(unknown, inner(cxr))]), 3.31e-4, atol = 0.01e-4)
    @test isapprox(σ(frbar_res[NeXLMatrixCorrection.FLabel(unknown, inner(cxr))]), 2.56e-5, atol = 0.1e-5)
    @test isapprox(σ(frbar_res[NeXLMatrixCorrection.RbarLabel(unknown, inner(cxr))]), 5.36e-6, atol = 0.1e-6)

    @test isapprox(value(frbar_res[NeXLMatrixCorrection.FLabel(unknown, inner(cxr))]), zaf[1].za.F, rtol = 1.0e-7)

    maintain = MaintainInputs(
        [NeXLMatrixCorrection.ϕ0Label, NeXLMatrixCorrection.RbarLabel, NeXLMatrixCorrection.FLabel],
        frbar_res,
    )
    pb = NeXLMatrixCorrection.StepPb(unknown, inner(cxr)) | maintain
    pb_res = pb(frbar_res)
    pb_model = pb ∘ frbar_model
    pb_mcres = mcpropagate(pb_model, input_uvs, 1000, parallel = false, rng = rgen)

    @test isapprox(value(pb_res[NeXLMatrixCorrection.PLabel(unknown, inner(cxr))]), 1.42e4, atol = 0.01e4)
    @test isapprox(value(pb_res[NeXLMatrixCorrection.bLabel(unknown, inner(cxr))]), 7.78e3, atol = 0.01e3)
    @test isapprox(σ(pb_res[NeXLMatrixCorrection.PLabel(unknown, inner(cxr))]), 2.85e2, atol = 0.4e2)
    @test isapprox(σ(pb_res[NeXLMatrixCorrection.bLabel(unknown, inner(cxr))]), 1.25e2, atol = 0.2e2)

    maintain = MaintainInputs(
        [
            NeXLMatrixCorrection.ϕ0Label,
            NeXLMatrixCorrection.FLabel,
            NeXLMatrixCorrection.PLabel,
            NeXLMatrixCorrection.bLabel,
        ],
        pb_res,
    )
    aϵ = NeXLMatrixCorrection.Stepaϵ(unknown, inner(cxr)) | maintain
    aϵ_res = aϵ(pb_res)
    aϵ_model = aϵ ∘ pb_model
    aϵ_mcres = mcpropagate(aϵ_model, input_uvs, 1000, parallel = false, rng = rgen)

    @test isapprox(value(aϵ_res[NeXLMatrixCorrection.aLabel(unknown, inner(cxr))]), 6.75e3, atol = 0.01e3)
    @test isapprox(value(aϵ_res[NeXLMatrixCorrection.aLabel(unknown, inner(cxr))]), zaf[1].za.a, rtol = 1.0e-7)
    @test isapprox(value(aϵ_res[NeXLMatrixCorrection.ϵLabel(unknown, inner(cxr))]), -0.133, atol = 0.001)
    @test isapprox(σ(aϵ_res[NeXLMatrixCorrection.aLabel(unknown, inner(cxr))]), 100, atol = 10)
    @test isapprox(σ(aϵ_res[NeXLMatrixCorrection.ϵLabel(unknown, inner(cxr))]), 0.0013, atol = 0.0001)

    maintain = MaintainInputs(
        [
            NeXLMatrixCorrection.bLabel,
            NeXLMatrixCorrection.ϕ0Label,
            NeXLMatrixCorrection.ϵLabel,
            NeXLMatrixCorrection.FLabel,
        ],
        aϵ_res,
    )
    AB = NeXLMatrixCorrection.StepAB(unknown, inner(cxr)) | maintain
    AB_res = AB(aϵ_res)
    AB_model = AB ∘ aϵ_model
    AB_mcres = mcpropagate(AB_model, input_uvs, 1000, parallel = false, rng = rgen)

    @test isapprox(value(AB_res[NeXLMatrixCorrection.ALabel(unknown, inner(cxr))]), 393, atol = 1)
    @test isapprox(value(AB_res[NeXLMatrixCorrection.BLabel(unknown, inner(cxr))]), -3.79e5, atol = 0.01e5)
    @test isapprox(σ(AB_res[NeXLMatrixCorrection.ALabel(unknown, inner(cxr))]), 8.46, atol = 0.01)
    @test isapprox(σ(AB_res[NeXLMatrixCorrection.BLabel(unknown, inner(cxr))]), 4.0e3, atol = 4.0e3)

    @test isapprox(value(AB_res[NeXLMatrixCorrection.ALabel(unknown, inner(cxr))]), zaf[1].za.A, rtol = 1.0e-7)
    @test isapprox(value(AB_res[NeXLMatrixCorrection.BLabel(unknown, inner(cxr))]), zaf[1].za.B, rtol = 1.0e-7)

    χdata = uvs(
        NeXLMatrixCorrection.μoρLabel(unknown, cxr) => uv(5813.0, 577.0), #
        NeXLMatrixCorrection.θLabel(unknown) => uv(θtoa, deg2rad(0.1)), #
        NeXLMatrixCorrection.dzLabel(unknown) => uv(0.0, 1.0e-6), # 10 nm
    )

    ABχ_res = cat(AB_res, χdata)

    maintain = MaintainInputs([NeXLMatrixCorrection.θLabel, NeXLMatrixCorrection.FLabel], ABχ_res)
    χFr = NeXLMatrixCorrection.StepχFr(unknown, cxr) | maintain
    χFr_res = χFr(ABχ_res)
    #using DataFrames, CSV
    #CSV.write("c:\\Users\\nritchie\\Desktop\\ABχ_res[1].csv",asa(DataFrame,ABχ_res))
    #χFr_model = χFr ∘ AB_model
    #χFr_mcres = mcpropagate(χFr, ABχ_res, 1000, parallel=false, rng=rgen)

    @test isapprox(value(χFr_res[NeXLMatrixCorrection.χLabel(unknown, cxr)]), 9.04e3, atol = 0.01e3)
    @test isapprox(value(χFr_res[NeXLMatrixCorrection.FrLabel(unknown, cxr)]), 2.85e-4, atol = 0.1e-4)
    @test isapprox(σ(χFr_res[NeXLMatrixCorrection.χLabel(unknown, cxr)]), 898, atol = 10)
    @test isapprox(σ(χFr_res[NeXLMatrixCorrection.FrLabel(unknown, cxr)]), 2.6e-5, atol = 0.1e-5)

    coatU = "10 nm C"
    coatingData = uvs( #
        NeXLMatrixCorrection.tcLabel(unknown) => uv(1.0e-6, 0.5e-6),
        NeXLMatrixCorrection.μoρLabel(coatU, cxr) => uv(11705.0, 2185.0),  # n"O K-L3" in n"C"
    )

    Frc_input = cat(χFr_res, coatingData)

    maintain = MaintainInputs([NeXLMatrixCorrection.θLabel, NeXLMatrixCorrection.FLabel], Frc_input)
    Frc = NeXLMatrixCorrection.StepFrc(unknown, coatU, cxr) | maintain
    Frcu_res = Frc(Frc_input)
    Frc_model = Frc
    Frc_mcres = mcpropagate(Frc, Frc_input, 1000, parallel = false, rng = rgen)

    @test isapprox(value(Frcu_res[NeXLMatrixCorrection.FrcLabel(unknown, coatU, cxr)]), 2.8e-4, atol = 0.1e-4)
    @test isapprox(σ(Frcu_res[NeXLMatrixCorrection.FrcLabel(unknown, coatU, cxr)]), 2.5e-5, atol = 0.1e-5)

    # Now calculate the standard

    input_uvs = uvs(
        (MassFractionLabel(standard, elm) => uv(std_mat[elm], 0.01 * std_mat[elm]) for elm in keys(std_mat))...,
        (AtomicWeightLabel(standard, elm) => uv(a(elm, std_mat), 0.001 * a(elm, std_mat)) for elm in keys(std_mat))...,
        (JzLabel(elm) => Ju(elm) for elm in keys(std_mat))...,
        NeXLMatrixCorrection.E0Label(standard) => uv(e0, 0.1e3),
        NeXLMatrixCorrection.mLabel(inner(cxr)) => uv(m, 0.01 * m),
    )

    maintain = MaintainInputs([NeXLMatrixCorrection.mLabel], input_uvs)
    mjz = StepMJZbarb(standard, [n"O", n"Si"]) | maintain

    mjz_res = mjz(input_uvs)
    mjz_mcres = mcpropagate(mjz, input_uvs, 10000, parallel = false, rng = rgen)

    # println("Analytical Result")
    # print(mjz_res)
    # println("MC Result")
    # print(mjz_mcres)

    dpt = StepDPT(standard, inner(cxr)) | AllInputs()

    dpt_res = dpt(mjz_res)
    dpt_mcres = mcpropagate(dpt ∘ mjz, input_uvs, 1000, parallel = false, rng = rgen)

    # println("Analytical Result")
    # print(dpt_res)
    # println("MC Result")
    # print(dpt_mcres)

    retain = MaintainInputs([NeXLMatrixCorrection.E0keVLabel, NeXLMatrixCorrection.ZbarbLabel], dpt_res)
    qla = StepQlaOoS(standard, inner(cxr)) | retain
    qla_res = qla(dpt_res)
    qla_model = qla ∘ dpt ∘ mjz
    qla_mcres = mcpropagate(qla_model, input_uvs, 10000, parallel = false, rng = rgen)

    #println("Analytical Result")
    #print(qla_res)
    #println("MC Result")
    #print(qla_mcres)

    maintain = MaintainInputs(
        [
            NeXLMatrixCorrection.E0keVLabel,
            NeXLMatrixCorrection.ZbarbLabel,
            NeXLMatrixCorrection.OoSLabel,
            NeXLMatrixCorrection.QlaLabel,
        ],
        qla_res,
    )
    rp = NeXLMatrixCorrection.StepRPhi0(standard, inner(cxr)) | maintain
    rp_res = rp(qla_res)
    rp_model = rp ∘ qla_model
    rp_mcres = mcpropagate(rp_model, input_uvs, 1000, parallel = false, rng = rgen)

    #println("Analytical Result")
    #print(rp_res)
    #println("MC Result")
    #print(rp_mcres)

    maintain = MaintainInputs(
        [NeXLMatrixCorrection.ZbarbLabel, NeXLMatrixCorrection.ϕ0Label, NeXLMatrixCorrection.E0keVLabel],
        rp_res,
    )
    frbar = NeXLMatrixCorrection.StepFRBar(standard, inner(cxr)) | maintain
    frbar_res = frbar(rp_res)
    frbar_model = frbar ∘ rp_model
    frbar_mcres = mcpropagate(frbar_model, input_uvs, 1000, parallel = false, rng = rgen)

    #println("Analytical Result")
    #print(frbar_res)
    #println("MC Result")
    #print(frbar_mcres)

    maintain = MaintainInputs(
        [NeXLMatrixCorrection.ϕ0Label, NeXLMatrixCorrection.RbarLabel, NeXLMatrixCorrection.FLabel],
        frbar_res,
    )
    pb = NeXLMatrixCorrection.StepPb(standard, inner(cxr)) | maintain
    pb_res = pb(frbar_res)
    pb_model = pb ∘ frbar_model
    pb_mcres = mcpropagate(pb_model, input_uvs, 1000, parallel = false, rng = rgen)

    maintain = MaintainInputs(
        [
            NeXLMatrixCorrection.ϕ0Label,
            NeXLMatrixCorrection.FLabel,
            NeXLMatrixCorrection.PLabel,
            NeXLMatrixCorrection.bLabel,
        ],
        pb_res,
    )
    aϵ = NeXLMatrixCorrection.Stepaϵ(standard, inner(cxr)) | maintain
    aϵ_res = aϵ(pb_res)
    aϵ_model = aϵ ∘ pb_model
    aϵ_mcres = mcpropagate(aϵ_model, input_uvs, 1000, parallel = false, rng = rgen)

    maintain = MaintainInputs(
        [
            NeXLMatrixCorrection.bLabel,
            NeXLMatrixCorrection.ϕ0Label,
            NeXLMatrixCorrection.ϵLabel,
            NeXLMatrixCorrection.FLabel,
        ],
        aϵ_res,
    )
    AB = NeXLMatrixCorrection.StepAB(standard, inner(cxr)) | maintain
    AB_res = AB(aϵ_res)
    AB_model = AB ∘ aϵ_model
    AB_mcres = mcpropagate(AB_model, input_uvs, 1000, parallel = false, rng = rgen)

    χdata = uvs(
        NeXLMatrixCorrection.μoρLabel(standard, cxr) => uv(4123.0, 668.0), #
        NeXLMatrixCorrection.θLabel(standard) => uv(θtoa, deg2rad(0.1)), #
        NeXLMatrixCorrection.dzLabel(standard) => uv(0.0, 1.0e-6), # 10 nm
    )

    ABχ_res = cat(AB_res, χdata)

    maintain = MaintainInputs([NeXLMatrixCorrection.θLabel, NeXLMatrixCorrection.FLabel], ABχ_res)
    χFr = NeXLMatrixCorrection.StepχFr(standard, cxr) | maintain
    χFr_res = χFr(ABχ_res)
    #@test isposdef(ABχ_res)
    #χFr_mcres = mcpropagate(χFr, ABχ_res, 1000, parallel = false, rng = rgen)

    coatS = "11 nm C"
    coatingData = uvs( #
        NeXLMatrixCorrection.tcLabel(standard) => uv(1.1e-6, 0.5e-6),
        NeXLMatrixCorrection.μoρLabel(coatS, cxr) => uv(mac(pure(n"C"), cxr), 2185.0),  # n"O K-L3" in n"C"
    )
    Frc_input = cat(χFr_res, coatingData)

    maintain = MaintainInputs([NeXLMatrixCorrection.θLabel, NeXLMatrixCorrection.FLabel], Frc_input)
    Frc = NeXLMatrixCorrection.StepFrc(standard, coatS, cxr) | maintain
    Frcs_res = Frc(Frc_input)
    Frc_model = Frc
    Frc_mcres = mcpropagate(Frc, Frc_input, 1000, parallel = false, rng = rgen)

    za_inputs = cat(Frcu_res, Frcs_res)

    za = NeXLMatrixCorrection.StepZA(unknown, standard, cxr, coatU, coatS)
    za_res = za(za_inputs)
    za_mcres = mcpropagate(za, za_inputs, 1000, parallel = false, rng = rgen)

    # Compare to DTSA-II and MC
    @test isapprox(value(za_res[NeXLMatrixCorrection.ZLabel(unknown, standard, inner(cxr))]), 1.1267, atol = 0.001)
    @test isapprox(
        value(za_res[NeXLMatrixCorrection.AbsLabel(unknown, standard, cxr, coatU, coatS)]),
        0.762,
        atol = 0.01,
    )
    @test isapprox(σ(za_res[NeXLMatrixCorrection.ZLabel(unknown, standard, inner(cxr))]), 0.024, atol = 0.001)
    @test isapprox(σ(za_res[NeXLMatrixCorrection.AbsLabel(unknown, standard, cxr, coatU, coatS)]), 0.14, atol = 0.01)

    zaf_za = ZAFc(zaf..., cxr, θtoa, θtoa)
    zaf_a = ZAFc(zaf..., cxr, θtoa, θtoa) / Z(zaf...)

    @test isapprox(value(za_res[NeXLMatrixCorrection.ZLabel(unknown, standard, inner(cxr))]), Z(zaf...), atol = 1e-5)
    @test_broken isapprox(
        value(za_res[NeXLMatrixCorrection.AbsLabel(unknown, standard, cxr, coatU, coatS)]),
        zaf_a,
        atol = 1e-5,
    )
    @test_broken isapprox(
        value(za_res[NeXLMatrixCorrection.ZALabel(unknown, standard, cxr, coatU, coatS)]),
        zaf_za,
        atol = 1e-5,
    )
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
using NPZ
using CodecZlib
using Random

"""
Represents an `Iterable` dataset with the following fields:

# Fields
- `dataset::Dict{Symbol, Any}`: representation of the dataset as a Dictionary with style as `style`.
- `epochs::Vector{Int}`: list of epochs loaded.
- `repo::String`: the repository from which the dataset is taken.
- `length::Int`: the length of the dataset.
- `batch_size::Int`: the size of the batches returned by `iterate`.
- `style::Tuple{Symbol}`: the style of the `Iterator` that is returned, check out: [`SARTS`](@ref), [`SART`](@ref) and [`SA`](@ref) for types supported out of the box.
- `rng<:AbstractRNG`.
- `meta::Dict`: the metadata provided along with the dataset.
- `is_shuffle::Bool`: determines if the batches returned by `iterate` are shuffled.
"""
struct AtariDataSet{T<:AbstractRNG} <:RLDataSet
    dataset::Dict{Symbol, Any}
    epochs::Vector{Int}
    repo::String
    length::Int
    batch_size::Int
    style::Tuple
    rng::T
    meta::Dict
    is_shuffle::Bool
end

const samples_per_epoch = Int(1e6)
const atari_frame_size = 84
const epochs_per_game = 50

"""
    dataset(dataset, index, epochs; <keyword arguments>)

Create a dataset enclosed in a [`AtariDataSet`](@ref) `Iterable` type. Contain other related metadata
for the `dataset` that is passed. The returned type is an infinite or a finite `Iterator`
respectively depending upon whether is_shuffle is `true` or `false`. For more information regarding
the dataset, refer to [google-research/batch_rl](https://github.com/google-research/batch_rl). Check out `atari_params()` for more info on arguments.

# Arguments

- `dataset::String`: name of the datset.
- `index::Int`: analogous to `v` and different values correspond to different `seed`s that are used for data collection. can be between `[1:5]`.
- `epochs::Vector{Int}`: list of epochs to load. included epochs should be between `[0:50]`.
- `style::NTuple=SARTS`: the style of the `Iterator` that is returned. can be [`SARTS`](@ref), [`SART`](@ref) or [`SA`](@ref).
- `repo::String="atari-replay-datasets"`: name of the repository of the dataset.
- `rng::AbstractRNG=StableRNG(123)`.
- `is_shuffle::Bool=true`: determines if the dataset is shuffled or not.
- `batch_size::Int=256` batch_size that is yielded by the iterator.

!!! warning

    The dataset takes up significant amount of space in RAM. Therefore it is advised to
    load even one epoch with 20GB of RAM. We are looking for ways to use lazy data loading here
    and any contributions are welcome.
"""
function dataset(
    game::String,
    index::Int,
    epochs::Vector{Int};
    style::NTuple=SARTS,
    repo::String="atari-replay-datasets",
    rng::AbstractRNG=MersenneTwister(123),
    is_shuffle::Bool=true,
    batch_size::Int=256
)

    try
        @datadep_str "$repo-$game-$index"
    catch e
        if isa(e, KeyError)
            throw("Invalid params, check out `atari_params()`")
        end
    end

    path = @datadep_str "$repo-$game-$index"

    @assert length(readdir(path)) == 1
    folder_name = readdir(path)[1]

    folder_path = "$path/$folder_name"
    files = readdir(folder_path)
    file_prefixes = collect(Set(map(x->join(split(x,"_")[1:2], "_"), files)))
    fields = map(collect(file_prefixes)) do x
        if split(x, "_")[1] == "\$store\$"
            x = split(x, "_")[2]
        else
            x = x
        end
    end

    s_epochs = Set(epochs)

    dataset = Dict()

    for (prefix, field) in zip(file_prefixes, fields)
        for epoch in s_epochs
            @assert epoch <= epochs_per_game
            data = open("$folder_path/$(prefix)_ckpt.$epoch.gz") do file
                stream = GzipDecompressorStream(file)
                NPZ.npzreadarray(stream)
            end

            if field == "observation"
                data = permutedims(data, [3, 2, 1]) # not sure about the orientation of the frame
            end

            if haskey(dataset, field)
                if field == "observation"
                    dataset[field] = cat(dataset[field], data, dims=3)
                else
                    dataset[field] = cat(dataset[field], data, dims=1)
                end
            else
                dataset[field] = data
            end
        end
    end

    num_epochs = length(s_epochs)

    atari_verify(dataset, num_epochs)

    N_samples = size(dataset["observation"])[3]

    final_dataset = Dict{Symbol, Any}()
    meta = Dict{String, Any}()

    for (key, d_key) in zip(["observation", "action", "reward", "terminal"], Symbol.(["state", "action", "reward", "terminal"]))
            final_dataset[d_key] = dataset[key]
    end

    for key in keys(dataset)
        if !(key in ["observation", "action", "reward", "terminal"])
            meta[key] = dataset[key]
        end
    end

    return AtariDataSet(final_dataset, epochs, repo, N_samples, batch_size, style, rng, meta, is_shuffle)

end

function iterate(ds::AtariDataSet, state = 0)
    rng = ds.rng
    batch_size = ds.batch_size
    length = ds.length
    is_shuffle = ds.is_shuffle
    style = ds.style

    if is_shuffle
        inds = rand(rng, 1:length-1, batch_size)
    else
        if (state+1) * batch_size <= length
            inds = state*batch_size+1:(state+1)*batch_size
        else
            return nothing
        end
        state += 1
    end

    batch = (state = view(ds.dataset[:state], :, :, inds),
    action = view(ds.dataset[:action], inds),
    reward = view(ds.dataset[:reward], inds),
    terminal = view(ds.dataset[:terminal], inds))

    if style == SARTS
        batch = merge(batch, (next_state = view(ds.dataset[:state], :, :, (1).+(inds)),))
    end

    return batch, state
end

take(ds::AtariDataSet, n::Integer) = take(ds.dataset, n)
length(ds::AtariDataSet) = ds.length
IteratorEltype(::Type{AtariDataSet}) = EltypeUnknown() # see if eltype can be known (not sure about carla and adroit)

function atari_verify(dataset::Dict, num_epochs::Int)
    @assert size(dataset["observation"]) == (atari_frame_size, atari_frame_size, num_epochs*samples_per_epoch)
    @assert size(dataset["action"]) == (num_epochs * samples_per_epoch,)
    @assert size(dataset["reward"]) == (num_epochs * samples_per_epoch,)
    @assert size(dataset["terminal"]) == (num_epochs * samples_per_epoch,)
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

include("bickleyjet_dg.jl")
include("bickleyjet_dg_reference.jl")

using BenchmarkTools

n1, n2 = 4096, 4096
Nq = 4

mesh = Meshes.RectilinearMesh(domain, n1, n2)
grid_topology = Topologies.Topology2D(mesh)
quad = Spaces.Quadratures.GLL{Nq}()
space = Spaces.SpectralElementSpace2D(grid_topology, quad)

X = coordinates(Val(Nq), n1, n2)
y0_ref = init_st_ref(X, Val(Nq), parameters)
dydt_ref = similar(y0_ref)
tendency_states = init_tendency_states(n1, n2, Val(Nq))
volume_ref!(
    dydt_ref,
    y0_ref,
    (n1, n2, parameters, Val(Nq), tendency_states),
    0.0,
)

eltm = @belapsed volume_ref!(
    $dydt_ref,
    $y0_ref,
    ($n1, $n2, $parameters, $(Val(Nq)), $tendency_states),
    0.0,
)

println(eltm)


text:
#!/usr/bin/env julia

using TrainRuns

train = Train("test/data/trains/freight.yaml")
path  = Path("test/data/paths/const.yaml")

runtime = trainrun(train, path)

println("The train needs $runtime seconds for the running path.")


text:
"""
    DockerContainer()
"""
Base.@kwdef struct DockerContainer
    label::String = Random.randstring(10)
end

"""
    DockerConfig(; kwargs...)

## Required Keyword Arguments:
- `base_image::String`

## Optional Keyword Arguments:
- `verbose::Bool = false`
- `env::Union{Dict{String, String}, Nothing} = nothing`
- `platform::Symbol = :linux`
- `read_only_maps::Union{Dict{String, String}, Nothing} = nothing`
- `read_write_maps::Union{Dict{String, String}, Nothing} = nothing`
- `stdin::IO = Base.devnull`
- `stdout::IO = Base.stdout`
- `stderr::IO = Base.stderr`
- `docker_build_stdout::Union{IO, Nothing} = nothing`
- `docker_build_stderr::Union{IO, Nothing} = nothing`
- `add_capabilities::Union{Vector{String}, Nothing} = nothing`
"""
Base.@kwdef struct DockerConfig
    # Required
    base_image::String

    # Optional
    verbose::Bool = false
    env::Union{Dict{String, String}, Nothing} = nothing
    platform::Symbol = :linux
    read_only_maps::Union{Dict{String, String}, Nothing} = nothing
    read_write_maps::Union{Dict{String, String}, Nothing} = nothing
    stdin::IO = Base.devnull
    stdout::IO = Base.stdout
    stderr::IO = Base.stderr
    docker_build_stdout::Union{IO, Nothing} = nothing
    docker_build_stderr::Union{IO, Nothing} = nothing
    add_capabilities::Union{Vector{String}, Nothing} = nothing
end


text:
using Pkg
Pkg.activate("JuliaTutor")
using JuliaTutor

JuliaTutor.greet()

"""this is some sort of prompt or something.
maybe some long instruction set!
maybe someone droning on and on about julia?""" |> inform

request_read_evaluate( "Define a variable `a` to be the integer 1234",
                    "a=1234")


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

# k : r-degree
# l : θ-degree
# m : z-degree

import ACE: alloc_B, alloc_dB
import Base: ==
import JuLIP: evaluate!, evaluate_d!,
              alloc_temp, alloc_temp_d,
              read_dict, write_dict

import JuLIP.MLIPs: IPBasis

"""
`struct EnvPairBasis`

The basis type for environment-dependent pair potentials. This models functions
of the form
```
   V(R; {Rⱼ})
```
where `R` is the bond-direction and `{Rⱼ}` the environment. V has the following
symmetries:
 - permutations of the environment
 - rotation of the environment about the R axis
 - rotation or reflection of the entire configuration

### Construction

Use `ACE.Bonds.envpairbasis` to construct a basis.

### Evaluation

To evaluate a basis with bond-direction `R` and environment
`Renv = [R1, R2,...]` use
```
evaluate(basis, R, Renv)            # allocating
evaluate!(B, tmp, basis, R, Renv)   # non-allocating
```

### After fitting parameters, to evaluate the resulting potential

<!!!TODO!!!>
"""
struct EnvPairBasis{T0, TR, TZ, TT, TI} <: IPBasis
   P0::T0                # basis for the bond-length coordinate / m0 = k0
   Pr::TR                # specifies the radial basis  / n = kr
   Pθ::TT                # the angular basis           / l = kθ
   Pz::TZ                # specifies the z-basis       / m = kz
   aalist::BondAAList{TI}    # datastructure specifying the basis
end

Base.length(basis::EnvPairBasis) = length(basis.aalist)

fltype(basis::EnvPairBasis) = Complex{Float64}   # fltype(basis.Pr)

alloc_B(basis::EnvPairBasis, args...) = zeros(fltype(basis), length(basis))

alloc_temp(basis::EnvPairBasis, args...) =
   ( P0 = alloc_B(basis.P0),
     tmp_P0 = alloc_temp(basis.P0),
     Pr = alloc_B(basis.Pr),
     tmp_Pr = alloc_temp(basis.Pr),
     Pθ = alloc_B(basis.Pθ),
     tmp_Pθ = alloc_temp(basis.Pθ),
     Pz = alloc_B(basis.Pz),
     tmp_Pz = alloc_temp(basis.Pz),
     A = alloc_A(basis.aalist.alist)
    )

function precompute_A!(A, tmp, basis::EnvPairBasis, R0, Rs)
   alist = basis.aalist.alist
   fill!(A, 0)
   # construct the coordinate system, and convert
   cylcoords = CylindricalCoordinateSystem(R0, R0/2)
   # loop through the environment to assemble the As
   for R in Rs
      rθz = cylcoords(R)
      evaluate!(tmp.Pr, tmp.tmp_Pr, basis.Pr, rθz.r)
      evaluate!(tmp.Pz, tmp.tmp_Pz, basis.Pz, rθz.z)
      evaluate!(tmp.Pθ, tmp.tmp_Pθ, basis.Pθ, rθz)
      for i = 1:length(alist)
         krθz = alist[i]
         A[i] += tmp.Pr[krθz.kr+1] *
                 tmp.Pθ[cyl_l2i(krθz.kθ, basis.Pθ)] *
                 tmp.Pz[krθz.kz+1]
      end
      # iz = z2i(ship, Z)
      # for i = alist.firstz[iz]:(alist.firstz[iz+1]-1)
      #    zklm = alist[i]
      #    A[i] += tmp.Pr[zklm.k+1] * tmp.Pθ[cyl_l2i(zklm.l)] * tmp.Pz[zklm.m+1]
      # end
   end
   return A
end

# R0   : typically SVector{T}
# Renv : typically Vector{SVector{T}} or a view into Vector{SVector{T}}
function evaluate!(B::AbstractVector{Complex{T}},
                   tmp, basis::EnvPairBasis, R0, Renv)  where {T}
   aalist = basis.aalist
   # construct the basis for the r1-variable
   r0 = norm(R0)
   P0 = evaluate!(tmp.P0, tmp.tmp_P0, basis.P0, r0)

   # evaluate the A-basis, i.e. the density projections of the environment
   A = precompute_A!(tmp.A, tmp, basis, R0, Renv)

   # loop over all basis functions
   for i = 1:length(aalist)
      B[i] = P0[aalist.i2k0[i] + 1]
      for α = 1:aalist.len[i]
         B[i] *= A[aalist.i2Aidx[i, α]]  # TODO: reverse i, α
      end
   end

   return B
end

# ------------------------------------------------------------

using JuLIP: @pot, AbstractCalculator
import JuLIP.MLIPs: combine

"""
A rudimentary implementation of an environment-dependent pair potential.
"""
struct EnvPairPot{T, T0, TR, TZ, TT, TI} <: AbstractCalculator
   basis::EnvPairBasis{T0, TR, TZ, TT, TI}
   c::Vector{T}
end

@pot EnvPairPot

combine(b::EnvPairBasis, c::AbstractVector{<: Number}) = EnvPairPot(b, c)

alloc_temp(V::EnvPairPot, args...) = alloc_temp(V.basis, args...)

# R0   : typically SVector{T}
# Renv : typically Vector{SVector{T}} or a view into Vector{SVector{T}}
function evaluate!(tmp, V::EnvPairPot, R0, Renv)
   val = zero(ComplexF64)

   aalist = V.basis.aalist
   # construct the basis for the r1-variable
   r0 = norm(R0)
   P0 = evaluate!(tmp.P0, tmp.tmp_P0, V.basis.P0, r0)
   # evaluate the A-basis, i.e. the density projections of the environment
   A = precompute_A!(tmp.A, tmp, V.basis, R0, Renv)

   # loop over all basis functions
   for i = 1:length(aalist)
      b = P0[aalist.i2k0[i] + 1]
      for α = 1:aalist.len[i]
         b *= A[aalist.i2Aidx[i, α]]
      end
      val += V.c[i] * b
   end

   return real(val)
end

# --------------(de-)serialisation----------------------------------------

import JuLIP: write_dict, read_dict

write_dict(b::EnvPairBasis) =
      Dict( "__id__" => "ACE_EnvPairBasis",
            "P0"     => write_dict(b.P0),
            "Pr"     => write_dict(b.Pr),
            "Ptheta" => write_dict(b.Pθ),
            "Pz"     => write_dict(b.Pz),
            "aalist" => write_dict(b.aalist) )

read_dict(::Val{:ACE_EnvPairBasis}, D::Dict) =
      EnvPairBasis(read_dict(D["P0"]),
                   read_dict(D["Pr"]),
                   read_dict(D["Ptheta"]),
                   read_dict(D["Pz"]),
                   read_dict(D["aalist"]))

write_dict(V::EnvPairPot) =
      Dict( "__id__" => "ACE_EnvPairPot",
            "basis" => write_dict(V.basis),
            "cr" => real.(V.c),
            "ci" => imag.(V.c) )

read_dict(::Val{:ACE_EnvPairPot}, D::Dict) =
      EnvPairPot( read_dict(D["basis"]),
                  D["cr"] + im * D["ci"] )

==(B1::EnvPairBasis, B2::EnvPairBasis) =
         all(getfield(B1, x) == getfield(B2, x) for x in fieldnames(EnvPairBasis))

==(B1::EnvPairPot, B2::EnvPairPot) =
         all(getfield(B1, x) == getfield(B2, x) for x in fieldnames(EnvPairPot))


text:
"""
`@=>(expr...)`

Create a selector based on expressions `expr`. Symbols are used to select columns and infer an appropriate anonymous function.
In this context, `_` refers to the whole row. To use actual symbols, escape them with `^`, as in `^(:a)`.
Use `cols(c)` to refer to field `c` where `c` is a variable that evaluates to a symbol. `c` must be available in
the scope where the macro is called.

## Examples

```jldoctest pair
julia> t = table((a = [1,2,3], b = [4,5,6]));

julia> select(t, @=>(:a, :a + :b))
Table with 3 rows, 2 columns:
a  a + b
────────
1  5
2  7
3  9
```
"""
macro =>(args...)
    esc(pair_helper(args...))
end

function pair_helper(expr; name = nothing)
    isquotenode(expr) && (name == nothing) && return expr
    ispair(expr) && (name == nothing) && return pair_helper(expr.args[3]; name = expr.args[2]) #Expr(:call, :(=>), expr.args[2], pair_helper(expr.args[3]))
    istuple(expr)  && return Expr(:tuple, map(pair_helper, expr.args)...)
    anon_func, syms = extract_anonymous_function(expr, replace_field)
    fields = (:(_) in syms || isempty(syms)) ? Expr(:call, :(JuliaDBMeta.All)) : Expr(:call, :(JuliaDBMeta.distinct_tuple), syms...)
    name == nothing && (name = Expr(:quote, Symbol(filter(t -> t != ':', string(expr)))))
    :($name => $fields => $anon_func)
end

pair_helper(expr...) = pair_helper(Expr(:tuple, expr...))

"""
`@select(d, x)`

Short-hand for `select(d, @=>(x))`

## Examples

```jldoctest select
julia> t = table((a = [1,2,3], b = [4,5,6]));

julia> @select(t, (:a, :a + :b))
Table with 3 rows, 2 columns:
a  a + b
────────
1  5
2  7
3  9
```
"""
macro select(args...)
    esc(select_helper(args...))
end

function select_helper(x)
    i = gensym()
    Expr(:(->), i, select_helper(i, x))
end

function select_helper(d, x)
    :(JuliaDBMeta.select($d, $(pair_helper(x))))
end


text:
module WWUMeshes

export read_obj, write_obj
export gen_mesh, est_normals
export cube_mesh, cylinder_mesh, sphere_mesh, estimate_normals
export OBJTriangle, OBJMesh

using FileIO
using LinearAlgebra

push!(LOAD_PATH, pwd())

include("GfxBase.jl")
using .GfxBase

""" OBJTriangle
A struct that represents a single triangle in a mesh. """
mutable struct OBJTriangle
    positions::Array{Int, 1} # vertex position indices
    uvs::Array{Int, 1} # vertex texture coordinate indices
    normals::Array{Int, 1} # normal vector indices
end

""" OBJMesh
A struct that represents an indexed triangle mesh for reading from
or writing to OBJ format. """
mutable struct OBJMesh
    positions::Array{Vec3, 1} # all vertex positions
    uvs::Array{Vec2, 1} # all texture coordinates
    normals::Array{Vec3, 1} # all vertex normals
    triangles::Array{OBJTriangle, 1} # the OBJTriangles belonging to the mesh
end

""" read_obj(obj_filename)
Read a mesh in OBJ format from file obj_filename."""
function read_obj(obj_filename)
    m = OBJMesh([], [], [], []) # create a mesh
    open(obj_filename) do f
        for (line_number, line) in enumerate(eachline(f))
            if line == "" || line[1] == "#"
                continue # skip comments
            end
            # Read the line and add its contents to the correct field of m:
            tokens = split(strip(line))
            if tokens[1] == "v" # vertex
                push!(m.positions, Vec3([parse(Float32, x) for x in tokens[2:end]]...))
            elseif tokens[1] == "vt" # vertex texture
                push!(m.uvs, Vec2([parse(Float32, x) for x in tokens[2:3]]...))
            elseif tokens[1] == "vn" # vertex normal
                push!(m.normals, Vec3([parse(Float32, x) for x in tokens[2:end]]...))
            elseif tokens[1] == "f"
                if length(tokens) == 4
                    # create a OBJTriangle face:
                    points = []
                    uvs = []
                    normals = []
                    # handle faces with no texture and/or normals
                    for corner in tokens[2:4]
                        indices = split(corner, '/')
                        if length(indices) == 3 # all 3 present, third is a normal
                            push!(normals, parse(Int, indices[3]))
                        end
                        if length(indices) >= 2 && indices[2] != ""
                            # if there are 2 or more and the second isn't blank, it's a texture
                            push!(uvs, parse(Int, indices[2]))
                        end
                        if length(indices) >= 1 # first value is the position
                            push!(points, parse(Int, indices[1]))
                        else # unless it has none, in which case it's not valid
                            error("in line $line_number: face vertex $corner could not be parsed")
                        end
                    end
                    # create the triangle and add it to the triangles array
                    push!(m.triangles, OBJTriangle(points, uvs, normals))
                else
                    # create a OBJTriangle face:
                    points = []
                    uvs = []
                    normals = []
                    # handle faces with no texture and/or normals
                    for corner in tokens[2:4]
                        indices = split(corner, '/')
                        if length(indices) == 3 # all 3 present, third is a normal
                            push!(normals, parse(Int, indices[3]))
                        end
                        if length(indices) >= 2 && indices[2] != ""
                            # if there are 2 or more and the second isn't blank, it's a texture
                            push!(uvs, parse(Int, indices[2]))
                        end
                        if length(indices) >= 1 # first value is the position
                            push!(points, parse(Int, indices[1]))
                        else # unless it has none, in which case it's not valid
                            error("in line $line_number: face vertex $corner could not be parsed")
                        end
                    end
                    # create the triangle and add it to the triangles array
                    push!(m.triangles, OBJTriangle(points, uvs, normals))

                    points = []
                    uvs = []
                    normals = []
                    # handle faces with no texture and/or normals
                    for corner in [tokens[4], tokens[5], tokens[2]]
                        indices = split(corner, '/')
                        if length(indices) == 3 # all 3 present, third is a normal
                            push!(normals, parse(Int, indices[3]))
                        end
                        if length(indices) >= 2 && indices[2] != ""
                            # if there are 2 or more and the second isn't blank, it's a texture
                            push!(uvs, parse(Int, indices[2]))
                        end
                        if length(indices) >= 1 # first value is the position
                            push!(points, parse(Int, indices[1]))
                        else # unless it has none, in which case it's not valid
                            error("in line $line_number: face vertex $corner could not be parsed")
                        end
                    end
                    # create the triangle and add it to the triangles array
                    push!(m.triangles, OBJTriangle(points, uvs, normals))
                end
            end
        end
    end
    return m
end

""" write_obj(obj_filename)
Write the given mesh in OBJ format to file obj_filename."""
function write_obj(obj_filename, mesh::OBJMesh)
    open(obj_filename, "w") do f
        # write all positions:
        for v in mesh.positions
            write(f, "v $(v[1]) $(v[2]) $(v[3])\n")
        end

        # write all texture coords:
        for v in mesh.uvs
            write(f, "vt $(v[1]) $(v[2])\n")
        end
        # write all normals:
        for v in mesh.normals
            write(f, "vn $(v[1]) $(v[2]) $(v[3])\n")
        end

        # write all triangles:
        for tri in mesh.triangles
            write(f, "f $(tri_vertex_str(tri))\n")
        end

    end

end

""" tri_vertex_str(triangle)
Return a string with the indices of applicable positions, texture coordinates,
and normals for a given triangle according to the OBJ specification.
In particular, if p, u, and n are position, vertex and normal, each corner
of the triangle is represented as one of the following:
    p       (position only)
    p/u     (position and texture)
    p//n    (position and normal)
    p/u/n   (position, texture, and normal) """
function tri_vertex_str(triangle::OBJTriangle)
    # determine whether textures and normals are present:
    write_uv = length(triangle.uvs) == length(triangle.positions)
    write_normals = length(triangle.normals) == length(triangle.positions)
    corners = []
    for i = 1:3
        output = "$(triangle.positions[i])"
        if write_uv && !write_normals
            output = output * "/$(triangle.uvs[i])" # * does concatenation(!)
        elseif !write_uv && write_normals
            output = output * "//$(triangle.normals[i])"
        elseif write_uv && write_normals
            output = output * "/$(triangle.uvs[i])/$(triangle.normals[i])"
        end
        push!(corners, output)
    end
    join(corners, " ")
end

""" gen_mesh(outfile, geom, divisionsU, divisionsV)
Generate a mesh and save the result in a file with name outfile.
geom may be "cube", "cylinder", or "sphere".
Cylinder requires divisionsU; sphere requires divisionsU and divisionsV. """
function gen_mesh(outfile, geom, divisionsU=0, divisionsV=0)
    if geom == "cube"
        mesh = cube_mesh()
    elseif geom == "cylinder"
        mesh = cylinder_mesh(divisionsU)
    elseif geom == "sphere"
        mesh = sphere_mesh(divisionsU, divisionsV)
    end
    write_obj(outfile, mesh)
end

""" est_normals(outfile, infile)
Estimate normals of the mesh stored in infile, saving the result in outfile."""
function est_normals(outfile, infile)
    input_mesh = read_obj(infile)
    mesh = estimate_normals(input_mesh)
    write_obj(outfile, mesh)
end

""" cube_mesh()
Return a new OBJMesh representing a 2x2x2 cube centered at the origin and
axis-aligned. """
function cube_mesh()
    positions = []
    uvs = []
    normals = []
    triangles = []
    # key to comments:
    # L/R = x = right/left
    # B/T = y = top/bottom
    # C/F = z = close/far
    push!(positions, Vec3( 1, -1, -1)) # 1 RBC
    push!(positions, Vec3( 1, -1,  1)) # 2 RBF
    push!(positions, Vec3(-1, -1,  1)) # 3 LBF
    push!(positions, Vec3(-1, -1, -1)) # 4 LBC
    push!(positions, Vec3( 1,  1, -1)) # 5 RTC
    push!(positions, Vec3( 1,  1,  1)) # 6 RTF
    push!(positions, Vec3(-1,  1,  1)) # 7 LTF
    push!(positions, Vec3(-1,  1, -1)) # 8 LTC

    # texture coordinates:
    push!(uvs, Vec2(1, 1)) # TR
    push!(uvs, Vec2(0, 1)) # TL
    push!(uvs, Vec2(0, 0)) # BL
    push!(uvs, Vec2(1, 0)) # BR

    # normals:
    push!(normals, Vec3( 1, 0, 0)) # R
    push!(normals, Vec3(-1, 0, 0)) # L
    push!(normals, Vec3( 0, 1, 0)) # U
    push!(normals, Vec3( 0,-1, 0)) # D
    push!(normals, Vec3( 0, 0, 1)) # C
    push!(normals, Vec3( 0, 0,-1)) # F

    # 8 faces, 2 triangles each
    push!(triangles, OBJTriangle([1,2,3], [1,2,3], [4,4,4])) # bottom face 1
    push!(triangles, OBJTriangle([1,3,4], [1,3,4], [4,4,4])) # bottom face 2
    push!(triangles, OBJTriangle([1,5,6], [4,1,2], [1,1,1])) # right face 1
    push!(triangles, OBJTriangle([1,6,2], [4,2,3], [1,1,1])) # right face 2
    push!(triangles, OBJTriangle([2,6,7], [4,1,2], [5,5,5])) # far face 1
    push!(triangles, OBJTriangle([2,7,3], [4,2,3], [5,5,5])) # far face 2
    push!(triangles, OBJTriangle([3,7,8], [2,3,4], [2,2,2])) # left face 1
    push!(triangles, OBJTriangle([3,8,4], [2,4,1], [2,2,2])) # left face 2
    push!(triangles, OBJTriangle([4,8,5], [2,3,4], [6,6,6])) # far face 1
    push!(triangles, OBJTriangle([4,5,1], [2,4,1], [6,6,6])) # far face 2
    push!(triangles, OBJTriangle([5,8,7], [1,2,3], [3,3,3])) # top face 1
    push!(triangles, OBJTriangle([5,7,6], [1,3,4], [3,3,3])) # top face 2

    # julia automatically returns the last value in the function:
    OBJMesh(positions, uvs, normals, triangles)

end

""" cylinder_mesh(n)
Return a new OBJMesh object approximation of a cylinder with radius 1 and
height 2, centered at the origin. The logitudinal axis is aligned with y, and
it is tesselated with n divisions arranged radially around the outer surface.
The ends of the cylinder are disc-shaped caps parallel to the xz plane. See the
assignment writeup for a diagram and details.
"""
function cylinder_mesh(divisionsU)
    positions = []
    uvs = []
    normals = []
    triangles = []

    # bottom center vertex
    push!(positions, Vec3(0, -1, 0))
    push!(uvs, Vec2(0.25, 0.75))
    push!(normals, Vec3(0, -1, 0))

    # top center vertex
    push!(positions, Vec3(0, 1, 0))
    push!(uvs, Vec2(0.75, 0.75))
    push!(normals, Vec3(0, 1, 0))

    # partially unwind loop where i = 0
    push!(positions, Vec3(1, 1, 0))
    push!(positions, Vec3(1, -1, 0))

    push!(uvs, Vec2(0.5, 0.75))
    push!(uvs, Vec2(1, 0.75))
    push!(uvs, Vec2(1, 0.5))
    push!(uvs, Vec2(1, 0))

    push!(normals, Vec3(1, 0, 0))

    delta = (2.0 * pi) / divisionsU
    for i in 1:divisionsU
        theta = delta * i
        x = cos(theta)
        z = sin(theta)

        push!(positions, Vec3(x, 1, z))
        push!(positions, Vec3(x, -1, z))

        push!(uvs, Vec2(0.25 + x / 4, 0.75 + z / 4))
        push!(uvs, Vec2(0.75 + x / 4, 0.75 + z / 4))
        push!(uvs, Vec2(1 - theta / (2.0 * pi), 0.5))
        push!(uvs, Vec2(1 - theta / (2.0 * pi), 0))

        push!(normals, Vec3(x, 0, z))

        # sides
        push!(triangles, OBJTriangle(
            [i * 2 + 2, i * 2 + 1, i * 2 + 3],
            [i * 4 + 2, i * 4 + 1, i * 4 + 5],
            [i + 2, i + 2, i + 3]))

        push!(triangles, OBJTriangle(
            [i * 2 + 3, i * 2 + 4, i * 2 + 2],
            [i * 4 + 5, i * 4 + 6, i * 4 + 2],
            [i + 3, i + 3, i + 2]))

        # caps
        push!(triangles, OBJTriangle(
            [1, i * 2 + 2, i * 2 + 4],
            [1, i * 4 - 1, i * 4 + 3],
            [1, 1, 1]))

        push!(triangles, OBJTriangle(
            [i * 2 + 1, 2, i * 2 + 3],
            [i * 4, 2, i * 4 + 4],
            [2, 2, 2]))
    end

    OBJMesh(positions, uvs, normals, triangles)
end

""" sphere_mesh(n, m)
Create a Latitude-Longitude-tesselated approximation of a sphere with radius 1
centered at the origin. There are n divisions around the equator and m
divisions from pole to pole along each line of longitude. The North pole is at
(0,1,0), the South pole at (0,-1,0), and points on the Greenwich meridian are
in the x = 0 plane with z > 0. The u texture coordinate depends on longitude,
with u=0 at 180 degrees West and u=1 at 180 degrees East. The v texture
coordinate varies with latitude with v=0 at the South pole and v=1 at the North
pole. Normals should be normal to the ideal sphere surface. See the assignment
for a diagram and further details. """
function sphere_mesh(n, m)
    positions = []
    uvs = []
    normals = []
    triangles = []

    deltaM = pi / m
    deltaN = (2.0 * pi) / n
    deltaV = 1 / m
    deltaU = 1 / n
    i = 0

    # partially unwind loop where j = 0
    for k in 1:n
        azimuth = deltaN * k

        push!(positions, Vec3(0, 1, 0))
        push!(uvs, Vec2(deltaU * k, 1))
        push!(normals, Vec3(0, 1, 0))
    end

    for j in 1:m + 1
        inclination = deltaM * j

        y = cos(inclination)

        # partially unwind loop where k = 0
        x = sin(inclination)

        push!(positions, Vec3(x, y, 0))
        push!(uvs, Vec2(0, 1 - deltaV * j))
        push!(normals, normalize(Vec3(x, y, 0)))

        for k in 1:n
            azimuth = deltaN * k

            x = sin(inclination) * cos(azimuth)
            z = sin(inclination) * sin(azimuth)

            i += 1
            push!(positions, Vec3(x, y, z))
            push!(uvs, Vec2(deltaU * k, 1 - deltaV * j))
            push!(normals, normalize(Vec3(x, y, z)))

            push!(triangles, OBJTriangle(
                [i, i + 1, i + n + 1],
                [i, i + 1, i + n + 1],
                [i, i + 1, i + n + 1]))

            push!(triangles, OBJTriangle(
                [i + n, i, i + n + 1],
                [i + n, i, i + n + 1],
                [i + n, i, i + n + 1]))
        end
    end

    OBJMesh(positions, uvs, normals, triangles)
end

"""
    estimate_normals(mesh::OBJMesh)
Estimates normals for the given mesh. Overwrites any existing normals and returns a new OBJMesh object.
"""
function estimate_normals(mesh::OBJMesh)
    mesh.normals = zeros(Vec3, length(mesh.positions))

    for t in 1:length(mesh.triangles)
        triangle = mesh.triangles[t]
        triangle.normals = zeros(3)

        v = mesh.positions[triangle.positions[2]] - mesh.positions[triangle.positions[1]]
        w = mesh.positions[triangle.positions[3]] - mesh.positions[triangle.positions[1]]

        normal = Vec3(
            (v[2] * w[3]) - (v[3] * w[2]),
            (v[3] * w[1]) - (v[1] * w[3]),
            (v[1] * w[2]) - (v[2] * w[1]))

        for v in 1:3
            mesh.normals[triangle.positions[v]] += normal
            triangle.normals[v] = triangle.positions[v]
        end
    end

    for i in 1:length(mesh.normals)
        mesh.normals[i] = normalize(mesh.normals[i])
    end

    mesh
end

end # module WWUMeshes

max_stars_repo_name: aedans/raytracer

max_issues_repo_name: aedans/raytracer

max_forks_repo_name: aedans/raytracer