text:
export TracerProvider

using Base.Threads
using Dates: time

struct Tracer{
    S<:AbstractSampler,
    SP<:AbstractSpanProcessor,
    G<:AbstractIdGenerator
} <: API.AbstractTracer
    sampler::S
    resource::Resource
    span_processor::SP
    id_generator::G
    instrumentation::InstrumentationInfo
end

function API.create_span(
    name::String,
    t::Tracer;
    parent_context=current_context(),
    kind=SPAN_KIND_INTERNAL,
    attributes=Attributes(;is_mutable=true),
    links=[],
    start_time=time()
)
    parent_span_ctx = parent_context |> current_span |> span_context
    if parent_span_ctx === INVALID_SPAN_CONTEXT
        trace_id = generate_trace_id(t.id_generator)
    else
        trace_id = parent_span_ctx.trace_id
    end

    sampling_result = should_sample(
        t.sampler,
        parent_context,
        trace_id,
        name,
        kind,
        attributes,
        links,
        parent_span_ctx.trace_state
    )

    span_ctx = SpanContext(
        trace_id = trace_id,
        span_id = generate_span_id(t.id_generator),
        is_remote=false,
        trace_flag = TraceFlag(sampled=is_sampled(sampling_result)),
        trace_state=sampling_result.trace_state
    )

    if is_recording(sampling_result)
        Span(
            ;name=name,
            span_context=span_ctx,
            parent_context=parent_context,
            kind=kind,
            attributes=attributes,
            links=links,
            start_time=start_time
        )
    else
        NonRecordingSpan(span_ctx)
    end
end

#####
# TracerProvider
#####

Base.@kwdef struct TracerProvider{
    S<:AbstractSampler,
    IDG<:AbstractIdGenerator
} <: API.AbstractTracerProvider
    sampler::S=DEFAULT_ON
    resource::Resource=Resource()
    span_processor::MultiSpanProcessor=MultiSpanProcessor()
    id_generator::IDG=RandomIdGenerator()
end

shut_down!(p::TracerProvider) = shut_down!(p.span_processor)

force_flush!(p::TracerProvider, args...) = force_flush!(p.span_processor, args...)

function API.get_tracer(p::TracerProvider, instrumentation_name, instrumentation_version=nothing)
    Tracer(
        p.sampler,
        p.resource,
        p.span_processor,
        p.id_generator,
        isnothing(instrumentation_version) ? InstrumentationInfo(instrumentation_name) : InstrumentationInfo(instrumentation_name, instrumentation_version)
    )
end


text:
struct SubDiskArray{T,N} <: AbstractDiskArray{T,N}
  v::SubArray{T,N}
end
replace_colon(s,::Colon) = Base.OneTo(s)
replace_colon(s,r) = r
function Base.view(a::AbstractDiskArray,i...)
  i2 = replace_colon.(size(a),i)
  SubDiskArray(SubArray(a,i2))
end
Base.view(a::AbstractDiskArray, i::CartesianIndices) = view(a,i.indices...)
function Base.view(a::SubDiskArray,i...)
  SubDiskArray(view(a.v,i...))
end
Base.view(a::SubDiskArray, i::CartesianIndices) = view(a,i.indices...)
function readblock!(a::SubDiskArray,aout,i::OrdinalRange...)
  pinds = parentindices(view(a.v,i...))
  inds,resh = interpret_indices_disk(parent(a.v),pinds)
  readblock!(parent(a.v),reshape(aout,map(length,inds)...),inds...)
end
function writeblock!(a::SubDiskArray,v,i::OrdinalRange...)
  pinds = parentindices(view(a.v,i...))
  inds,resh = interpret_indices_disk(parent(a.v),pinds)
  writeblock!(parent(a.v),reshape(v,map(length,inds)...),inds...)
end
Base.size(a::SubDiskArray) = size(a.v)
eachchunk(a::SubDiskArray) = eachchunk_view(haschunks(a.v.parent),a.v)
Base.parent(a::SubDiskArray) = a.v.parent
function eachchunk_view(::Chunked, vv)
  pinds = parentindices(vv)
  iomit = findints(pinds)
  chunksparent = eachchunk(parent(vv))
  #TODO this currently only works if the parent chunks are GridChunks, should be generalized in the future
  parentoffset, parentsize = chunksparent.offset, chunksparent.chunksize
  offsets = ([mod1(first(pinds[i]),parentsize[i])+parentoffset[i]-1 for i in 1:length(pinds) if !in(i,iomit)]...,)
  newcs = ([parentsize[i] for i in 1:length(pinds) if !in(i,iomit)]...,)
  GridChunks(vv,newcs,offset = offsets)
end
eachchunk_view(::Unchunked, a) = GridChunks(a,estimate_chunksize(a))
haschunks(a::SubDiskArray) = haschunks(parent(a.v))

ext: jl


text:
"""
Allometry. Scaling rules that relate size to mass
"""
abstract type AbstractAllometry end

update_height!(o, u) = update_height!(allometry_pars(o), o, u)
update_height!(f::Nothing, o, u) = nothing
update_height!(f::AbstractAllometry, o, u) =
    set_height!(o, allometric_height(f, mass(o, u)))

"""
    SqrtAllometry(β0, β)

Height is given by the square root of mass above the initial mass `β0`,
multiplied by the scalar `β1`.

$(FIELDDOCTABLE)
"""
@flattenable @columns struct SqrtAllometry{B0,B1} <: AbstractAllometry
    # Field       | Flatn | Default | Unit | Bounds         | Log  | Description
    β0::B0        | false | 1e-4*24 | g    | (1e-6, 10.00)  | true | "Intercept. Mass at Om"
    β1::B1        | true  | 0.1     | m    | (1e-2, 1.0)    | true | "Scalar for conversion to meters"
end

allometric_height(f::SqrtAllometry, mass) =
    sqrt((mass - f.β0) / unit(f.β0)) * f.β1

"""
    Allometry(β0, β, α)

Simple allometric relationship between mass and height

$(FIELDDOCTABLE)
"""
@flattenable @columns struct Allometry{B0,B1,A} <: AbstractAllometry
    β0::B0        | false | 1e-4*24 | g    | (1e-6, 10.00)  | true | "Intercept. Mass at Om"
    β1::B1        | true  | 0.1     | m    | (1e-5, 10.00)  | true | "Scalar for conversion to meters"
    α::A          | true  | 0.1     | _    | (1e-3, 100.00) | true | "Exponent relating mass to vertical dimension"
end

allometric_height(f::Allometry, mass) =
    f.β1 * ((max(zero(mass), mass - f.β0)) / unit(f.β0))^f.α

"""
    FixedAllometry(height)

Height is fixed at `height`, independent of mass.

$(FIELDDOCTABLE)
"""
@flattenable @columns struct FixedAllometry{M} <: AbstractAllometry
    height::M     | true  | 1.0     | m    | (1e-3, 100.00) | true | "Fixed height or depth"
end

allometric_height(f::FixedAllometry, mass) = f.height

mass(o, u) = u[:V] * w_V(o)


text:
## hashing a single value ##

hash(x::Any) = hash(x, zero(Uint))
hash(w::WeakRef, h::Uint) = hash(w.value, h)

## hashing general objects ##

hash(x::ANY, h::Uint) = hash(object_id(x), h)

## core data hashing functions ##

function hash_64_64(n::Uint64)
    local a::Uint64 = n
    a = ~a + a << 21
    a =  a $ a >> 24
    a =  a + a << 3 + a << 8
    a =  a $ a >> 14
    a =  a + a << 2 + a << 4
    a =  a $ a >> 28
    a =  a + a << 31
    return a
end

function hash_64_32(n::Uint64)
    local a::Uint64 = n
    a = ~a + a << 18
    a =  a $ a >> 31
    a =  a * 21
    a =  a $ a >> 11
    a =  a + a << 6
    a =  a $ a >> 22
    return uint32(a)
end

function hash_32_32(n::Uint32)
    local a::Uint32 = n
    a = a + 0x7ed55d16 + a << 12
    a = a $ 0xc761c23c $ a >> 19
    a = a + 0x165667b1 + a << 5
    a = a + 0xd3a2646c $ a << 9
    a = a + 0xfd7046c5 + a << 3
    a = a $ 0xb55a4f09 $ a >> 16
    return a
end

if Uint == Uint64
    hash_uint64(x::Uint64) = hash_64_64(x)
    hash_uint(x::Uint)     = hash_64_64(x)
else
    hash_uint64(x::Uint64) = hash_64_32(x)
    hash_uint(x::Uint)     = hash_32_32(x)
end

## hashing small, built-in numeric types ##

hx(a::Uint64, b::Float64, h::Uint) = hash_uint64((3a + reinterpret(Uint64,b)) - h)
const hx_NaN = hx(uint64(0), NaN, uint(0  ))

hash(x::Uint64,  h::Uint) = hx(x, float64(x), h)
hash(x::Int64,   h::Uint) = hx(reinterpret(Uint64,x), float64(x), h)
hash(x::Float64, h::Uint) = isnan(x) ? (hx_NaN $ h) : hx(box(Uint64,fptosi(unbox(Float64,x))), x, h)

hash(x::Union(Bool,Char,Int8,Uint8,Int16,Uint16,Int32,Uint32), h::Uint) = hash(int64(x), h)
hash(x::Float32, h::Uint) = hash(float64(x), h)

## hashing complex numbers ##

const h_imag = uint(0x32a7a07f3e7cd1f9)
const hash_0_imag = hash(0, h_imag)

function hash(z::Complex, h::Uint)
    # TODO: with default argument specialization, this would be better:
    # hash(real(z), h $ hash(imag(z), h $ h_imag) $ hash(0, h $ h_imag))
    hash(real(z), h $ hash(imag(z), h_imag) $ hash_0_imag)
end

## symbol & expression hashing ##

hash(x::Symbol, h::Uint) = hash(object_id(x), h)
hash(x::Expr, h::Uint) = hash(x.args, hash(x.head, h + uint(0x83c7900696d26dc6)))

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"""
```
Simulate
```
A Julia package for discrete event simulation.

The current stable, registered version is installed with
```julia
pkg> add Simulate
```

The development version is installed with:
```julia
pkg> add("https://github.com/pbayer/Simulate.jl")
```
"""
module Simulate

"""
    version

Gives the package version:

```jldoctest
julia> using Simulate

julia> Simulate.version
v"0.3.0"
```
"""
const version = v"0.3.0"

using Unitful, Random, DataStructures, Logging, .Threads
import Unitful: FreeUnits, Time

include("types.jl")
include("components.jl")
include("events.jl")
include("schedule.jl")
include("clock.jl")
include("process.jl")
include("threads.jl")
include("timer.jl")
include("printing.jl")

export  Clock, ActiveClock, RTClock, RTC, setUnit!, 𝐶,
        Timing, at, after, every, before, until, tau, sample_time!,
        Action, Fun, event!, periodic!,
        incr!, run!, stop!, resume!, reset!, sync!,
        Prc, process!, interrupt!, delay!, wait!, now!,
        PClock, fork!, collapse!, pclock, diag

Random.seed!(123)
rng = MersenneTwister(2020)
𝐶.state == Undefined() ? init!(𝐶) : nothing

end # module


text:
@testset "Array-related tests" begin show_testset_section() #Scope for test data

@testset "pad*()" begin show_testset_description()
	v = [4,5,6,10]
	@test padt(v, 7) == [0,0,0, 4,5,6,10]
	@test padl(v, 7) == [0,0,0, 4,5,6,10]
	@test padb(v, 7) == [4,5,6,10, 0,0,0]
	@test padr(v, 7) == [4,5,6,10, 0,0,0]
	v = [4 5 6; 8 10 11]
	@test padt(v, 5) == [0 0 0; 0 0 0; 0 0 0; 4 5 6; 8 10 11]
	@test padb(v, 5) == [4 5 6; 8 10 11; 0 0 0; 0 0 0; 0 0 0]
	@test padl(v, 7) == [0 0 0 0 4 5 6; 0 0 0 0 8 10 11]
	@test padr(v, 7) == [4 5 6 0 0 0 0; 8 10 11 0 0 0 0]
end

end


text:
function DSFactory(arg0::jint, arg1::jint)
    return DSFactory((jint, jint), arg0, arg1)
end

function build(obj::DSFactory, arg0::Vector{jdouble})
    return jcall(obj, "build", DerivativeStructure, (Vector{jdouble},), arg0)
end

function constant(obj::DSFactory, arg0::jdouble)
    return jcall(obj, "constant", DerivativeStructure, (jdouble,), arg0)
end

function get_compiler(obj::DSFactory)
    return jcall(obj, "getCompiler", DSCompiler, ())
end

function get_derivative_field(obj::DSFactory)
    return jcall(obj, "getDerivativeField", Field, ())
end

function variable(obj::DSFactory, arg0::jint, arg1::jdouble)
    return jcall(obj, "variable", DerivativeStructure, (jint, jdouble), arg0, arg1)
end


text:
using SimpleSDMDataSources, Test, Dates

include("worldclim-bioclim.jl")
include("worldclim-weather.jl")
include("earthenv-landcover.jl")
include("earthenv-heterogeneity.jl")
include("chelsa-bioclim.jl")
include("alwb.jl")
include("awap.jl")


text:
using Flux:epseltype
using Flux.Losses:crossentropy, kldivergence
using Statistics:mean
using Transformers:logcrossentropy, logkldivergence

"""
Return the loss for a neural machine translation model,
where ŷ is given as log().
If α ∈ (0,1), return KL-divergence between ŷ and label-smoothed y.
Otherwise, return cross entropy between ŷ and y without label smoothing.
"""
function logtranslationloss(ŷ, y; α=0.0, dims=1, agg=mean, ϵ=epseltype(ŷ))
    return translationloss(exp.(ŷ), y, α=α, dims=dims, agg=agg, ϵ=ϵ)
end

"""
Return the loss for a neural machine translation model.
If α ∈ (0,1), return KL-divergence between ŷ and label-smoothed y.
Otherwise, return cross entropy between ŷ and y without label smoothing.
"""
function translationloss(ŷ, y; α=0.0, dims=1, agg=mean, ϵ=epseltype(ŷ))
    if 0 < α < 1
        y = label_smoothing(y, α, dims=dims)
        return kldivergence(ŷ, y, dims=dims, agg=agg, ϵ=ϵ)
    else
        return crossentropy(y, ŷ, dims=dims, agg=agg, ϵ=ϵ)
    end
end

"""
Smooth ground proof labels.
Copied from https://github.com/FluxML/Flux.jl/blob/7a9a5eef2532d7029a960f7209b9a2e2fc76c29d/src/losses/functions.jl#L135-L147
which was released under MIT License.
"""
function label_smoothing(y::Union{AbstractArray,Number}, α::Number; dims::Int=1)
    if !(0 < α < 1)
        throw(ArgumentError("α must be between 0 and 1"))
    end
    if dims == 0
        y_smoothed = y .* (1 - α) .+ α * 1 // 2
    elseif dims == 1
        y_smoothed = y .* (1 - α) .+ α * 1 // size(y, 1)
    else
        throw(ArgumentError("`dims` should be either 0 or 1"))
    end
    return y_smoothed
end


text:
@testset "chainrules.jl" begin

    function FiniteDifferences.to_vec(k::AxisKeys.KaNda)
        v, b = to_vec(k.data)
        back(x) = wrapdims(b(x); AxisKeys.named_axiskeys(k)...)
        return v, back
    end

    function FiniteDifferences.to_vec(k::KeyedArray)
        v, b = to_vec(k.data)
        back(x) = wrapdims(b(x), axiskeys(k)...)
        return v, back
    end

    @testset "ProjectTo" begin
        data = rand(3)
        ka = wrapdims(data, a=1:3)
        p = ProjectTo(ka)
        @test p(data) == ka

        data = rand(3, 4)
        ka = wrapdims(data, a=1:3, b='a':'d')
        p = ProjectTo(ka)
        @test p(data) == ka
    end

    @testset "KeyedVector" begin
        data = rand(3)
        test_rrule(AxisKeys.keyless_unname, wrapdims(data, a=1:3); check_inferred=false)
        test_rrule(AxisKeys.keyless_unname, wrapdims(data, 1:3); check_inferred=false)
        test_rrule(AxisKeys.keyless_unname, data; check_inferred=false)

        # with matrix output tangent
        test_rrule(AxisKeys.keyless_unname, wrapdims(data, a=1:3); output_tangent=rand(3, 1), check_inferred=false)
        test_rrule(AxisKeys.keyless_unname, wrapdims(data, 1:3); output_tangent=rand(3, 1), check_inferred=false)
        test_rrule(AxisKeys.keyless_unname, data; output_tangent=rand(3, 1), check_inferred=false)
    end

    @testset "KeyedMatrix" begin
        data = rand(3, 4)
        test_rrule(AxisKeys.keyless_unname, wrapdims(data, a=1:3, b='a':'d'); check_inferred=false)
        test_rrule(AxisKeys.keyless_unname, wrapdims(data, 1:3, 'a':'d'); check_inferred=false)
        test_rrule(AxisKeys.keyless_unname, data; check_inferred=false)
    end

end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    get_juldate() -> julian_days

### Purpose ###

Return the number of Julian days for current time.

### Explanation ###

Return for current time the number of Julian calendar days since epoch
`-4713-11-24T12:00:00` as a floating point.

### Example ###

``` julia
get_juldate()
daycnv(get_juldate())
```

### Notes ###
Use `jdcnv` to get the number of Julian days for a different date.
"""
get_juldate() = jdcnv(now(Dates.UTC))


text:
using Test
using LinearAlgebra
using SparseIR
using SparseIR._LinAlg
using MultiFloats

@testset "_linalg.jl" begin
    @testset "jacobi with T = $T" for T in (Float64, Float64x2)
        A = T.(randn(20, 10))
        U, S, V = svd_jacobi(A)
        @test U * Diagonal(S) * V' ≈ A
    end

    @testset "rrqr with T = $T" for T in (Float64, Float64x2)
        A = T.(randn(40, 30))
        A_eps = norm(A) * eps(eltype(A))
        A_qr, A_rank = rrqr(A)
        A_rec = A_qr.Q * A_qr.R * A_qr.P'
        @test isapprox(A_rec, A; rtol=0, atol=4 * A_eps)
        @test A_rank == 30
    end

    @testset "rrqr_trunc with T = $T" for T in (Float64, Float64x2)
        # Vandermonde matrix
        A = Vector{T}(-1:0.02:1) .^ Vector(0:20)'
        m, n = size(A)
        A_qr, k = SparseIR._LinAlg.rrqr(A; rtol=1e-5)
        @test k < min(m, n)

        Q, R = SparseIR._LinAlg.truncate_qr_result(A_qr, k)
        A_rec = Q * R * A_qr.P'
        @test isapprox(A, A_rec, rtol=0, atol=1e-5 * norm(A))
    end

    @testset "tsvd with T = $T" for T in (Float64, Float64x2), tol in (1e-14, 1e-13)
        A = Vector{T}(-1:0.01:1) .^ Vector(0:50)'
        U, S, V = SparseIR._LinAlg.tsvd(A; rtol=tol)
        k = length(S)

        @test U * Diagonal(S) * V' ≈ A rtol = 0 atol = tol * norm(A)
        @test U'U ≈ I
        @test V'V ≈ I
        @test issorted(S; rev=true)
        @test k < minimum(size(A))

        A_svd = svd(Float64.(A))
        @test S ≈ A_svd.S[1:k]
    end
end


text:
################################################################################
#                                 Euler Angles
################################################################################

export angle_to_angleaxis, angle_to_angle, angle_to_dcm, angle_to_quat
export smallangle_to_dcm, smallangle_to_quat
export angle_to_rot, smallangle_to_rot
export inv

################################################################################
#                                  Operations
################################################################################

"""
    *(Θ₂::EulerAngles, Θ₁::EulerAngles)

Compute the composed rotation of `Θ₁ -> Θ₂`. Notice that the rotation will be
represented by Euler angles (see `EulerAngles`) with the same rotation sequence
as `Θ₂`.

"""
@inline function *(Θ₂::EulerAngles, Θ₁::EulerAngles)
    # Convert to quaternions, compute the composition, and convert back to Euler
    # angles.
    q₁ = angle_to_quat(Θ₁)
    q₂ = angle_to_quat(Θ₂)

    quat_to_angle(q₁*q₂, Θ₂.rot_seq)
end

"""
    inv(Θ::EulerAngles)

Return the Euler angles that represent the inverse rotation of `Θ`. Notice that
the rotation sequence of the result will be the inverse of the input. Hence, if
the input rotation sequence is, for example, `:XYZ`, then the result will be
represented using `:ZYX`.

"""
function inv(Θ::EulerAngles)
    # Check what will be the inverse rotation.
    if Θ.rot_seq == :XYZ
        inv_rot_seq = :ZYX
    elseif Θ.rot_seq == :XZY
        inv_rot_seq = :YZX
    elseif Θ.rot_seq == :YXZ
        inv_rot_seq = :ZXY
    elseif Θ.rot_seq == :YZX
        inv_rot_seq = :XZY
    elseif Θ.rot_seq == :ZXY
        inv_rot_seq = :YXZ
    elseif Θ.rot_seq == :ZYX
        inv_rot_seq = :XYZ
    else
        inv_rot_seq = Θ.rot_seq
    end

    # Return the Euler angle that represented the inverse rotation.
    EulerAngles(-Θ.a3, -Θ.a2, -Θ.a1, inv_rot_seq)
end

################################################################################
#                                      IO
################################################################################

"""
    show(io::IO, Θ::EulerAngles{T}) where T
    show(io::IO, mime::MIME"text/plain", Θ::EulerAngles{T}) where T

Print the Euler angles `Θ` to the IO `io`.

"""
function show(io::IO, Θ::EulerAngles{T}) where T
    str     = @sprintf "%8.4f %8.4f %8.4f rad" Θ.a1 Θ.a2 Θ.a3
    rot_seq = String(Θ.rot_seq)

    print(io, "EulerAngles{$T}:")
    print(io,   " R($rot_seq): " * str)

    nothing
end

function show(io::IO, mime::MIME"text/plain", Θ::EulerAngles{T}) where T
    # Check if the user wants colors.
    color = get(io, :color, false)

    d = (color) ? _d : ""
    g = (color) ? _g : ""
    y = (color) ? _y : ""
    u = (color) ? _u : ""

    str_θ₁  = @sprintf "%8.4f rad (%9.4f deg)" Θ.a1 rad2deg(Θ.a1)
    str_θ₂  = @sprintf "%8.4f rad (%9.4f deg)" Θ.a2 rad2deg(Θ.a2)
    str_θ₃  = @sprintf "%8.4f rad (%9.4f deg)" Θ.a3 rad2deg(Θ.a3)
    rot_seq = String(Θ.rot_seq)

    println(io, "EulerAngles{$T}:")
    println(io,   "$g  R($(rot_seq[1])): $d" * str_θ₁)
    println(io,   "$y  R($(rot_seq[2])): $d" * str_θ₂)
      print(io,   "$u  R($(rot_seq[3])): $d" * str_θ₃)

    nothing
end

################################################################################
#                                 Conversions
################################################################################

# Euler Angle and Axis
# ==============================================================================

"""
    angle_to_angleaxis(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq::Symbol = :ZYX)
    angle_to_angleaxis(Θ::EulerAngles)

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq` to an Euler angle and axis representation.  Those values can also be
passed inside the structure `Θ` (see `EulerAngles`).

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```julia-repl
julia> angle_to_angleaxis(1,0,0,:XYZ)
EulerAngleAxis{Float64}:
  Euler angle:   1.0000 rad ( 57.2958 deg)
   Euler axis: [  1.0000,   0.0000,   0.0000]

julia> Θ = EulerAngles(1,1,1,:XYZ);

julia> angle_to_angleaxis(Θ)
EulerAngleAxis{Float64}:
  Euler angle:   1.9391 rad (111.1015 deg)
   Euler axis: [  0.6924,   0.2031,   0.6924]

```
"""
@inline function angle_to_angleaxis(θ₁::Number, θ₂::Number, θ₃::Number,
                                    rot_seq::Symbol = :ZYX)
    # First convert to quaternion and then to Euler angle and axis.
    quat_to_angleaxis( angle_to_quat( θ₁, θ₂, θ₃, rot_seq ) )
end

@inline angle_to_angleaxis(Θ::EulerAngles) =
    angle_to_angleaxis(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)

# Euler Angles
# ==============================================================================

"""
    angle_to_angle(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq_orig::Symbol, rot_seq_dest::Symbol)
    angle_to_angle(Θ::EulerAngles, rot_seq_dest::Symbol)

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq_orig` to a new set of Euler angles with rotation sequence
`rot_seq_dest`. The input values of the origin Euler angles can also be passed
inside the structure `Θ` (see `EulerAngles`).

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```julia-repl
julia> angle_to_angle(-pi/2, -pi/3, -pi/4, :ZYX, :XYZ)
EulerAngles{Float64}:
  R(X):  -1.0472 rad ( -60.0000 deg)
  R(Y):   0.7854 rad (  45.0000 deg)
  R(Z):  -1.5708 rad ( -90.0000 deg)

julia> angle_to_angle(-pi/2, 0, 0, :ZYX, :XYZ)
EulerAngles{Float64}:
  R(X):   0.0000 rad (   0.0000 deg)
  R(Y):   0.0000 rad (   0.0000 deg)
  R(Z):  -1.5708 rad ( -90.0000 deg)

julia> Θ = EulerAngles(1,2,3,:XYX)
EulerAngles{Int64}:
  R(X):   1.0000 rad (  57.2958 deg)
  R(Y):   2.0000 rad ( 114.5916 deg)
  R(X):   3.0000 rad ( 171.8873 deg)

julia> angle_to_angle(Θ,:ZYZ)
EulerAngles{Float64}:
  R(Z):  -2.7024 rad (-154.8356 deg)
  R(Y):   1.4668 rad (  84.0393 deg)
  R(Z):  -1.0542 rad ( -60.3984 deg)
```
"""
@inline angle_to_angle(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq_orig::Symbol,
                       rot_seq_dest::Symbol) =
    quat_to_angle(angle_to_quat(θ₁, θ₂, θ₃, rot_seq_orig), rot_seq_dest)

@inline angle_to_angle(Θ::EulerAngles, rot_seq_dest::Symbol) =
    angle_to_angle(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq, rot_seq_dest)

# Direction Cosine Matrix
# ==============================================================================

"""
    angle_to_dcm(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq::Symbol = :ZYX)

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq` to a direction cosine matrix.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Remarks

This function assigns `dcm = A3 * A2 * A1` in which `Ai` is the DCM related with
the *i*-th rotation, `i Є [1,2,3]`.

# Example

```julia-repl
dcm = angle_to_dcm(pi/2, pi/3, pi/4, :ZYX)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553
```

"""
function angle_to_dcm(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq::Symbol = :ZYX)

    # Compute the sines and cosines.
    s₁, c₁ = sincos(θ₁)
    s₂, c₂ = sincos(θ₂)
    s₃, c₃ = sincos(θ₃)

    if rot_seq == :ZYX
        return DCM(       c₂*c₁    ,        c₂*s₁    , -s₂  ,
                   s₃*s₂*c₁ - c₃*s₁, s₃*s₂*s₁ + c₃*c₁, s₃*c₂,
                   c₃*s₂*c₁ + s₃*s₁, c₃*s₂*s₁ - s₃*c₁, c₃*c₂)'

    elseif rot_seq == :XYX
        return DCM( c₂  ,         s₁*s₂    ,       -c₁*s₂    ,
                   s₂*s₃, -s₁*c₂*s₃ + c₁*c₃, c₁*c₂*s₃ + s₁*c₃,
                   s₂*c₃, -s₁*c₃*c₂ - c₁*s₃, c₁*c₃*c₂ - s₁*s₃)'

    elseif rot_seq == :XYZ
        return DCM( c₂*c₃,  s₁*s₂*c₃ + c₁*s₃, -c₁*s₂*c₃ + s₁*s₃,
                   -c₂*s₃, -s₁*s₂*s₃ + c₁*c₃,  c₁*s₂*s₃ + s₁*c₃,
                      s₂ ,        -s₁*c₂    ,       c₁*c₂)'

    elseif rot_seq == :XZX
        return DCM(   c₂ ,         c₁*s₂    ,         s₁*s₂    ,
                   -s₂*c₃,  c₁*c₃*c₂ - s₁*s₃,  s₁*c₃*c₂ + c₁*s₃,
                    s₂*s₃, -c₁*c₂*s₃ - s₁*c₃, -s₁*c₂*s₃ + c₁*c₃)'

    elseif rot_seq == :XZY
        return DCM(c₃*c₂, c₁*c₃*s₂ + s₁*s₃, s₁*c₃*s₂ - c₁*s₃,
                    -s₂ ,        c₁*c₂    ,        s₁*c₂    ,
                   s₃*c₂, c₁*s₂*s₃ - s₁*c₃, s₁*s₂*s₃ + c₁*c₃)'

    elseif rot_seq == :YXY
        return DCM(-s₁*c₂*s₃ + c₁*c₃, s₂*s₃ , -c₁*c₂*s₃ - s₁*c₃,
                          s₁*s₂    ,   c₂  ,         c₁*s₂    ,
                   s₁*c₃*c₂ + c₁*s₃, -s₂*c₃,  c₁*c₃*c₂ - s₁*s₃)'

    elseif rot_seq == :YXZ
        return DCM( c₁*c₃ + s₂*s₁*s₃, c₂*s₃, -s₁*c₃ + s₂*c₁*s₃,
                   -c₁*s₃ + s₂*s₁*c₃, c₂*c₃,  s₁*s₃ + s₂*c₁*c₃,
                        s₁*c₂       ,  -s₂ ,      c₂*c₁       )'

    elseif rot_seq == :YZX
        return DCM(        c₁*c₂    ,    s₂ ,        -s₁*c₂    ,
                   -c₃*c₁*s₂ + s₃*s₁,  c₂*c₃,  c₃*s₁*s₂ + s₃*c₁,
                    s₃*c₁*s₂ + c₃*s₁, -s₃*c₂, -s₃*s₁*s₂ + c₃*c₁)'

    elseif rot_seq == :YZY
        return DCM(c₁*c₃*c₂ - s₁*s₃, s₂*c₃, -s₁*c₃*c₂ - c₁*s₃,
                         -c₁*s₂    ,   c₂ ,         s₁*s₂    ,
                   c₁*c₂*s₃ + s₁*c₃, s₂*s₃, -s₁*c₂*s₃ + c₁*c₃)'

    elseif rot_seq == :ZXY
        return DCM(c₃*c₁ - s₂*s₃*s₁, c₃*s₁ + s₂*s₃*c₁, -s₃*c₂,
                      -c₂*s₁       ,     c₂*c₁       ,    s₂ ,
                   s₃*c₁ + s₂*c₃*s₁, s₃*s₁ - s₂*c₃*c₁,  c₂*c₃)'

    elseif rot_seq == :ZXZ
        return DCM(-s₁*c₂*s₃ + c₁*c₃, c₁*c₂*s₃ + s₁*c₃, s₂*s₃,
                   -s₁*c₃*c₂ - c₁*s₃, c₁*c₃*c₂ - s₁*s₃, s₂*c₃,
                           s₁*s₂    ,       -c₁*s₂    ,  c₂)'

    elseif rot_seq == :ZYZ
        return DCM( c₁*c₃*c₂ - s₁*s₃,  s₁*c₃*c₂ + c₁*s₃, -s₂*c₃,
                   -c₁*c₂*s₃ - s₁*c₃, -s₁*c₂*s₃ + c₁*c₃,  s₂*s₃,
                           c₁*s₂    ,         s₁*s₂    ,    c₂)'
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

"""
    angle_to_dcm(Θ::EulerAngles)

Convert the Euler angles `Θ` (see `EulerAngles`) to a direction cosine matrix.

# Returns

The direction cosine matrix.

# Remarks

This function assigns `dcm = A3 * A2 * A1` in which `Ai` is the DCM related with
the *i*-th rotation, `i Є [1,2,3]`.

# Example

```julia-repl
julia> angle_to_dcm(EulerAngles(pi/2, pi/3, pi/4, :ZYX))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553
```

"""
angle_to_dcm(Θ::EulerAngles) = angle_to_dcm(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)

"""
    smallangle_to_dcm(θx::Number, θy::Number, θz::Number; normalize = true)

Create a direction cosine matrix from three small rotations of angles `θx`,
`θy`, and `θz` [rad] about the axes X, Y, and Z, respectively. If the keyword
`normalize` is `true`, then the matrix will be normalized using the function
`orthonormalize`.

# Example

```julia-repl
julia> smallangle_to_dcm(+0.01, -0.01, +0.01)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.9999     0.00989903  0.010098
 -0.009999   0.999901    0.00989802
 -0.009999  -0.009998    0.9999

julia> smallangle_to_dcm(+0.01, -0.01, +0.01; normalize = false)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  1.0    0.01  0.01
 -0.01   1.0   0.01
 -0.01  -0.01  1.0
```

"""
@inline function smallangle_to_dcm(θx::Number, θy::Number, θz::Number;
                                   normalize = true)
    D = DCM(  1, +θz, -θy,
            -θz,   1, +θx,
            +θy, -θx,   1)'

    normalize ? orthonormalize(D) : D
end

# Quaternion
# ==============================================================================

"""
    angle_to_quat(θ₁::T1, θ₂::T2, θ₃::T3, rot_seq::Symbol = :ZYX) where {T1<:Number, T2<:Number, T3<:Number}

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq` to a quaternion.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Remarks

This function assigns `q = q1 * q2 * q3` in which `qi` is the quaternion related
with the *i*-th rotation, `i Є [1,2,3]`.

# Example

```julia-repl
julia> angle_to_quat(pi/2, pi/3, pi/4, :ZYX)
Quaternion{Float64}:
  + 0.7010573846499779 - 0.09229595564125714.i + 0.5609855267969309.j + 0.43045933457687935.k
```

"""
function angle_to_quat(θ₁::T1, θ₂::T2, θ₃::T3, rot_seq::Symbol = :ZYX) where
    {T1<:Number, T2<:Number, T3<:Number}

    T = promote_type(T1,T2,T3)

    # Compute the sines and cosines of half angle.
    s₁, c₁ = sincos(T(θ₁)/2)
    s₂, c₂ = sincos(T(θ₂)/2)
    s₃, c₃ = sincos(T(θ₃)/2)

    if rot_seq == :ZYX
        q0 = c₁*c₂*c₃ + s₁*s₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*c₂*s₃ - s₁*s₂*c₃),
                          s*(c₁*s₂*c₃ + s₁*c₂*s₃),
                          s*(s₁*c₂*c₃ - c₁*s₂*s₃))
    elseif rot_seq == :XYX
        q0 = c₁*c₂*c₃ - s₁*c₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*c₂*s₃ + s₁*c₂*c₃),
                          s*(c₁*s₂*c₃ + s₁*s₂*s₃),
                          s*(s₁*s₂*c₃ - c₁*s₂*s₃))
    elseif rot_seq == :XYZ
        q0 = c₁*c₂*c₃ - s₁*s₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(s₁*c₂*c₃ + c₁*s₂*s₃),
                          s*(c₁*s₂*c₃ - s₁*c₂*s₃),
                          s*(c₁*c₂*s₃ + s₁*s₂*c₃))
    elseif rot_seq == :XZX
        q0 = c₁*c₂*c₃ - s₁*c₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*c₂*s₃ + s₁*c₂*c₃),
                          s*(c₁*s₂*s₃ - s₁*s₂*c₃),
                          s*(c₁*s₂*c₃ + s₁*s₂*s₃))
    elseif rot_seq == :XZY
        q0 = c₁*c₂*c₃ + s₁*s₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(s₁*c₂*c₃ - c₁*s₂*s₃),
                          s*(c₁*c₂*s₃ - s₁*s₂*c₃),
                          s*(c₁*s₂*c₃ + s₁*c₂*s₃))
    elseif rot_seq == :YXY
        q0 = c₁*c₂*c₃ - s₁*c₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*s₂*c₃ + s₁*s₂*s₃),
                          s*(c₁*c₂*s₃ + s₁*c₂*c₃),
                          s*(c₁*s₂*s₃ - s₁*s₂*c₃))
    elseif rot_seq == :YXZ
        q0 = c₁*c₂*c₃ + s₁*s₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*s₂*c₃ + s₁*c₂*s₃),
                          s*(s₁*c₂*c₃ - c₁*s₂*s₃),
                          s*(c₁*c₂*s₃ - s₁*s₂*c₃))
    elseif rot_seq == :YZX
        q0 = c₁*c₂*c₃ - s₁*s₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*c₂*s₃ + s₁*s₂*c₃),
                          s*(s₁*c₂*c₃ + c₁*s₂*s₃),
                          s*(c₁*s₂*c₃ - s₁*c₂*s₃))
    elseif rot_seq == :YZY
        q0 = c₁*c₂*c₃ - s₁*c₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(s₁*s₂*c₃ - c₁*s₂*s₃),
                          s*(c₁*c₂*s₃ + s₁*c₂*c₃),
                          s*(c₁*s₂*c₃ + s₁*s₂*s₃))
    elseif rot_seq == :ZXY
        q0 = c₁*c₂*c₃ - s₁*s₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*s₂*c₃ - s₁*c₂*s₃),
                          s*(c₁*c₂*s₃ + s₁*s₂*c₃),
                          s*(s₁*c₂*c₃ + c₁*s₂*s₃))
    elseif rot_seq == :ZXZ
        q0 = c₁*c₂*c₃ - s₁*c₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*s₂*c₃ + s₁*s₂*s₃),
                          s*(s₁*s₂*c₃ - c₁*s₂*s₃),
                          s*(c₁*c₂*s₃ + s₁*c₂*c₃))
    elseif rot_seq == :ZYZ
        q0 = c₁*c₂*c₃ - s₁*c₂*s₃

        s = (q0 < 0) ? -1 : +1

        return Quaternion(s*q0,
                          s*(c₁*s₂*s₃ - s₁*s₂*c₃),
                          s*(c₁*s₂*c₃ + s₁*s₂*s₃),
                          s*(c₁*c₂*s₃ + s₁*c₂*c₃))
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

"""
    angle_to_quat(eulerang::EulerAngles)

Convert the Euler angles `eulerang` (see `EulerAngles`) to a quaternion.

# Remarks

This function assigns `q = q1 * q2 * q3` in which `qi` is the quaternion related
with the *i*-th rotation, `i Є [1,2,3]`.

# Example

```julia-repl
julia> angle_to_quat(pi/2, pi/3, pi/4, :ZYX)
Quaternion{Float64}:
  + 0.7010573846499779 - 0.09229595564125714.i + 0.5609855267969309.j + 0.43045933457687935.k
```

"""
angle_to_quat(Θ::EulerAngles) = angle_to_quat(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)

"""
    smallangle_to_quat(θx::Number, θy::Number, θz::Number)

Create a quaternion from three small rotations of angles `θx`, `θy`, and `θz`
[rad] about the axes X, Y, and Z, respectively.

# Remarks

The quaternion is normalized.

# Example

```julia-repl
julia> smallangle_to_quat(+0.01, -0.01, +0.01)
Quaternion{Float64}:
  + 0.9999625021092433 + 0.004999812510546217.i - 0.004999812510546217.j + 0.004999812510546217.k
```

"""
function smallangle_to_quat(θx::Number, θy::Number, θz::Number)
    q0     = 1
    q1     = θx/2
    q2     = θy/2
    q3     = θz/2
    norm_q = sqrt(q0^2 + q1^2 + q2^2 + q3^2)

    Quaternion(q0/norm_q, q1/norm_q, q2/norm_q, q3/norm_q)
end

################################################################################
#                                     API
################################################################################

"""
    angle_to_rot([T,] θx::Number, θy::Number, θz::Number, rot_seq::Symbol)

Convert the Euler angles `Θx`, `Θy`, and `Θz` [rad] with the rotation sequence
`rot_seq` to a rotation description of type `T`, which can be `DCM` or
`Quaternion`. If the type `T` is not specified, then it defaults to `DCM`.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```julia-repl
julia> dcm = angle_to_rot(pi/2, pi/3, pi/4, :ZYX)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553

julia> q   = angle_to_rot(Quaternion,pi/2, pi/3, pi/4, :ZYX)
Quaternion{Float64}:
  + 0.7010573846499779 - 0.09229595564125714.i + 0.5609855267969309.j + 0.43045933457687935.k
```

"""
@inline angle_to_rot(θx::Number, θy::Number, θz::Number, rot_seq::Symbol) =
    angle_to_dcm(θx, θy, θz, rot_seq)

@inline angle_to_rot(::Type{DCM},
                     θx::Number, θy::Number, θz::Number, rot_seq::Symbol) =
    angle_to_dcm(θx, θy, θz, rot_seq)

@inline angle_to_rot(::Type{Quaternion},
                     θx::Number, θy::Number, θz::Number, rot_seq::Symbol) =
    angle_to_quat(θx, θy, θz, rot_seq)

"""
    angle_to_rot([T,] Θ::EulerAngles)

Convert the Euler angles `Θ` (see `EulerAngles`) to a rotation description of
type `T`, which can be `DCM` or `Quaternion`. If the type `T` is not specified,
then it defaults to `DCM`.

# Example

```julia-repl
julia> dcm = angle_to_rot(EulerAngles(pi/2, pi/3, pi/4, :ZYX))
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553

julia> q   = angle_to_rot(Quaternion,EulerAngles(pi/2, pi/3, pi/4, :ZYX))
Quaternion{Float64}:
  + 0.7010573846499779 - 0.09229595564125714.i + 0.5609855267969309.j +
  0.43045933457687935.k
```

"""
@inline angle_to_rot(Θ::EulerAngles) = angle_to_dcm(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)

@inline angle_to_rot(::Type{DCM}, Θ::EulerAngles) =
    angle_to_dcm(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)

@inline angle_to_rot(::Type{Quaternion}, Θ::EulerAngles) =
    angle_to_quat(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)

"""
    smallangle_to_rot([T,] θx::Number, θy::Number, θz::Number[; normalize = true])

Create a rotation description of type `T` from three small rotations of angles
`θx`, `θy`, and `θz` [rad] about the axes X, Y, and Z, respectively.

The type `T` of the rotation description can be `DCM` or `Quaternion`. If the
type `T` is not specified, then if defaults to `DCM`.

If `T` is `DCM`, then the resulting matrix will be orthonormalized using the
`orthonormalize` function if the keyword `normalize` is `true`.

# Example

```julia-repl
julia> dcm = smallangle_to_rot(+0.01, -0.01, +0.01)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.9999     0.00989903  0.010098
 -0.009999   0.999901    0.00989802
 -0.009999  -0.009998    0.9999

julia> dcm = smallangle_to_rot(+0.01, -0.01, +0.01; normalize = false)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  1.0    0.01  0.01
 -0.01   1.0   0.01
 -0.01  -0.01  1.0

julia> q   = smallangle_to_rot(Quaternion,+0.01, -0.01, +0.01)
Quaternion{Float64}:
  + 0.9999625021092433 + 0.004999812510546217.i - 0.004999812510546217.j + 0.004999812510546217.k
```
"""
@inline smallangle_to_rot(θx::Number, θy::Number, θz::Number; normalize = true) =
    smallangle_to_dcm(θx, θy, θz; normalize = normalize)

@inline smallangle_to_rot(::Type{DCM}, θx::Number, θy::Number, θz::Number;
                          normalize = true) =
    smallangle_to_dcm(θx, θy, θz; normalize = normalize)

@inline smallangle_to_rot(::Type{Quaternion}, θx::Number, θy::Number, θz::Number) =
    smallangle_to_quat(θx, θy, θz)


text:
module GeometricMultigrid

include("FieldVector.jl")
export FieldVector, @loop

include("FieldMatrix.jl")
export FieldMatrix, Poisson

include("SolveState.jl")
export SolveState

include("MultiGrid.jl")
export mg_state, mg, mg!

end


text:
import JSON

# ------------------------------------
# Commenting out the json file
# ------------------------------------
function make_json(PARAMDAT)
    fff = []
    open(PARAMDAT, "r") do f
        fff = read(f,String)
    end

    fff = replace(fff,r"#(.+)\n" => "\n")

    read_PARAMDAT = "read_"*PARAMDAT
    open(read_PARAMDAT, "w") do f
        write(f,fff)
    end
    return read_PARAMDAT
end

# ------------------------------------
# prepare reading the json file
# ------------------------------------
function read_json(read_PARAMDAT)
    dict=1
    open(read_PARAMDAT, "r") do f
        dicttxt = read(f,String)  # file information to string
        dict = JSON.parse(dicttxt)  # parse and transform data
    end
    return dict
end

# ------------------------------------
# read PARAMDAT
# ------------------------------------
function read_para(dict)
    # file name
    out_file_front = dict["out_file_front"]
    out_ext        = dict["out_ext"]

    # restart
    restartnum   = dict["restartnum"]
    restart_file = dict["Restart"] * restartnum * dict["in_ext"]
    restartnum   = Int(parse(Float64,restartnum))

    # output
    init_small　=　parse(Float64,dict["init_small"])
    norm_ok   　=　parse(Float64,dict["norm_ok"])

    # time step
    time_integ   = dict["time_integ"]
    nt           = Int(parse(Float64,dict["nt"]))
    dt           = parse(Float64,dict["dt"])
    every_outnum = Int(parse(Float64,dict["every_outnum"]))
    in_nt        = Int(parse(Float64,dict["inner_step"]))
    dtau         = parse(Float64,dict["dtau"])
    cfl          = parse(Float64,dict["CFL"])

    # scheme
    ad_scheme = Int(parse(Float64,dict["advection"]))

    # Initial values
    init_rho = parse(Float64,dict["init_rho"])
    init_u   = parse(Float64,dict["init_u"])
    init_v   = parse(Float64,dict["init_v"])
    init_w   = parse(Float64,dict["init_w"])
    init_p   = parse(Float64,dict["init_p"])
    init_T   = parse(Float64,dict["init_T"])

    # constant
    specific_heat_ratio = parse(Float64,dict["specific_heat_ratio"])
    Rd = parse(Float64,dict["Rd"])

    # boundary conditions
    #=
    bdcon[i][j]
    i:bd number (1,2,3..)
    j=1:bd number (0~4)
    j=2:rho
    j=3:u
    j=4:v
    j=5:p
    =#
    bdcon = []
    k = 1
    while true
        try
            bd=[parse(Int,dict["bd"*string(k)*"_con"]),   parse(Float64,dict["bd"*string(k)*"_rho"]),
                parse(Float64,dict["bd"*string(k)*"_u"]), parse(Float64,dict["bd"*string(k)*"_v"]), parse(Float64,dict["bd"*string(k)*"_w"]),
                parse(Float64,dict["bd"*string(k)*"_p"]), parse(Float64,dict["bd"*string(k)*"_T"])]
            push!(bdcon,bd)
        catch
            break
        end
        k+=1
    end

    return out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
            time_integ, nt, dt, every_outnum, in_nt, dtau, cfl, ad_scheme,
            init_rho, init_u, init_v, init_w, init_p, init_T, specific_heat_ratio, Rd, bdcon
end

function input_para(PARAMDAT)
    read_PARAMDAT = make_json(PARAMDAT)
    dict          = read_json(read_PARAMDAT)
    out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
    time_integ, nt, dt, every_outnum, in_nt, dtau, cfl, ad_scheme,
    init_rho, init_u, init_v, init_w, init_p, init_T, specific_heat_ratio, Rd, bdcon = read_para(dict)
    println("fin reading para")

    return out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
            time_integ, nt, dt, every_outnum, in_nt, dtau, cfl, ad_scheme,
            init_rho, init_u, init_v, init_w, init_p, init_T, specific_heat_ratio, Rd, bdcon
end


text:
#===============================================================================

    Defines:
        - DiffusionGlobalWorkspace      : global workspace for diffusions
        - DiffusionGlobalSubworkspace   : convenience object with containers
        - DiffusionLocalWorkspace       : local workspace for diffusions
        - DiffusionLocalSubworkspace    : convenience object with containers
        - LocalUpdtParamNames           : convenience object for setting params
    and methods for them.

===============================================================================#

#===============================================================================
                        Global workspace for diffusions
===============================================================================#

"""
    struct DiffusionGlobalSubworkspace{T,TGP,TW,TWn,TX} <: eMCMC.GlobalWorkspace{T}
        P::TGP
        WW::TW
        Wnr::TWn
        XX::TX
    end

Convenience sub-workspace that holds
- `P::Vector{<:Vector{<:GuidProp}}`: the `GuidProp` laws [TODO change to PP or
                                     PPP]
- `WW::Vector{<:Vector{<:Trajectory}}`: the containers for sampling the Wiener
                                        noise
- `Wnr::Vector{<:Wiener}`: the flags for the Wiener noise
- `XX::Vector{<:Vector{<:Trajectory}}`: the containers for sampling the
                                        diffusion paths
The outer array corresponds to successive `recording`s, the inner arrays
correspond to successive observations within each recording. Each
`DiffusionSubworkspace` has two `DiffusionGlobalSubworkspace`s, one for the
accepted law/path another for the proposal.

    DiffusionGlobalSubworkspace{T}(
        aux_laws, data::AllObservations, tts, args,
    ) where T

A standard constructor. `aux_laws` is a list of auxiliary laws (one for each
recording), `data` holds the observations, `tts` are the time-grids on which
to do imputation when sampling and `args` are additional positional arguments
passed to initializers of `GuidProp`.
"""
struct DiffusionGlobalSubworkspace{T,TGP,TW,TWn,TX} <: eMCMC.GlobalWorkspace{T}
    P::TGP
    WW::TW
    Wnr::TWn
    XX::TX

    function DiffusionGlobalSubworkspace{T}(
            aux_laws, data::AllObservations, tts, args,
        ) where T
        N = num_recordings(data)
        P = map(1:N) do i
            build_guid_prop(aux_laws[i], data.recordings[i], tts[i], args[i])
        end

        traj = map(1:N) do i
            trajectory(tts[i], data.recordings[i].P)
        end
        XX = map(x->x.process, traj)
        WW = map(x->x.wiener, traj)
        Wnr = [Wiener(data.recordings[i].P) for i=1:N]
        init_paths!(P, WW, Wnr, XX, data)
        new{T,typeof(P),typeof(WW),typeof(Wnr),typeof(XX)}(P, WW, Wnr, XX)
    end
end

"""
    init_paths!(P, WW, Wnr, XX, data)

Sample paths of guided proposals without using the preconditioned Crank–Nicolson
scheme.
"""
function init_paths!(P, WW, Wnr, XX, data)
    num_recordings = length(P)

    # this loop should be entirely parallelizable, as recordings
    # are—by design—conditionally independent
    for i in 1:num_recordings
        success = false
        while !success
            y1 = rand(data.recordings[i].x0_prior)
            success = forward_guide!(P[i], XX[i], WW[i], y1; Wnr=Wnr[i])
        end
    end
end

"""
    struct DiffusionGlobalWorkspace{T,SW,SWD,TD,TX} <: eMCMC.GlobalWorkspace{T}
        sub_ws::SW
        sub_ws_diff::SWD
        sub_ws_diff°::SWD
        data::TD
        stats::GenericChainStats{T}
        XX_buffer::TX
        pnames::Vector{Symbol}
    end

Global workspace for MCMC problems concerning diffusions.

# Arguments
- `sub_ws`: `GenericGlobalSubworkspace` from `ExtensibleMCMC` that stores info
            about the parameter chain
- `sub_ws_diff`: `DiffusionGlobalSubworkspace` containers for the accepted
                 diffusion law/paths
- `sub_ws_diff°`: `DiffusionGlobalSubworkspace` containers for the proposed
                  diffusion law/paths
- `data`: observations
- `stats::GenericChainStats`: generic online statistics for the param chain
- `XX_buffer`: cache for saving a chain of paths
- `pnames`: a list of parameter names that are being updated
"""
struct DiffusionGlobalWorkspace{T,SW,SWD,TD,TX} <: eMCMC.GlobalWorkspace{T}
    sub_ws::SW
    sub_ws_diff::SWD
    sub_ws_diff°::SWD
    data::TD
    stats::GenericChainStats{T}
    XX_buffer::TX
    pnames::Vector{Symbol}
end

function eMCMC.init_global_workspace(
        backend::DiffusionMCMCBackend,
        num_mcmc_steps,
        updates::Vector{<:eMCMC.MCMCUpdate},
        data,
        θinit::Nothing;
        kwargs...
    ) where T
    dummy = OrderedDict(
        :none => -Inf, # no parameters needed for smoothing
    )
    eMCMC.init_global_workspace(
        backend, num_mcmc_steps, updates, data, dummy; kwargs...
    )
end

function eMCMC.init_global_workspace(
        ::DiffusionMCMCBackend,
        num_mcmc_steps,
        updates::Vector{<:eMCMC.MCMCUpdate},
        data,
        θinit::OrderedDict{Symbol,T};
        kwargs...
    ) where T
    dkwargs = Dict(kwargs)
    M, NU = num_mcmc_steps, length(updates)

    _θ_init = collect(values(θinit))
    _pnames = collect(keys(θinit))

    sub_ws = eMCMC.StandardGlobalSubworkspace(M, NU, data, _θ_init)

    sub_ws_diff = DiffusionGlobalSubworkspace{T}(
        extract_aux_laws(updates, data),
        data,
        OBS.setup_time_grids(
            data,
            get(dkwargs, :dt, 0.01),
            get(dkwargs, :τ, _DEFAULT_TIME_TRANSFORM),
            get(dkwargs, :Ttime, Float64),
            get(dkwargs, :tt, nothing),
        ),
        package(get(dkwargs, :guid_prop_args, tuple()), data), # from ObservationSchemes
    )

    path_saving_buffer = PathSavingBuffer(
        sub_ws_diff.XX,
        get(dkwargs, :path_buffer_size, 1)
    )

    SW, SWD, TD = typeof(sub_ws), typeof(sub_ws_diff), typeof(data)
    TX = typeof(path_saving_buffer)

    DiffusionGlobalWorkspace{T,SW,SWD,TD,TX}(
        sub_ws,
        sub_ws_diff,
        deepcopy(sub_ws_diff),
        data,
        GenericChainStats(_θ_init, NU, M),
        path_saving_buffer,
        _pnames,
    )
end

function extract_aux_laws(updates::Vector{<:eMCMC.MCMCUpdate}, data)
    imps = filter(x->typeof(x)<:PathImputation, updates)
    N = length(imps)
    @assert N > 0
    aux_laws = map(x->OBS.package(x.aux_laws, data), imps)
    @assert N == 1 || all([aux_laws[1]==aux_laws[i] for i=2:N])
    aux_laws[1]
end

#===============================================================================
                        Local workspace for diffusions
===============================================================================#
"""
    struct DiffusionLocalSubworkspace{T,TP,TW,TWn,TX,Tz0} <: eMCMC.LocalWorkspace{T}
        ll::Vector{Float64}
        ll_history::Vector{Vector{Float64}}
        P::TP
        WW::TW
        Wnr::TWn
        XX::TX
        z0s::Tz0
    end

Convenience sub-workspace that holds
- `ll::Vector{Float64}`: computational buffer for log-likelihoods (one for each
                         recording)
- `ll_history::Vector{Vector{Float64}}`:  history of all log-likelihoods
- `P`: appropriately shaped views to containers with `GuidProp`s
- `WW`: appropriately shaped views to containers for sampled Wiener process
- `Wnr`: list of Wiener flags for each recording
- `XX`: appropriately shaped views to containers for sampled diffusion paths
- `z0s`: list of white noise for starting points (one for each recording)
"""
struct DiffusionLocalSubworkspace{T,TP,TW,TWn,TX,Tz0} <: eMCMC.LocalWorkspace{T}
    ll::Vector{Float64}
    ll_history::Vector{Vector{Float64}}
    P::TP
    WW::TW
    Wnr::TWn
    XX::TX
    z0s::Tz0

    function DiffusionLocalSubworkspace{T}(
            ll, ll_hist, P::TP, WW::TW, Wnr::TWn, XX::TX, z0s::Tz0
        ) where {T,TP,TW,TWn,TX,Tz0}
        new{T,TP,TW,TWn,TX,Tz0}(ll, ll_hist, P, WW, Wnr, XX, z0s)
    end
end

function DiffusionLocalSubworkspace{T}(
        sws::DiffusionGlobalSubworkspace,
        block_layout,
        num_mcmc_steps
    ) where T
    N, M = num_mcmc_steps, length(block_layout)
    DiffusionLocalSubworkspace{T}(
        Float64[-Inf for _=1:M],
        [zeros(Float64, M) for _=1:N],
        block_view(sws.P, block_layout),
        block_view(sws.WW, block_layout),
        block_wiener(sws.Wnr, block_layout),
        block_view(sws.XX, block_layout),
        block_start_pts(sws.XX, block_layout)
    )
end

function block_view(v, layout)
    [ view(v[bl_i], bl_i_range) for (bl_i, bl_i_range) in layout ]
end

function block_view(v, layout::Nothing)
    [ view(v[i], :) for i in 1:length(v) ]
end

function block_start_pts(XX, layout::Nothing)
    [ deepcopy(X[1].x[1]) for X in XX ]
end

function block_start_pts(XX, layout)
    [ deepcopy(XX[bl_i][bl_i_range[1]].x[1]) for (bl_i, bl_i_range) in layout ]
end

block_wiener(Wnr, layout::Nothing) = Wnr

block_wiener(Wnr, layout) = [ deepcopy(Wnr[bl_i]) for (bl_i, _) in layout ]

const _SYMS = Tuple{Vararg{Pair{Int64,Symbol},N} where N}
const _OBSIDX = Tuple{Vararg{Pair{Int64,Int64},N} where N}

struct LocalUpdtParamNames
    vp_names::Vector{Tuple{Vararg{Symbol,N} where N}}
    vp_names_aux::Vector{Vector{Tuple{Vararg{Symbol,N} where N}}}
    θ_local_names::Vector{_SYMS}
    θ_local_aux_names::Vector{Vector{_SYMS}}
    θ_local_obs::Vector{Vector{_OBSIDX}}

    function LocalUpdtParamNames(
            PPP::AbstractArray{<:AbstractArray{<:GuidProp}},
            data::AllObservations,
            θinit::Vector{Symbol},
        )
        psym, oind = OBS.local_symbols(data, PPP, x->x.P_target, θinit)
        pauxsym, _ = OBS.local_symbols(data, PPP, x->x.P_aux, θinit)

        new(
            [DD.var_parameter_names(PP) for PP in PPP],
            [[DD.var_parameter_names(P.P_aux) for P in PP] for PP in PPP],
            first.(psym), # target laws are the same on each obs of a single rec
            pauxsym,
            oind,
        )
    end

    LocalUpdtParamNames() = new()
end

"""
    struct DiffusionLocalWorkspace{T,SW,SWD,Tpr,AH} <: eMCMC.LocalWorkspace{T}
        sub_ws::SW
        sub_ws°::SW
        sub_ws_diff::SWD
        sub_ws_diff°::SWD
        xx0_priors::Tpr
        acceptance_history::Vector{AH}
        loc_pnames::LocalUpdtParamNames
        critical_param_change::Vector{Bool}
    end

Local workspace for MCMC problems concerning diffusions.

# Arguments
- `sub_ws`: `GenericLocalSubworkspace` from `ExtensibleMCMC` that acts as a
            local cache for some light parameter computations
- `sub_ws°`: same as `sub_ws`, but concerns the proposed parameter
- `sub_ws_diff`: `DiffusionLocalSubworkspace`, views to containers for the
                 accepted diffusion law/paths
- `sub_ws_diff°`: `DiffusionLocalSubworkspace`, views to containers for the
                  proposed diffusion law/paths
- `xx0_priors`: priors over starting points
- `acceptance_history`: history of results from the accept/reject
                        Metropolis–Hastings step
- `loc_pnames`: various lists of parameter names helping in performing the
                update in a clean way
- `critical_param_change`: a list of flags for whether a given update by itself
                           prompts for recomputation of the guiding term
"""
struct DiffusionLocalWorkspace{T,SW,SWD,Tpr,CGH,AH} <: eMCMC.LocalWorkspace{T}
    sub_ws::SW
    sub_ws°::SW
    sub_ws_diff::SWD
    sub_ws_diff°::SWD
    xx0_priors::Tpr
    acceptance_history::Vector{AH}
    loc_pnames::LocalUpdtParamNames
    cu_helper::CGH
    critical_param_change::Vector{Bool}
end

function create_workspace(
        backend::DiffusionMCMCBackend,
        mcmcupdate,
        global_ws,
        block_layout,
        num_mcmc_steps,
    )
    _state = (
        typeof(mcmcupdate) <: eMCMC.MCMCParamUpdate ?
        copy(state(global_ws, mcmcupdate)) :
        Float64[]
    )
    sub_ws = eMCMC.StandardLocalSubworkspace(_state, num_mcmc_steps)
    xx0_priors = block_setup_x0_priors(
        global_ws.data, block_layout, global_ws.sub_ws_diff
    )
    T = eltype(state(global_ws))
    sub_ws_diff = DiffusionLocalSubworkspace{T}(
        global_ws.sub_ws_diff, block_layout, num_mcmc_steps
    )
    init_ll!(sub_ws_diff)

    sub_ws_diff° = DiffusionLocalSubworkspace{T}(
        global_ws.sub_ws_diff°, block_layout, num_mcmc_steps
    )
    a_h = (
        typeof(mcmcupdate) <: eMCMC.MCMCParamUpdate ?
        fill(false, num_mcmc_steps) :
        [similar(v, Bool) for v in sub_ws_diff.ll_history]
    )

    init_update!(mcmcupdate, block_layout)

    lpn = (
        typeof(mcmcupdate) <: eMCMC.MCMCParamUpdate ?
        LocalUpdtParamNames(
            sub_ws_diff.P,
            global_ws.data,
            global_ws.pnames[eMCMC.coords(mcmcupdate)]
        ) :
        LocalUpdtParamNames()
    )

    cgu_helper = (
        typeof(mcmcupdate) <: DiffusionConjugGsnUpdate ?
        ConjugGsnHelper(
            mcmcupdate,
            global_ws,
            sub_ws_diff.P,
            lpn,
            #NOTE finally, computation type, stay with outofplace for now
        ) :
        ConjugGsnHelper()
    )

    crit_updt = (
        typeof(mcmcupdate) <: eMCMC.MCMCParamUpdate ?
        [
            GP.is_critical_update(
                PP, lpn.θ_local_aux_names[i], lpn.θ_local_obs[i]
            ) for (i,PP) in enumerate(sub_ws_diff.P)
        ] :
        [false for _ in eachindex(sub_ws_diff.P)]
    )

    DiffusionLocalWorkspace{
        T, typeof(sub_ws), typeof(sub_ws_diff), typeof(xx0_priors),
        typeof(cgu_helper), eltype(a_h)
    }(
        sub_ws, deepcopy(sub_ws), sub_ws_diff, sub_ws_diff°, xx0_priors, a_h,
        lpn, cgu_helper, crit_updt,
    )
end

function init_ll!(sub_ws_diff)
    XXX, PPP = sub_ws_diff.XX, sub_ws_diff.P
    for i in 1:length(sub_ws_diff.ll)
        recompute_guiding_term!(PPP[i])
        sub_ws_diff.ll[i] = loglikhd(PPP[i], XXX[i])
    end
end

function block_setup_x0_priors(data, block_layout, sub_ws_diff)
    #println(block_layout)
    map(block_layout) do b_l
        bl_i, bl_i_range = b_l
        #println(bl_i, bl_i_range)
        i1 = bl_i_range[1]
        pr = (
            i1 == 1 ?
            data.recordings[bl_i].x0_prior :
            KnownStartingPt(sub_ws_diff.XX[bl_i][i1].x[1])
        )
        pr
    end
end

function eMCMC.create_workspaces(backend::DiffusionMCMCBackend, mcmc::MCMC)
    map(1:length(mcmc.updates)) do i
        create_workspace(
            backend,
            mcmc.updates[i],
            mcmc.workspace,
            mcmc.schedule.extra_info.blocking_layout[i],
            mcmc.schedule.num_mcmc_steps
        )
    end
end

function eMCMC.update_workspaces!(
        local_updt::MCMCDiffusionImputation,
        global_ws::eMCMC.GlobalWorkspace,
        local_ws::eMCMC.LocalWorkspace,
        step,
        prev_ws,
    )
    for i in 1:length(step.same_layout)
        if !step.same_layout[i]
            change_laws_for_blocking!(local_ws, i) # TODO
            recompute_loglikhd!(local_ws, i)
        else
            prev_ws===nothing || ( ll(local_ws)[i] = ll(prev_ws)[i] )
        end
    end
    local_ws, local_updt
end

accepted(ws::DiffusionLocalWorkspace, i::Int) = ws.acceptance_history[i]
#set_accepted!(ws::DiffusionLocalWorkspace, i::Int, accepted, j::Int=1) = (ws.acceptance_history[i][j] = accepted)
ll(ws::DiffusionLocalWorkspace) = ws.sub_ws_diff.ll
ll°(ws::DiffusionLocalWorkspace) = ws.sub_ws_diff°.ll
ll(ws::DiffusionLocalWorkspace, i::Int) = ws.sub_ws_diff.ll_history[i]
ll°(ws::DiffusionLocalWorkspace, i::Int) = ws.sub_ws_diff°.ll_history[i]


text:
# Collection of functions to compute transit times, impact parameters, sky-velocity, and derivatives.

function detect_transits!(s::State{T},d::Derivatives{T},tt::TransitOutput{T},intr::Integrator{T}; grad::Bool=true) where T<:AbstractFloat
    rstar::T = 1e12 # Could this be removed?
    # Save current state as prior state
    set_state!(tt.s_prior, s)

    # Check to see if a transit may have occured before current state.
    # Sky is x-y plane; line of sight is z.
    # Body being transited is tt.ti, tt.occs is list of occultors:
    for i in tt.occs
        # Compute the relative sky velocity dotted with position:
        gi = g!(i,tt.ti,s.x,s.v)
        ri = sqrt(s.x[1,i]^2+s.x[2,i]^2+s.x[3,i]^2)  # orbital distance
        # See if sign of g switches, and if planet is in front of star (by a good amount):
        if gi > 0 && tt.gsave[i] < 0 && -s.x[3,i] > 0.25*ri && ri < rstar
            # A transit has occurred between the time steps - integrate ahl21!
            tt.count[i] += 1
            if tt.count[i] <= tt.ntt
                dt0 = -tt.gsave[i]*intr.h/(gi-tt.gsave[i]) # Starting estimate
                set_state!(s,tt.s_prior)
                findtransit!(tt.ti,i,dt0,s,d,tt,intr;grad=grad) # Search for transit time (integrating 'backward')
            end
        end
        tt.gsave[i] = gi
    end
    set_state!(s,tt.s_prior)
    return
end

function findtransit!(i::Int64,j::Int64,dt0::T,s::State{T},d::Derivatives{T},tt::TransitOutput{T},intr::Integrator;grad::Bool=true) where T<:AbstractFloat
    # Computes the transit time approximating the motion as a fraction of a AHL21 step backward in time.
    # Computes the impact parameter and sky-velocity, if passed a TransitParameters.
    # Also computes the Jacobian of the transit time with respect to the initial parameters, dtbvdq[1-3,7,n].
    # Initial guess using linear interpolation:

    s.dqdt .= 0.0
    set_state!(tt.s_transit,s)

    dt = one(T)
    iter = 0
    r3 = zero(T)
    gdot = zero(T)
    gsky = zero(T)
    stmp = zero(T)
    TRANSIT_TOL = 10*eps(dt)
    tt1 = dt0 + 1
    tt2 = dt0 + 2
    ITMAX = 20
    while abs(dt) > TRANSIT_TOL && iter < 20
        tt2 = tt1
        tt1 = dt0
        set_state!(s,tt.s_transit)
        # Advance planet state at start of step to estimated transit time:
        zero_out!(d)
        intr.scheme(s,d,dt0)
        # Compute time offset:
        gsky = g!(i,j,s.x,s.v)
        #  # Compute derivative of g with respect to time:
        gdot = gd!(i,j,s.x,s.v,s.dqdt)
        # Refine estimate of transit time with Newton's method:
        dt = -gsky/gdot
        # Add refinement to estimated time:
        #dt0 += dt
        dt0,stmp = comp_sum(dt0,stmp,dt)
        iter += 1
        # Break out if we have reached maximum iterations, or if
        # current transit time estimate equals one of the prior two steps:
        if (iter >= ITMAX) || (dt0 == tt1) || (dt0 == tt2)
            break
        end
    end

    if grad
        # Compute derivatives with updated time step
        set_state!(s,tt.s_transit)
        zero_out!(d)
        intr.scheme(s,d,dt0)
    end

    if tt isa TransitParameters
        # Compute the impact parameter and sky velocity, save to tt along with transit time.
        tt.ttbv[1,j,tt.count[j]] = s.t[1] + dt0

        if grad
            # Compute derivative of transit time, impact parameter, and sky velocity.
            vsky,bsky2 = dtbvdq!(i,j,s.x,s.v,s.jac_step,s.dqdt,tt.dtbvdq)
            tt.ttbv[2,j,tt.count[j]] = vsky
            tt.ttbv[3,j,tt.count[j]] = bsky2
            for itbv=1:3, k=1:7, p=1:s.n
                tt.dtbvdq0[itbv,j,tt.count[j],k,p] = tt.dtbvdq[itbv,k,p]
            end
            return
        end
        tt.ttbv[2,j,tt.count[j]] = calc_vsky(s.v,i,j)
        tt.ttbv[3,j,tt.count[j]] = calc_bsky2(s.x,i,j)
        return
    end

    tt.tt[j,tt.count[j]] = s.t[1] + dt0
    if grad
        # Compute derivative of transit time
        dtbvdq!(i,j,s.x,s.v,s.jac_step,s.dqdt,tt.dtdq)
        for k=1:7, p=1:s.n
            tt.dtdq0[j,tt.count[j],k,p] = tt.dtdq[1,k,p]
        end
    end
    return
end

function calc_dtdelements!(s::State{T},tt::TransitTiming{T}) where T <: AbstractFloat
    for i=1:s.n, j=1:tt.count[i]
        if j <= tt.ntt
            # Now, multiply by the initial Jacobian to convert time derivatives to orbital elements:
            for k=1:s.n, l=1:7
                tt.dtdelements[i,j,l,k] = zero(T)
                for p=1:s.n, q=1:7
                    tt.dtdelements[i,j,l,k] += tt.dtdq0[i,j,q,p]*s.jac_init[(p-1)*7+q,(k-1)*7+l]
                end
            end
        end
    end
end

function calc_dtdelements!(s::State{T},ttbv::TransitParameters{T}) where T <: AbstractFloat
    for itbv = 1:3, i=1:s.n, j = 1:ttbv.count[i]
        if j <= ttbv.ntt
            # Now, multiply by the initial Jacobian to convert time derivatives to orbital elements:
            for k=1:s.n, l=1:7
                ttbv.dtbvdelements[itbv,i,j,l,k] = zero(T)
                for p=1:s.n, q=1:7
                    ttbv.dtbvdelements[itbv,i,j,l,k] += ttbv.dtbvdq0[itbv,i,j,q,p]*s.jac_init[(p-1)*7+q,(k-1)*7+l]
                end
            end
        end
    end
end

"""Used in computing transit time inside `findtransit3`."""
function g!(i::Int64,j::Int64,x::Array{T,2},v::Array{T,2}) where {T <: Real}
    # See equation 8-10 Fabrycky (2008) in Seager Exoplanets book
    g = (x[1,j]-x[1,i])*(v[1,j]-v[1,i])+(x[2,j]-x[2,i])*(v[2,j]-v[2,i])
    return g
end

function gd!(i::Int64,j::Int64,x::Matrix{T},v::Matrix{T},dqdt::Array{T}) where T<:AbstractFloat
    return ((x[1,j]-x[1,i])*(dqdt[(j-1)*7+4]-dqdt[(i-1)*7+4])+(x[2,j]-x[2,i])*(dqdt[(j-1)*7+5]-dqdt[(i-1)*7+5])
            +(v[1,j]-v[1,i])*(dqdt[(j-1)*7+1]-dqdt[(i-1)*7+1])+(v[2,j]-v[2,i])*(dqdt[(j-1)*7+2]-dqdt[(i-1)*7+2]))
end

calc_bsky2(x::Matrix{T},i::Int64,j::Int64) where T<:AbstractFloat = (x[1,j]-x[1,i])^2 + (x[2,j]-x[2,i])^2
calc_vsky(v::Matrix{T},i::Int64,j::Int64) where T<:AbstractFloat = sqrt((v[1,j]-v[1,i])^2 + (v[2,j]-v[2,i])^2)

function dtbvdq!(i::Int64,j::Int64,x::Matrix{T},v::Matrix{T},jac_step::Matrix{T},dqdt::Vector{T},dtbvdq::Array{T}) where T<:AbstractFloat
    n = size(x)[2]
    # Compute time offset:
    gsky = g!(i,j,x,v)
    # Compute derivative of g with respect to time:
    gdot = gd!(i,j,x,v,dqdt)

    fill!(dtbvdq,zero(typeof(x[1])))
    indj = (j-1)*7+1
    indi = (i-1)*7+1
    for p=1:n
        indp = (p-1)*7
        for k=1:7
            # Compute derivatives:
            dtbvdq[1,k,p] = -((jac_step[indj  ,indp+k]-jac_step[indi  ,indp+k])*(v[1,j]-v[1,i])+(jac_step[indj+1,indp+k]-jac_step[indi+1,indp+k])*(v[2,j]-v[2,i])+
                          (jac_step[indj+3,indp+k]-jac_step[indi+3,indp+k])*(x[1,j]-x[1,i])+(jac_step[indj+4,indp+k]-jac_step[indi+4,indp+k])*(x[2,j]-x[2,i]))/gdot
        end
    end
    ntbv = size(dtbvdq)[1]
    if ntbv == 3
        # Compute the impact parameter and sky velocity:
        vsky = calc_vsky(v,i,j)
        bsky2 = calc_bsky2(x,i,j)
        # partial derivative v_{sky} with respect to time:
        dvdt = ((v[1,j]-v[1,i])*(dqdt[(j-1)*7+4]-dqdt[(i-1)*7+4])+(v[2,j]-v[2,i])*(dqdt[(j-1)*7+5]-dqdt[(i-1)*7+5]))/vsky
        # (note that \partial b/\partial t = 0 at mid-transit since g_{sky} = 0 mid-transit).
        for p=1:n
            indp = (p-1)*7
            for k=1:7
                # Compute derivatives:
                #v_{sky} = sqrt((v[1,j]-v[1,i])^2+(v[2,j]-v[2,i])^2)
                dtbvdq[2,k,p] = ((jac_step[indj+3,indp+k]-jac_step[indi+3,indp+k])*(v[1,j]-v[1,i])+(jac_step[indj+4,indp+k]-jac_step[indi+4,indp+k])*(v[2,j]-v[2,i]))/vsky + dvdt*dtbvdq[1,k,p]
                #b_{sky}^2 = (x[1,j]-x[1,i])^2+(x[2,j]-x[2,i])^2
                dtbvdq[3,k,p] = 2*((jac_step[indj  ,indp+k]-jac_step[indi  ,indp+k])*(x[1,j]-x[1,i])+(jac_step[indj+1,indp+k]-jac_step[indi+1,indp+k])*(x[2,j]-x[2,i]))
            end
        end
        return vsky,bsky2
    end
    return
end


text:
using ContextManagers: @with, SharedResource
using Julio
using Test

function test_channel()
    Julio.withtaskgroup() do tg
        send_endpoint, receive_endpoint = Julio.channel()
        @with(shared_send = SharedResource(send_endpoint)) do
            for i in 1:5
                Julio.spawn!(tg, shared_send) do send_endpoint
                    put!(send_endpoint, i)
                end
            end
        end
        @with(receive_endpoint) do
            @test sort!(collect(receive_endpoint)) == 1:5
        end
    end
end


text:
mutable struct RINEXData
    satellite_system::String
    prn::Int
    epoch::Int
    sv_clock_bias::Int
    sv_clock_drift::Int
    sv_clock_drift_rate::Int
    IODE::Float64
    C_rs::Float64
    Δn::Float64
    M_0::Float64
    C_uc::Float64
    e::Float64
    Cus::Float64
    sqrt_A::Float64
    t_oe::Float64
    C_ic::Float64
    ω::Float64
    Ω_dot::Float64
    IDOT::Float64
    l2codes::Float64
    gps_week::Float64
    l2pflag::Float64
    ura::Float64
    sv_health::Float64
    T_GD::Float64
    IODC::Float64
    ttm::Float64
    fit_interval::Float64
end


text:
using DataFrames
using Distributions
using PyPlot
using LinearAlgebra
function makeData()
    groupOne = rand(MvNormal([10.0, 10.0], 10.0 * Matrix{Float64}(I, 2,2)), 100)
    groupTwo = rand(MvNormal([0.0, 0.0], 10 * Matrix{Float64}(I, 2,2)), 100)
    groupThree = rand(MvNormal([15.0, 0.0], 10.0 * Matrix{Float64}(I, 2,2)), 100)
    return hcat(groupOne, groupTwo, groupThree)'
end

data = makeData()
scatter(data[1:100, 1], data[1:100, 2], color="blue")
scatter(data[101:200, 1], data[101:200, 2], color="red")
scatter(data[201:300, 1], data[201:300, 2], color="green")

using StatsBase

mutable struct kMeans
    x::DataFrames.DataFrame
    k::Int
end

function euclidean(sourcePoint, destPoint)
    sum = 0
    for i in 1:length(sourcePoint)
        sum += (destPoint[i] - sourcePoint[i]) ^ 2
    end
    dist = sqrt(sum)
    return dist
end

function minkowski(sourcePoint, destPoint)
    sum = 0
    for i in 1:length(sourcePoint)
        sum += abs(destPoint[i] - sourcePoint[i])
    end
    return sum
end

function calcDist(sourcePoint::Array, destPoint::Array; method="euclidean")

    if length(sourcePoint) != length(destPoint)
        error("The lengths of two arrays are different.")
        return
    end

    if method == "euclidean"
        return euclidean(sourcePoint, destPoint)
    elseif method == "minkowski"
        return minkowski(sourcePoint, destPoint)
    end
end

function classify(kMeans::kMeans)
    dataPointsNum = size(kMeans.x, 1)
    estimatedClass = Array{Int, dataPointsNum}
    direct_sample(1:kMeans.k, estimatedClass)

    while true
        # update representative points
        representativePoints = []
        for representativeIndex in 1:kMeans.k
            groupIndex = find(estimatedClass .== representativeIndex)
            groupData = kMeans.x[groupIndex, :]

            # TODO: check the return type of colwise
            representativePoint = [ valArray[1] for valArray in colwise(mean, groupData) ]
            push!(representativePoints, representativePoint)
        end

        # update group belonging
        tempEstimatedClass = Array{Int}(dataPointsNum)
        for dataIndex in 1:dataPointsNum
            dataPoint = Array(kMeans.x[dataIndex, :])
            distances = Array{Float64}(kMeans.k)
            for representativeIndex in 1:kMeans.k
                distances[representativeIndex] = calcDist(dataPoint, representativePoints[representativeIndex])
            end

            # TODO: check the existence of argmin
            tempEstimatedClass[dataIndex] = sortperm(distances)[1]
        end

        if estimatedClass == tempEstimatedClass
            break
        end
        estimatedClass = tempEstimatedClass
    end
    return estimatedClass
end

    kmeansElbow = kMeans(DataFrame(data), 3)

    predictedClassElbow = classify(kmeansElbow)


text:
using MAT

using MultivariateStats
using GLMakie
using ProperOrthogonalDecomposition
using Statistics
using Random
using Distributions

include("utils.jl")
include("decompositions.jl")

data = matread("../spod/jet_data/jetLES.mat")
p = data["p"]
x = data["x"][1,:]
r = data["r"][:,1]

iftest = true
ifica = false
ifpod = false
ifplotdata = false

n_latent = 5
n_out = 1
n_samples = 1000
s_noise = 0.1

Atest = rand(n_out,n_latent)
data = A*rand(Exponential(0.1),n_latent,n_samples) + randn(n_latent,n_samples)*s_noise
data = Matrix(data')

modes = fit(ICA, data, n_latent; maxiter=100000, tol=0.00001).W

# x = (1:length(x))*maximum(x)/length(x)
# r = (1:length(r))*maximum(r)/length(r)

# switches the spatial dimensions to the leading dimension
p = permutedims(p, [2, 3, 1])

# Strides in i (x) and j (r) and k (t)
ds_i = 1
ds_j = 1
ds_k = 1

# Get arrays
i_arr = 1 : ds_i : size(p, 1)
j_arr = 1 : ds_j : size(p, 2)
k_arr = 1 : ds_k : size(p, 3)

i_length = length(i_arr)
j_length = length(j_arr)
k_length = length(k_arr)

#
# p = p[:,:,k_arr]

# Remove mean
p = p .- mean(p[:])

# Normalize p
amp_p = mean(abs.(p))
p = p/amp_p

@show maximum(p)
@show minimum(p)

# Add inverted jet
size_fac = 1.0
num_subsamples = 2
# p = p[end:-1:1,end:-1:1,randperm(size(p,3))]*size_fac + p

# p = p[:,:,rand(randperm(size(p,3))[1:num_subsamples],size(p,3))] +
#     p[end:-1:1,end:-1:1,rand(randperm(size(p,3))[1:num_subsamples],size(p,3))]

p1 = repeat(p[:,:,1234],1,1,size(p,3))
p2 = repeat(p[end:-1:1,end:-1:1,2345],1,1,size(p,3))

p1[:,:,randperm(size(p,3))[1:1000]] .= 0.0
p2[:,:,randperm(size(p,3))[1:1000]] .= 0.0
p = p1 + p2

# p = p[end:-1:1,end:-1:1,:]*size_fac

# Try adding random component
# rand_mat = rand(size(p)...) * amp_p
# p = p + rand_mat

# Number of POD/ICA components to compute and/or plot
num_comp = 4

if ifplotdata
    println("Display Data")
    display_ts = 5000
    fig_data = heatmap(x[j_arr], r[i_arr], Matrix(p[i_arr, j_arr, display_ts]'),colormap =:summer)
    display(fig_data)
    save("figures/data.png",fig_data)
end

if ifica
    # The fastica algorithm used by the packages doesn't seem to converge
    # reliably. Convergence seems to depend on the initialization, so we try
    # it until it works
    println("ICA")

    # When to stop ICA if it never gets a valid one
    stop_ica_it = 5

    ica = nothing
    i_indices = nothing
    j_indices = nothing
    k = 1

    while isnothing(ica)
        # @code_warntype compute_ica(p, num_comp, ds_i, ds_j)
        @time try
            global ica, i_indices, j_indices = compute_ica(p, num_comp, ds_i, ds_j)
        catch
            @warn "ica not converged in $k-th attempt"
            global k = k + 1
        end
        if k > stop_ica_it break end
    end

    # This becomes a three-way tensor; the k-th component is given by components[:, :, k]
    # now included in the ica implementation
    # independent_components = reshape(ica, (length(i_indices), length(j_indices), size(ica, 2)))

    # Plots the independent components
    figs_per_row = 4
    fig_ica = Figure()
    for (k, ind) in enumerate(CartesianIndices((figs_per_row, num_comp))[:])
        if k <= num_comp
            heatmap(fig_ica[ind[2], ind[1]], x[j_arr], r[i_arr], Matrix(ica[i_arr, j_arr, k]'),colormap =:bluesreds)
        end
    end
    display(fig_ica)
    save("figures/ica.png",fig_ica)
end

if ifpod
    println("POD")
    pod, i_indices, j_indices = compute_pod_svd(p, ds_i, ds_j)
    # This becomes a three-way tensor; the k-th component is given by components[:, :, k]
    #proper_components = reshape(pod, (length(i_indices), length(j_indices), size(pod, 2)))
    pod = pod[:, :, 1:num_comp]

    # Plots the independent components
    figs_per_row = 4
    fig_pod = Figure()
    for (k, ind) in enumerate(CartesianIndices((figs_per_row, num_comp))[:])
        if k <= num_comp
            heatmap(fig_pod[ind[2], ind[1]], x[j_arr], r[i_arr], Matrix(pod[i_arr, j_arr, k]'),colormap =:bluesreds)
        end
    end
    display(fig_pod)
    save("figures/pod.png",fig_pod)
end

max_stars_repo_name: comp-physics/tensor-modal-decomp

max_issues_repo_name: comp-physics/tensor-modal-decomp

max_forks_repo_name: comp-physics/tensor-modal-decomp


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from Euler angles to Euler angle and axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angle_to_angleaxis.jl
# =============================================

# Functions: angle_to_angleaxis
# -----------------------------

@testset "Euler angles => Euler angle and axis" begin
    for T in (Float32, Float64)
        # We do not need comprehensive test here because `angle_to_angleaxis`
        # first converts the Euler angles to DCM and then to Euler angle and
        # axis. Those two operations are already heavily tested.

        ea = EulerAngles(deg2rad(T(45)), 0, 0, :XYZ)
        av = angle_to_angleaxis(ea)
        @test eltype(av) === T
        @test av.a ≈ deg2rad(T(45))
        @test av.v ≈ T[1, 0, 0]

        ea = EulerAngles(deg2rad(T(45)), 0, 0, :YXZ)
        av = angle_to_angleaxis(ea)
        @test eltype(av) === T
        @test av.a ≈ deg2rad(T(45))
        @test av.v ≈ T[0, 1, 0]

        ea = EulerAngles(deg2rad(T(45)), 0, 0, :ZXY)
        av = angle_to_angleaxis(ea)
        @test eltype(av) === T
        @test av.a ≈ deg2rad(T(45))
        @test av.v ≈ T[0, 0, 1]

        ea = EulerAngles(0, deg2rad(T(45)), 0, :ZXY)
        av = angle_to_angleaxis(ea)
        @test eltype(av) === T
        @test av.a ≈ deg2rad(T(45))
        @test av.v ≈ T[1, 0, 0]
    end
end


text:
###################
### Point Types ###
###################

# See the REQUIRE file for details.
const FieldVector{T} = StaticArrays.FieldVector{3, T}

"""
    LLA(lat, lon, alt = 0.0)
    LLA(lat = ϕ, lon = Θ, alt = h)

Latitude, longitude, and alititude co-ordinates. *Note:* assumes degrees not radians
"""
struct LLA{T <: Number}
    lat::T
    lon::T
    alt::T
end
LLA(lat, lon, alt=zero(lat)) = LLA(promote(lat, lon, alt)...)
LLA(;lat=NaN,lon=NaN,alt=0.0) = LLA(lat,lon,alt) # Constructor that is idependent of storage order

Base.show(io::IO, lla::LLA) = print(io, "LLA(lat=$(lla.lat)°, lon=$(lla.lon)°, alt=$(lla.alt))") # Maybe show to nearest mm, or 6 decimel places?
Base.isapprox(lla1::LLA, lla2::LLA; atol = 1e-6, kwargs...) = isapprox(lla1.lat, lla2.lat; atol = 180*atol/6.371e6, kwargs...) & isapprox(lla1.lon, lla2.lon; atol = 180*atol/6.371e6, kwargs...) & isapprox(lla1.alt, lla2.alt; atol = atol, kwargs...) # atol in metres (1μm)

# For radians, we could export ° = pi/180, then output in degrees... ??

"""
    LatLon(lat, lon)
    LatLon(lat = ϕ, lon = Θ)

Latitude and longitude co-ordinates. *Note:* assumes degrees not radians
"""
struct LatLon{T <: Number}
    lat::T
    lon::T
end
LatLon(lat::Number, lon::Number) = LatLon(promote(lat, lon)...)
LatLon(;lat=NaN,lon=NaN) = LatLon(lat,lon) # Constructor that is independent of storage order
LatLon(lla::LLA) = LatLon(lla.lat, lla.lon)
function LatLon(x, datum::Union{Datum,Ellipsoid})
    lla = LLA(x, datum)
    return LatLon(lla.lat, lla.lon)
end
LLA(ll::LatLon) = LLA(ll.lat, ll.lon)
Base.show(io::IO, ll::LatLon) = print(io, "LatLon(lat=$(ll.lat)°, lon=$(ll.lon)°)") # Maybe show to nearest mm, or 6 decimel places?
Base.isapprox(ll1::LatLon, ll2::LatLon; atol = 1e-6, kwargs...) = isapprox(ll1.lat, ll2.lat; atol = 180*atol/6.371e6, kwargs...) & isapprox(ll1.lon, ll2.lon; atol = 180*atol/6.371e6, kwargs...) # atol in metres (1μm)

"""
    ECEF(x, y, z)

Earth-Centered-Earth-Fixed (ECEF) coordinates. A global Cartesian coordinate
system rotating with the Earth.
"""
struct ECEF{T <: Number} <: FieldVector{T}
    x::T
    y::T
    z::T
end
@inline function ECEF(x,y,z)
    T = promote_type(promote_type(typeof(x),typeof(y)), typeof(z))
    ECEF{T}(x,y,z)
end
Base.show(io::IO, ::MIME"text/plain", ecef::ECEF) = print(io, "ECEF($(ecef.x), $(ecef.y), $(ecef.z))")

"""
    ENU(e, n, u = 0.0)

East-North-Up (ENU) coordinates. A local Cartesian coordinate system, linearized about a reference point.
"""
struct ENU{T <: Number} <: FieldVector{T}
    e::T
    n::T
    u::T
end
ENU(x :: T, y :: T) where {T} = ENU(x, y, zero(T))
@inline function ENU(x,y,z)
    T = promote_type(promote_type(typeof(x),typeof(y)), typeof(z))
    ENU{T}(x,y,z)
end
Base.show(io::IO, ::MIME"text/plain", enu::ENU) = print(io, "ENU($(enu.e), $(enu.n), $(enu.u))")

"""
    NED(n, e, d = 0.0)

North-East-Down (NED) coordinates. A local Cartesian coordinate system, linearized about a reference point.
"""
struct NED{T <: Number} <: FieldVector{T}
    n::T
    e::T
    d::T
end
NED(x :: T, y :: T) where {T} = NED(x, y, zero(T))
@inline function NED(x,y,z)
    T = promote_type(promote_type(typeof(x),typeof(y)), typeof(z))
    NED{T}(x,y,z)
end
Base.show(io::IO, ::MIME"text/plain", ned::NED) = print(io, "NED($(ned.n), $(ned.e), $(ned.d))")

"""
    UTM(x, y, z = 0.0)

Universal transverse Mercator (UTM) coordinates. Common projection type for
world points. Zone not included in coordinates - it is a parameterized in the
relavant transformations `UTMfromLLA` and `LLAfromUTM` (see also the `UTMZ` type).

This type may be used to parameterize UPS coordinates (Universal Polar
Stereographic) to accurately represent the polar regions, in zone "0".
"""
struct UTM{T <: Number}
    x::T
    y::T
    z::T
end
UTM(x :: T, y :: T) where {T} = UTM(x, y, zero(T))
Base.isapprox(utm1::UTM, utm2::UTM; atol = 1e-6, kwargs...) = isapprox(utm1.x, utm2.x; atol = atol, kwargs...) & isapprox(utm1.y, utm2.y; atol = atol, kwargs...) & isapprox(utm1.z, utm2.z; atol = atol, kwargs...) # atol in metres (1μm)
Base.show(io::IO, utm::UTM) = print(io, "UTM($(utm.x), $(utm.y), $(utm.z))")

"""
    UTMZ(x, y, z = 0.0, zone::Integer, hemisphere::Bool)

Universal transverse Mercator (UTM) coordinates with zone number. Common
projection type for world points. The UTM zone is included in coordinates
(see also the `UTM` type).

This type may be used to parameterize UPS coordinates (Universal Polar
Stereographic) to accurately represent the polar regions, in zone "0".
"""
struct UTMZ{T <: Number}
    x::T
    y::T
    z::T
    zone::UInt8
    isnorth::Bool # which hemisphere
end
UTMZ(x::T, y::T, zone::Integer, isnorth::Bool) where {T} = UTMZ(x, y, zero(T), UInt8(zone), isnorth)
UTMZ(x, y, z, zone::Integer, isnorth::Bool) = UTMZ(promote(x, y, z)..., UInt8(zone), isnorth)
UTMZ(utm::UTM, zone::Integer, isnorth::Bool) = UTMZ(utm.x, utm.y, utm.z, UInt8(zone), isnorth)
UTM(utmz::UTMZ) = UTM(utmz.x, utmz.y, utmz.z)
Base.isapprox(utm1::UTMZ, utm2::UTMZ; atol = 1e-6, kwargs...) = (utm1.zone == utm2.zone) & (utm1.isnorth == utm2.isnorth) & isapprox(utm1.x, utm2.x; atol = atol, kwargs...) & isapprox(utm1.y, utm2.y; atol = atol, kwargs...) & isapprox(utm1.z, utm2.z; atol = atol, kwargs...) # atol in metres (1μm)
Base.show(io::IO, utm::UTMZ) = print(io, "UTMZ($(utm.x), $(utm.y), $(utm.z), zone=$(utm.zone == 0 ? "polar" : Int(utm.zone)) ($(utm.isnorth ? "north" : "south")))")


text:
function default_builtins()
    a = string()
    b = false
    c = zero(Int)
    @assert(a == "")
    @assert(b == false)
    @assert(c == 0)
end

if abspath(PROGRAM_FILE) == @__FILE__
    default_builtins()
    a = max(1, 2)
    @assert(a == 2)
    b = min(1, 2)
    @assert(b == 1)
    println("OK")
end


text:
"""

```
mddf_linkedcells(trajectory::Trajectory, options::Options,samples::Samples, RNG, mddf_compute!)
```

Computes the MDDF using linked cells, serial version.

"""
function mddf_linkedcells(trajectory::Trajectory, options::Options,
                          samples::Samples, RNG, mddf_compute!)

  # Initializing the structure that carries all results
  R = Result(trajectory,options)

  # Data structure to be passed to mddf_frame
  framedata = FrameData(trajectory,R)

  # Print some information about this run
  options.silent || title(R,trajectory.solute,trajectory.solvent)

  # Computing all minimum-distances
  if ! options.silent
    progress = Progress(R.nframes_read,1)
  end
  for iframe in 1:R.lastframe_read

    # reading coordinates of next frame
    nextframe!(trajectory)
    if iframe < options.firstframe
      continue
    end
    if iframe%options.stride != 0
      continue
    end
    mddf_compute!(iframe,framedata,options,RNG,R)

    options.silent || next!(progress)
  end # frames
  closetraj(trajectory)

  # Setup the final data structure with final values averaged over the number of frames,
  # sampling, etc, and computes final distributions and integrals
  finalresults!(R,options,trajectory,samples)
  options.silent || println(bars)

  return R

end


text:
module SimpleSymbolic

export S

"""
    S(x::Symbol)
    S(x::Expr)
    S(x::Number)

Create a symbolic wrapper `S <: Number` for the object `x`.
"""
struct S <: Number
    x
end

Base.show(io::IO, s::S) = print(io, "S($(s.x))")

# Promotion & conversion wraps types in `S`
Base.convert(::Type{S}, x::S) = x
Base.convert(::Type{S}, x) = S(x)
Base.promote_rule(::Type{S}, ::Type{<:Any}) = S

# Additive & multiplicative identity
struct Zero ; end
struct One  ; end

Base.zero(::Type{S}) = S(Zero())
Base.one(::Type{S}) = S(One())
Base.zero(::S) = zero(S)
Base.one(::S) = one(S)

# Arithmetic on `S` is implemented by unwrapping and rewrapping so that we can
# use the standard promotion machinery but also do some basic simplifications
# on the fly.
_funcname(f) = typeof(f).name.mt.name # See Base show(::IO, ::Function)
_op(f, es...) = Expr(:call, _funcname(f), es...)

unwrap() = ()
unwrap(s::S, ss::S...) = (s.x, unwrap(ss...)...)

Base.:+(s1::S, ss::S...) = S(_op(+, unwrap(s1, ss...)...))
Base.:*(s1::S, ss::S...) = S(_op(*, unwrap(s1, ss...)...))
Base.:-(s1::S, ss::S...) = S(_op(-, unwrap(s1, ss...)...))
Base.:/(s1::S, ss::S...) = S(_op(/, unwrap(s1, ss...)...))

for func in [:sqrt, :exp, :log, :sin, :cos, :conj]
    @eval Base.$func(s::S) = S(_op($func, s.x))
end

# Basic simplification
# Simplification for numbers
_op(f, xs::Number...) = f(xs...)

# Simplification for standard left-fold reductions
iscall(x::Expr, name) = x.head == :call && x.args[1] == name

_op(f::typeof(+), x::Zero, x2) = x2
_op(f::typeof(+), x::Zero, xs...) = _op(f, xs...)

_op(f::typeof(*), x::Zero, xs...) = Zero()
_op(f::typeof(*), x::One, x2) = x2
_op(f::typeof(*), x::One, xs...) = _op(f, xs...)

_op(f::typeof(conj), x::Expr) = iscall(x, :conj) ? x.args[2] : Expr(:call, _funcname(f), x)

end # module


text:
module TestSerd
using Test
using Serd, Serd.CSerd, Serd.RDF
import Serd: to_serd, from_serd

import ..TurtleEx1

# Reader
########

@test read_rdf_string(TurtleEx1.turtle) == TurtleEx1.statements
@test read_rdf_file(TurtleEx1.turtle_path) == TurtleEx1.statements

# Writer
########

@test strip(sprint(write_rdf, TurtleEx1.statements)) == strip(TurtleEx1.turtle_alt)

# Data types
############

# Node conversion
@test to_serd(ResourceCURIE("rdf", "type")) == SerdNode("rdf:type", SERD_CURIE)
@test from_serd(SerdNode("rdf:type", SERD_CURIE)) == ResourceCURIE("rdf", "type")

@test to_serd(ResourceURI("foo")) == SerdNode("foo", SERD_URI)
@test from_serd(SerdNode("foo", SERD_URI)) == ResourceURI("foo")

@test to_serd(Blank("b1")) == SerdNode("b1", SERD_BLANK)
@test from_serd(SerdNode("b1", SERD_BLANK)) == Blank("b1")

# Statement conversion
triple = Triple(Resource("bob"), Resource("rdf","type"), Resource("Person"))
stmt = SerdStatement(
  0,
  nothing,
  SerdNode("bob", SERD_URI),
  SerdNode("rdf:type", SERD_CURIE),
  SerdNode("Person", SERD_URI),
  nothing,
  nothing,
)
@test to_serd(triple) == stmt
@test from_serd(stmt) == triple

triple = Triple(Resource("bob"), Resource("name"), Literal("Bob"))
stmt = SerdStatement(
  0,
  nothing,
  SerdNode("bob", SERD_URI),
  SerdNode("name", SERD_URI),
  SerdNode("Bob", SERD_LITERAL),
  nothing,
  nothing,
)
@test to_serd(triple) == stmt
@test from_serd(stmt) == triple

triple = Triple(Resource("bob"), Resource("age"), Literal(50))
stmt = SerdStatement(
 0,
 nothing,
 SerdNode("bob", SERD_URI),
 SerdNode("age", SERD_URI),
 SerdNode("50", SERD_LITERAL),
 SerdNode(Serd.XSD_INT, SERD_URI),
 nothing,
)
@test to_serd(triple) == stmt
@test from_serd(stmt) == triple

quad = Quad(Resource("bob"), Resource("friendly"), Literal(true), Resource("people"))
stmt = SerdStatement(
 0,
 SerdNode("people", SERD_URI),
 SerdNode("bob", SERD_URI),
 SerdNode("friendly", SERD_URI),
 SerdNode("true", SERD_LITERAL),
 SerdNode(Serd.XSD_BOOLEAN, SERD_URI),
 nothing,
)
@test to_serd(quad) == stmt
@test from_serd(stmt) == quad

end


text:
# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/action_type_support_struct.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/message_initialization.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/message_type_support_struct.h
# Automatically generated using Clang.jl wrap_c

function get_message_typesupport_handle(handle, identifier)
    ccall((:get_message_typesupport_handle, librosidl_generator_c), Ptr{rosidl_message_type_support_t}, (Ptr{rosidl_message_type_support_t}, Cstring), handle, identifier)
end

function get_message_typesupport_handle_function(handle, identifier)
    ccall((:get_message_typesupport_handle_function, librosidl_generator_c), Ptr{rosidl_message_type_support_t}, (Ptr{rosidl_message_type_support_t}, Cstring), handle, identifier)
end
# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/primitives_sequence.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/primitives_sequence_functions.h
# Automatically generated using Clang.jl wrap_c

function rosidl_generator_c__bool__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__bool__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__bool__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__bool__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__bool__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__bool__Sequence},), sequence)
end

function rosidl_generator_c__byte__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__byte__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__byte__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__byte__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__byte__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__byte__Sequence},), sequence)
end

function rosidl_generator_c__char__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__char__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__char__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__char__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__char__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__char__Sequence},), sequence)
end

function rosidl_generator_c__float32__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__float32__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__float32__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__float32__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__float32__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__float32__Sequence},), sequence)
end

function rosidl_generator_c__float64__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__float64__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__float64__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__float64__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__float64__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__float64__Sequence},), sequence)
end

function rosidl_generator_c__int8__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__int8__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__int8__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__int8__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__int8__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__int8__Sequence},), sequence)
end

function rosidl_generator_c__uint8__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__uint8__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__uint8__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__uint8__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__uint8__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__uint8__Sequence},), sequence)
end

function rosidl_generator_c__int16__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__int16__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__int16__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__int16__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__int16__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__int16__Sequence},), sequence)
end

function rosidl_generator_c__uint16__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__uint16__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__uint16__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__uint16__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__uint16__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__uint16__Sequence},), sequence)
end

function rosidl_generator_c__int32__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__int32__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__int32__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__int32__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__int32__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__int32__Sequence},), sequence)
end

function rosidl_generator_c__uint32__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__uint32__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__uint32__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__uint32__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__uint32__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__uint32__Sequence},), sequence)
end

function rosidl_generator_c__int64__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__int64__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__int64__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__int64__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__int64__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__int64__Sequence},), sequence)
end

function rosidl_generator_c__uint64__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__uint64__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__uint64__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__uint64__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__uint64__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__uint64__Sequence},), sequence)
end
# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/service_type_support_struct.h
# Automatically generated using Clang.jl wrap_c

function get_service_typesupport_handle(handle, identifier)
    ccall((:get_service_typesupport_handle, librosidl_generator_c), Ptr{rosidl_service_type_support_t}, (Ptr{rosidl_service_type_support_t}, Cstring), handle, identifier)
end

function get_service_typesupport_handle_function(handle, identifier)
    ccall((:get_service_typesupport_handle_function, librosidl_generator_c), Ptr{rosidl_service_type_support_t}, (Ptr{rosidl_service_type_support_t}, Cstring), handle, identifier)
end
# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/string.h
# Automatically generated using Clang.jl wrap_c

# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/string_functions.h
# Automatically generated using Clang.jl wrap_c

function rosidl_generator_c__String__init(str)
    ccall((:rosidl_generator_c__String__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__String},), str)
end

function rosidl_generator_c__String__fini(str)
    ccall((:rosidl_generator_c__String__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__String},), str)
end

function rosidl_generator_c__String__assignn(str, value, n)
    ccall((:rosidl_generator_c__String__assignn, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__String}, Cstring, Csize_t), str, value, n)
end

function rosidl_generator_c__String__assign(str, value)
    ccall((:rosidl_generator_c__String__assign, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__String}, Cstring), str, value)
end

function rosidl_generator_c__String__Sequence__init(sequence, size)
    ccall((:rosidl_generator_c__String__Sequence__init, librosidl_generator_c), Bool, (Ptr{rosidl_generator_c__String__Sequence}, Csize_t), sequence, size)
end

function rosidl_generator_c__String__Sequence__fini(sequence)
    ccall((:rosidl_generator_c__String__Sequence__fini, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__String__Sequence},), sequence)
end

function rosidl_generator_c__String__Sequence__create(size)
    ccall((:rosidl_generator_c__String__Sequence__create, librosidl_generator_c), Ptr{rosidl_generator_c__String__Sequence}, (Csize_t,), size)
end

function rosidl_generator_c__String__Sequence__destroy(sequence)
    ccall((:rosidl_generator_c__String__Sequence__destroy, librosidl_generator_c), Cvoid, (Ptr{rosidl_generator_c__String__Sequence},), sequence)
end
# Julia wrapper for header: /opt/ros/crystal/include/rosidl_generator_c/visibility_control.h
# Automatically generated using Clang.jl wrap_c


text:
idx = [4,5,8,10]
one_to_n = collect(1:10)
trnidx = [2,4,8,10]
validx = [1,3,5,6,7,9]

# found=findall((in)(idx),trnidx)
# integersNew=trnidx[found]
# trnidxNEW=findall((in)(integersNew),idx)
found = findall((in)(idx), trnidx)
found = trnidx[found]
trnidxNEW = findall((in)(found), idx)

found = findall((in)(idx), validx)
integersNew = validx[found]
validxNEW = findall((in)(integersNew), idx)

# function Base.getindex(r::DTMTable,idx::Vector{T};compact_features=false) where T <: Number
	@assert !compact_features
	key = r.key[idx]
numerator = r.numerator[idx]
denominator = r.denominator[idx]
weight = r.weight[idx]
	features = r.features[idx,:]

	# trn and val idx:
	one_to_n = collect(1:length(r.weight))
	rows_to_keep = idx # one_to_n[idx]
	foundtrn = findall((in)(rows_to_keep), r.trnidx)
	foundval = findall((in)(rows_to_keep), r.validx)
	new_trnidx = r.trnidx[foundtrn]
new_validx = r.validx[foundval]

	(length(new_validx) != 0) || @warn("DTM: Validation index of the DTMTable has length zero! You may need to redefine validx. The function resample_trnvalidx! might be helpful.")
	(length(new_trnidx) != 0) || @warn("DTM: Training index of the DTMTable has length zero! You may need to redefine trnidx. The function resample_trnvalidx! might be helpful.")

#
#	dt=DTMTable(key,new_trnidx,new_validx,numerator,denominator,weight,features,deepcopy(r.candMatWOMaxValues),deepcopy(r.mappings))
#   return dt
# end