text:
# Define abstract interface
"""
$(TYPEDEF)

Abstract supertype for tracer distributions nodes.
"""
abstract type TracerDistributions{T} end

"""
$(SIGNATURES)

Return the tracer concentration [number of particles/m³] for a given location
r [m] and time t [s].
"""
function tracerConcentration(tracerDistribution::TracerDistributions{T},r::SVector{3,T},t::T) where {T}
    error("tracerConcentration is not implemented for tracerDistribution::$(typeof(tracerDistribution)).")
end

"""
$(SIGNATURES)

Return if (`Bool`) magnetization field is zero at a given location r [m] at all
times.
"""
function isZero(tracerDistribution::TracerDistributions{T},r::SVector{3,T}) where {T}
    error("isZero is not implemented for tracerDistribution::$(typeof(tracerDistribution)).")
end

"""
$(SIGNATURES)

Return `QuadratureNodes` required to integrate over a region containing the
support of the `tracerDistribution`
"""
function QuadratureNodes(tracerDistribution::TracerDistributions{T},density) where {T}
    error("QuadratureNodes is not implemented for tracerDistribution::$(typeof(tracerDistribution)).")
end

"""
$(TYPEDEF)

Abstract supertype for static tracer distributions nodes.
"""
abstract type StaticTracerDistributions{T} <: TracerDistributions{T} end

function tracerConcentration(tracerDistribution::StaticTracerDistributions{T},position::SVector{3,T},t::T) where {T}
    if isZero(tracerDistribution,position)
        return zero(T)
    else
        return tracerDistribution.tracerConcentration
    end
end

# Define specific subtypes
## homogeneous tracer distribution within axis aligned box
struct AxisOrientedBox{T} <: StaticTracerDistributions{T}
    tracerConcentration::T
    center::SVector{3,T}
    positiveHalfLengths::SVector{3,T}

    AxisOrientedBox(tracerConcentration::T,center::Vector{T},positiveHalfLengths::Vector{T}) where {T} = new{T}(tracerConcentration,SVector{3,T}(center),SVector{3,T}(positiveHalfLengths))
end

@inline function isZero(tracerDistribution::AxisOrientedBox{T},r::SVector{3,T}) where {T}
    return !all(abs.(tracerDistribution.center-r).<=tracerDistribution.positiveHalfLengths)
end

function QuadratureNodes(tracerDistribution::AxisOrientedBox{T},density=10/1e-3) where {T}
    nx,ny,nz = round.(Int,2*density*tracerDistribution.positiveHalfLengths)
    nx = max(nx,1)
    ny = max(ny,1)
    nz = max(nz,1)
    n = tuple(nx,ny,nz)

    nodesx, weightsx = gausslegendre(nx)
    nodesx = tracerDistribution.positiveHalfLengths[1]*nodesx .+ tracerDistribution.center[1]
    weightsx *= tracerDistribution.positiveHalfLengths[1]
    nodesy, weightsy = gausslegendre(ny)
    nodesy = tracerDistribution.positiveHalfLengths[2]*nodesy .+ tracerDistribution.center[2]
    weightsy *= tracerDistribution.positiveHalfLengths[2]
    nodesz, weightsz = gausslegendre(nz)
    nodesz = tracerDistribution.positiveHalfLengths[3]*nodesz .+ tracerDistribution.center[3]
    weightsz *= tracerDistribution.positiveHalfLengths[3]

    return GaussLegendre{T}(nodesx,weightsx,nodesy,weightsy,nodesz,weightsz,n)
end

## homogeneous tracer distribution within sphere
struct Sphere{T} <: StaticTracerDistributions{T}
    tracerConcentration::T
    center::SVector{3,T}
    radius::T

    Sphere(tracerConcentration::T,center::Vector{T},radius::T) where {T} = new{T}(tracerConcentration,SVector{3,T}(center),radius)
end

@inline function isZero(tracerDistribution::Sphere{T},r::SVector{3,T}) where {T}
    return norm(tracerDistribution.center-r)>tracerDistribution.radius
end

function QuadratureNodes(tracerDistribution::Sphere{T},density=10/1e-3) where {T}
    r = tracerDistribution.radius
    center = tracerDistribution.center
    lfb = center - SVector{3,T}(r,r,r)
    rrt = center + SVector{3,T}(r,r,r)
    return MidPoint(lfb,rrt,density)
end

# misc
"""
    `numberOfParticles(cmolFe=1e-3,d=25e-9)`

Most often, in experiments the iron density (mol(Fe)/L) is provided rather
than the particle concentration (number of particles per m³). This function
allows to convert the former `cmolFe` (mol(Fe)/L) to number of particles
per m³ assuming nanoparticles with a spherical magnetite core of diameter `d`.
"""
function numberOfParticles(cmolFe=1e-3,d=25e-9)
    molarMassFe3O4 = 0.231533 # kg/mol(Fe3O4)
    densityFe3O4 = 5200 # kg/m³
    volumeNanoparticle = pi/6 * d^3 # m³
    weightNanoparticle = volumeNanoparticle .* densityFe3O4 # kg
    conversionFactor = molarMassFe3O4 / (3 * weightNanoparticle) # Number of nanoparticles per mol(Fe)/m³ (factor 3 accounts for the fact that Fe3O4 contains 3 iron atoms)
    return cmolFe * conversionFactor # Number of nanoparticles per m³
end


text:
mutable struct LinearLeastSquares{BFVT}

    # Basis functions linear with respect to quantities to be estimated
    # h[i] ∀ i ∈ [0, n] should be a function which takes the independant
    # variable as an argument (i.e h[i](tₖ) where tₖ is the time at
    # the kth epoch.
    h::BFVT

    # Basis function matrix
    H::Matrix{Float64}

    # Estimate vector
    xhat::Vector{Float64}

    # Measurements
    y::Vector{Float64}

    # Independant variables
    t::Vector{Float64}

end

function LinearLeastSquares(h::AbstractVector)
    LinearLeastSquares(h, Matrix{Float64}(undef, (0, 0)), Vector{Float64}(undef, length(h)),
        Vector{Float64}(undef, 0), Vector{Float64}(undef, 0))
end

function FeedMeasurementBatch!(lls::LinearLeastSquares, ys::AbstractVector, ts::AbstractVector)
    m = length(ys)
    n = length(lls.h)
    if m != length(ts)
        throw(ArgumentError("Measurement vector and indepandant variable vector are not the same length."))
    end

    # Set measurements and indepandant variable vectors
    lls.y = ys
    lls.t = ts

    # Size and fill H matrix
    H = zeros(m, n)
    @inbounds for row in 1:m
        for col in 1:n
            H[row, col] = lls.h[col](lls.t[row])
        end
    end
    lls.H = H

    return nothing
end

function ComputeEstimate!(lls::LinearLeastSquares)
    # xhat = (HᵀH)⁻¹Hᵀy

    # Compute HᵀH and factorize result
    HTH = factorize(transpose(lls.H)*lls.H)

    # Compute estimate
    ldiv!(lls.xhat, HTH, transpose(lls.H)*lls.y)

    return nothing
end

function GetEstimate(lls::LinearLeastSquares)
    return lls.xhat
end

function GetMeasurementEstimate(lls::LinearLeastSquares, t::Real)
    y = 0.0
    @inbounds for i in 1:length(lls.h)
        y += lls.h[i](t)*lls.xhat[i]
    end
    return y
end

ext: jl


text:
vis = Visualizer()

if haskey(ENV, "CI")
    port = split(split(url(vis), ':')[end], '/')[1]
    @show port
    stream, proc = open(`julia $(joinpath(@__DIR__, "dummy_websocket_client.jl")) $port`)
else
    proc = nothing
    open(vis)
end

wait(vis)
delete!(vis)

@testset "self-contained visualizer" begin
    @testset "shapes" begin
        v = vis[:shapes]
        delete!(v)
        settransform!(v, Translation(1., 0, 0))
        @testset "box" begin
            setobject!(v[:box], HyperRectangle(Vec(0., 0, 0), Vec(0.1, 0.2, 0.3)))
            settransform!(v[:box], Translation(-0.05, -0.1, 0))
        end

        @testset "cylinder" begin
            setobject!(v[:cylinder], Mesh(
               Cylinder(Point(0., 0, 0), Point(0, 0, 0.2), 0.1),
               MeshLambertMaterial(color=colorant"lime")))
            settransform!(v[:cylinder], Translation(0, 0.5, 0.0))
        end

        @testset "sphere" begin
            setobject!(v[:sphere],
                HyperSphere(Point(0., 0, 0), 0.15),
                MeshLambertMaterial(color=colorant"maroon"))
            settransform!(v[:sphere], Translation(0, 1, 0.15))
        end

        @testset "ellipsoid" begin
            setobject!(v[:ellipsoid], HyperEllipsoid(Point(0., 1.5, 0), Vec(0.3, 0.1, 0.1)))
            settransform!(v[:ellipsoid], Translation(0, 0, 0.1))
        end

        @testset "cube" begin
            setobject!(v[:cube], HyperCube(Vec(-0.1, -0.1, 0), 0.2), MeshBasicMaterial())
            settransform!(v[:cube], Translation(0, 2.0, 0))
        end

        @testset "more complicated cylinder" begin
            setobject!(v[:cylinder2], Cylinder(Point(0, 2.5, 0), Point(0.1, 0.1, 0), 0.05))
            settransform!(v[:cylinder2], Translation(0, 0, 0.05))
        end
    end

    @testset "meshes" begin
        v = vis[:meshes]
        @testset "cat" begin
            mesh = load(joinpath(Pkg.dir("GeometryTypes"), "test", "data", "cat.obj"))
            setobject!(v[:cat], mesh)
            settransform!(v[:cat], Translation(0, -1, 0) ∘ LinearMap(RotZ(π)) ∘ LinearMap(RotX(π/2)))
        end

        @testset "textured valkyrie" begin
            head = Mesh(
                load(joinpath(MeshCat.VIEWER_ROOT, "data", "head_multisense.obj"), GLUVMesh),
                MeshLambertMaterial(
                    map=Texture(
                        image=PngImage(joinpath(MeshCat.VIEWER_ROOT, "data", "HeadTextureMultisense.png"))
                    )
                ))
            setobject!(v[:valkyrie, :head], head)
            settransform!(v[:valkyrie, :head], Translation(0, 0.5, 0.5))
        end
    end

    @testset "points" begin
        v = vis[:points]
        settransform!(v, Translation(-1, 0, 0))
        @testset "random points" begin
            verts = rand(Point3f0, 100000);
            colors = reinterpret(RGB{Float32}, verts);
            setobject!(v[:random], PointCloud(verts, colors))
            settransform!(v[:random], Translation(-0.5, -0.5, 0))
        end
    end

    @testset "Polyhedra" begin
        ext1 = SimpleVRepresentation([0 1 2 3;0 2 1 3; 1 0 2 3; 1 2 0 3; 2 0 1 3; 2 1 0 3;
                                    0 1 3 2;0 3 1 2; 1 0 3 2; 1 3 0 2; 3 0 1 2; 3 1 0 2;
                                    0 3 2 1;0 2 3 1; 3 0 2 1; 3 2 0 1; 2 0 3 1; 2 3 0 1;
                                    3 1 2 0;3 2 1 0; 1 3 2 0; 1 2 3 0; 2 3 1 0; 2 1 3 0])
        poly1 = CDDPolyhedron{4,Rational{BigInt}}(ext1)

        poly2 = project(poly1, [1 1 1; -1 1 1; 0 -2 1; 0 0 -3])
        setobject!(vis[:polyhedron], poly2)
        settransform!(vis[:polyhedron],  Translation(-0.5, 1.0, 0) ∘ LinearMap(UniformScaling(0.1)))
    end

end

close(vis)

if proc !== nothing
    kill(proc)
end


text:
# Types using the base types

"""
    Gain(value, frequency)
    Gain(; value, frequency)

Complex type for sensitivity and frequency ranges. This complex type
can be used to represent both overall sensitivities and individual stage gains. The
`FrequencyRangeGroup` is an optional construct that defines a pass band in Hertz
(`frequency_start` and `frequency_nd`) in which the `SensitivityValue` is valid within the
number of decibels specified in FrequencyDBVariation.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw mutable struct Gain
    "A scalar that, when applied to the data values, converts the data to different units (e.g. Earth units)."
    value::Float64
    "The frequency (in Hertz) at which the Value is valid."
    frequency::Float64
end

function parse_node(::Type{Gain}, node::EzXML.Node, warn::Bool=false)
    local value, frequency
    for elm in EzXML.eachelement(node)
        if elm.name == "Value"
            value = parse(Float64, elm.content)
        elseif elm.name == "Frequency"
            frequency = parse(Float64, elm.content)
        else
            warn && @warn("unexpected element name \"$(elm.name)\" in Gain")
        end
    end
    Gain(value, frequency)
end

"""
    Sensitivity(value, frequency, input_units, output_units)
    Sensitivity(; value, frequency, input_units, output_units, frequency_start=missing, frequency_end=missing, frequency_db_variation=missing)

Sensitivity and frequency ranges.

Optionally specify a frequency range that defines a pass band in Hertz
`frequency_start` and `frequency_end` in which the `value` is valid within the
number of decibels specified in `frequency_db_variation`.

Note that either all or none of `frequency_start`, `frequency_end` and
`frequency_db_variation` must be specified upon construction.
This is **not** enforced however when mutating the fields of a `Sensitivity`
object.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw mutable struct Sensitivity
    "A scalar that, when applied to the data values, converts the data to different units (e.g. Earth units).
     Although this must be present according to the specification, it can be
     `missing` here."
    value::M{Float64} = missing
    "The frequency (in Hertz) at which the `value` is valid.
     Although this must be present according to the specification, it can be
     `missing` here."
    frequency::M{Float64} = missing
    "The units of the data as input from the perspective
	 of data acquisition. After correcting data for this response, these
	 would be the resulting units."
    input_units::Units
    "The units of the data as output from the perspective
	 of data acquisition. These would be the units of the data prior to
	 correcting for this response."
    output_units::Units
    frequency_start::M{Float64} = missing
    frequency_end::M{Float64} = missing
    frequency_db_variation::M{Float64} = missing

    function Sensitivity(value, frequency, input_units, output_units,
            frequency_start, frequency_end, frequency_db_variation)
        freq_group = (frequency_start, frequency_end, frequency_db_variation)
        count(x->x===missing, freq_group) in (0, 3) ||
            throw(ArgumentError("either all or none of $freq_group must be given"))
        new(value, frequency, input_units, output_units,
            frequency_start, frequency_end, frequency_db_variation)
    end
end

"""
    Equipment(; kwargs...)

Sensor equipment installed at a Station.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw mutable struct Equipment
    type::M{String} = missing
    description::M{String} = missing
    manufacturer::M{String} = missing
    vendor::M{String} = missing
    model::M{String} = missing
    serial_number::M{String} = missing
    installation_date::M{DateTime} = missing
    removal_date::M{DateTime} = missing
    calibration_date::Vector{DateTime} = DateTime[]
    other::M{String} = missing
    "This field contains a string that should serve as a unique
	 resource identifier. This identifier can be interpreted differently depending on
	 the datacenter/software that generated the document. Also, we recommend to use
	 something like GENERATOR:Meaningful ID. As a common behaviour equipment with the
	 same ID should contains the same information/be derived from the same base
	 instruments."
    resource_id::M{String} = missing
end

attribute_fields(::Type{Equipment}) = (:resource_id,)

"""
    Log(; entry=Comment[])

Container for log entries.

!!! note
    The `Log` type is no longer present in StationXML v1.1 and later.
    It is included here for compatibility with v1.0 files.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Log
    entry::Vector{Comment} = Comment[]
end

@enumerated_struct(PZTransferFunction,
    ("LAPLACE (RADIANS/SECOND)", "LAPLACE (HERTZ)", "DIGITAL (Z-TRANSFORM)"))

"""
    PolesZeros

Response: complex poles and zeros. Corresponds to SEED blockette 53.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw mutable struct PolesZeros
    @BaseFilter
    pz_transfer_function_type::PZTransferFunction
    normalization_factor::Float64 = 1.0
    normalization_frequency::Float64
    zero::Vector{PoleZero} = PoleZero[]
    pole::Vector{PoleZero} = PoleZero[]
end

attribute_fields(::Type{PolesZeros}) = BASE_FILTER_ATTRIBUTES

"""
    NumeratorCoefficient

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct NumeratorCoefficient <: NumberType
    value::Float64
    i::M{Int} = missing
    NumeratorCoefficient(value, i=missing) = new(value, i)
end

attribute_fields(::Type{NumeratorCoefficient}) = (:i,)

function parse_node(::Type{NumeratorCoefficient}, node::EzXML.Node, warn::Bool=false)
    value = parse(Float64, node.content)
    i = missing
    for att in EzXML.eachattribute(node)
        if att.name == "i"
            i = parse(Int, att.content)
        else
            warn && @warn("unexpected attribute $(att.name) for NumeratorCoefficient")
        end
    end
    NumeratorCoefficient(value, i)
end

@enumerated_struct(FIRSymmetry,
    ("NONE", "EVEN", "ODD"))

"""
    FIR(; kwargs...)

Response: FIR filter. Corresponds to SEED blockette 61. FIR filters
are also commonly documented using the `Coefficients` type.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw mutable struct FIR
    @BaseFilter
    "Symmetry of filter.  Must be one of `\"NONE\"`, `\"EVEN\"` or `\"ODD\"`.
     See [`StationXML.FIR`](@ref)."
    symmetry::FIRSymmetry
    "Set of [`NumeratorCoefficient`](@ref StationXML.NumeratorCoefficient)s."
    numerator_coefficient::Vector{NumeratorCoefficient} = NumeratorCoefficient[]
end

attribute_fields(::Type{FIR}) = BASE_FILTER_ATTRIBUTES

@enumerated_struct(CfTransferFunction,
    ("ANALOG (RADIANS/SECOND)", "ANALOG (HERTZ)", "DIGITAL"))

"""
    Coefficient

Extension of [`FloatNoUnit`](@ref StationXML.FloatNoUnit) with an
additional `number` attribute, used as the `numerator` and `denominator`
of [`Coefficients`](@ref StationXML.Coefficients) and
[`Polynomial`](@ref StationXML.Polynomial).

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Coefficient <: NumberType
    "Value of the coefficient (no unit)"
    value::Float64
    "Absolute error in the positive direction."
    plus_error::M{Float64} = missing
    "Absolute error in the negative direction (i.e., this value should be positive)."
    minus_error::M{Float64} = missing
    "Method used to make measurement"
    measurement_method::M{String} = missing
    "Number of the coefficient"
    number::M{Int} = missing

    function Coefficient(value, plus_error=missing, minus_error=missing,
            measurement_method=missing, number=missing)
        number !== missing && number < 0 &&
            throw(ArgumentError("number must be 0 or greater"))
        new(value, plus_error, minus_error, measurement_method, number)
    end
end

attribute_fields(::Type{Coefficient}) = (:plus_error, :minus_error,
    :measurement_method, :number)

function parse_node(::Type{Coefficient}, node::EzXML.Node, warn::Bool=false)
    value = parse(Float64, node.content)
    plus_error = minus_error = measurement_method = number = missing
    for att in EzXML.eachattribute(node)
        name = att.name
        if name == "number"
            number = parse(Int, att.content)
        elseif name == "plusError"
            plus_error = parse(Float64, att.content)
        elseif name == "minusError"
            minus_error = parse(Float64, att.content)
        elseif name == "measurementMethod"
            measurement_method = att.content
        else
            warn && @warn("unexpected attribute $(att.name) for NumeratorCoefficient")
        end
    end
    Coefficient(value, plus_error, minus_error, measurement_method, number)
end

const Numerator = Coefficient
const Denominator = Coefficient

"""
    Coefficients(; kwargs...)

Response: coefficients for FIR filter. Laplace transforms or IIR
filters can be expressed using type as well but the PolesAndZerosType should be used
instead. Corresponds to SEED blockette 54.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw mutable struct Coefficients
    @BaseFilter
    "Type of the transfer function (see [`CfTransferFunction`](@ref)."
    cf_transfer_function_type::CfTransferFunction
    numerator::Vector{Numerator} = Numerator[]
    denominator::Vector{Denominator} = Denominator[]
end

attribute_fields(::Type{Coefficients}) = BASE_FILTER_ATTRIBUTES

"""
    ResponseListElement(frequency, amplitude, phase)

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct ResponseListElement
    "Frequency at which the `amplitude` and `phase` was measured."
    frequency::Frequency
    "Amplitude response at the given `frequency`."
    amplitude::Float
    "Phase angle response at the given `frequency`"
    phase::Angle
end

"""
    ResponseList

Response: list of frequency, amplitude and phase values. Corresponds to SEED blockette 55.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct ResponseList
    @BaseFilter
    "Set of [`ResponseListElement`](@ref)s."
    response_list_element::Vector{ResponseListElement} = ResponseListElement[]
end

attribute_fields(::Type{ResponseList}) = BASE_FILTER_ATTRIBUTES

@enumerated_struct(Approximation, ("MACLAURIN",))

"""
    Polynomial

Response: expressed as a polynomial (allows non-linear sensors to be
described). Corresponds to SEED blockette 62. Can be used to describe a stage of
acquisition or a complete system.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Polynomial
    @BaseFilter
    approximation_type::Approximation = Approximation("MACLAURIN")
    frequency_lower_bound::Frequency
    frequency_upper_bound::Frequency
    "The low frequency corner for which the sensor is valid."
    approximation_lower_bound::Float64
    "The high frequency corner for which the sensor is valid."
    approximation_upper_bound::Float64
    maximum_error::Float64
    coefficient::Vector{Coefficient} = Coefficient[]
end

"""
    Decimation

Description of the decimation stage of a response.

Corresponds to SEED blockette 57.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Decimation
    "Sample rate of the signal before decimation."
    input_sample_rate::Frequency
    "Number of times the signal is decimated.  The output sample rate
     is therefore `input_sample_rate ÷ factor`."
    factor::Int
    "Which sample is chosen for use when decimating, starting at 0 for
     the first sample in the window of `factor` samples.
     (This implies `0 <= offset < factor`.)"
    offset::Int
    "The estimated pure delay of the signal, where positive values
     mean a time delay or shift to later in time."
    delay::Float64
    "Time shift applied to correct for the delay introduced in this
     stage.  Positive numbers mean a time advance or shift to earlier in time."
    correction::Float64
end

"""
    ResponseStage

Represents a channel's response and covers SEED blockettes 53 to 56.

!!! note
    Although the `stage_gain` field is mandatory in v1.0 of the StationXML specification,
    it is common in StationXML files that the field is absent for certain
    types of channel, such as state-of-health channels.  Hence it is an
    optional field of `ResponseStage` in StationXML.jl.
    Version 1.1 of the specification removed the requirement for a `stage_gain`
    when a `polynomial` stage is defined, rather than one of `poles_zeros`,
    `coefficients`, `response_list` or `fir`.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct ResponseStage
    "A choice of `polynomial` versus all other response types.
     There should be one response per stage."
    poles_zeros::M{PolesZeros} = missing
    coefficients::M{Coefficients} = missing
    response_list::M{ResponseList} = missing
    fir::M{FIR} = missing
    polynomial::M{Polynomial} = missing
    decimation::M{Decimation} = missing
    "The gain at the stage of the encapsulating
     response element and corresponds to SEED blockette 58. In the SEED
     convention, stage 0 gain represents the overall sensitivity of the channel.
     In this schema, stage 0 gains are allowed but are considered deprecated.
     Overall sensitivity should be specified in the InstrumentSensitivity
     element."
    stage_gain::M{Gain} = missing
    "Stage sequence number. This is used in all the response SEED blockettes."
    number::Int
    "Same meaning as the `resource_id` field of [`Equipment`](@ref StationXML.Equipment)."
    resource_id::M{String} = missing
    function ResponseStage(poles_zeros, coefficients, response_list, fir, polynomial,
        decimation, stage_gain, number, resource_id)
        count(!ismissing, (poles_zeros, coefficients, response_list, fir, polynomial)) <= 1 ||
            throw(ArgumentError("none or only one response type can be specified"))
        number >= 0 || throw(ArgumentError("number must be 0 or more"))
        polynomial !== missing && stage_gain !== missing &&
            @warn("stage_gain cannot be specified for polynomial stages in StationXML v1.1")
        new(poles_zeros, coefficients, response_list, fir, polynomial,
            decimation, stage_gain, number, resource_id)
    end
end

attribute_fields(::Type{ResponseStage}) = (:number, :resource_id)

"""
    Response

Instrument sensitivities, or the complete system sensitivity,
can be expressed using either a sensitivity value or a polynomial. The
information can be used to convert raw data to real units at a specified
frequency or within a range of frequencies.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Response
    "The total sensitivity for a channel, representing the
     complete acquisition system expressed as a scalar. Equivalent to SEED
     stage 0 gain with (blockette 58) with the ability to specify a frequency
     range."
    instrument_sensitivity::M{Sensitivity} = missing
    "The total sensitivity for a channel, representing the
	complete acquisition system expressed as a
    [`Polynomial`](@ref StationXML.Polynomial). Equivalent to
	SEED stage 0 polynomial (blockette 62). "
    instrument_polynomial::M{Polynomial} = missing
    "Set of [`ResponseStage`](@ref StationXML.ResponseStage)s describing
     each stage in the "
    stage::Vector{ResponseStage} = ResponseStage[]
    resource_id::M{String} = missing
end

attribute_fields(::Type{Response}) = (:resource_id,)

"""
    Operator

An operating agency and associated contact persons. If
there are multiple operators, each one should be encapsulated within an
Operator tag. Since the Contact element is a generic type that
represents any contact person, it also has its own optional `agency`
field.  It is expected that typically the contact’s optional `agency`
field will match the `Operator` `agency`. Only contacts appropriate for the
enclosing element should be included in the `Operator` object.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct Operator
    "Only one agency is allowed in StationXML v1.1."
    agency::Vector{String} = String[]
    contact::Vector{Person} = Person[]
    web_site::M{String} = missing
end

"""
    DataAvailability

A description of time series data availability. This information should be
considered transient and is primarily useful as a guide for generating time
series data requests. The information for a DataAvailability:Span may be
specific to the time range used in a request that resulted in the document
or limited to the availability of data withing the request range. These details
may or may not be retained when synchronizing metadata between data centers.

# List of fields
$(DocStringExtensions.TYPEDFIELDS)
"""
@with_kw struct DataAvailability
    extent::M{DataAvailabilityExtent} = missing
    span::Vector{DataAvailabilitySpan} = DataAvailabilitySpan[]
end


text:
#factcheck is a basic testing package for Julia. Can be applied to property based testing.
#Using basic fact check funtions.
using FactCheck

facts("Messages") do
    x = [1, 2, 3, 4]
    y = [4, 2, 3, 1]
    for i in 1:4
        @fact x[i] --> y[i] "mismatch at i=$i"
    end
end
#fact test returs mismatch at i = 1 & 4

@fact 2 + 1e-5 --> roughly(2.0)
@fact 9.5 --> roughly(10; atol=1.0)
A = [2.0, 3.0]
B = (1 + 1e-6)*A
@fact A --> roughly(B)

#can test for rough estimates of outputs. Returns fact as true if in range of 1.0


text:
module Factotum

using LinearAlgebra
using Optim
using PrettyTables
using Printf
using Statistics
using StatsBase
using StatsFuns

abstract type AbstractFactorModel end

struct FactorModel{M <: AbstractMatrix, V <: AbstractVector} <: AbstractFactorModel
    "The matrix of factors"
    factors::M
    "The matrix of loadings"
    loadings::M
    "The eigenvalues of X'X"
    eigenvalues::V
    "Residual " # Why here? Maybe to use to view and update in place...
    residuals::M
    "To demean and scale the matrix X"
    center::M
    scale::M
    "The original matrix"
    X::M
    "The rescaled matrix"
    X̄::M
end

struct FactorModelView{M <: AbstractMatrix, S <: AbstractMatrix, V <: AbstractVector, R} <: AbstractFactorModel
    "The matrix of factors"
    factors::M
    "The matrix of loadings"
    loadings::S
    "The eigenvalues of X'X"
    eigenvalues::V
    "Residual " # Why here? Maybe to use to view and update in place...
    residuals::R
    "The rescaled matrix"
    X̄::R
end

function FactorModel(Z::AbstractMatrix{G}, numfactors; kwargs...) where G
    T, n = size(Z)
    (F, Λ, λ, ε, μ, σₓ, Z, X) = T > n ? extract_ΛΛ(Z, numfactors; kwargs...) : extract_FF(Z; kwargs...)
    FactorModel(F, Λ, λ, ε, μ, σₓ, Z, X)
end

function extract_ΛΛ(Z, numfactors; demean::Bool = true, scale::Bool = false, corrected::Bool = false)
    ## Estimate
    ## X = FΛ' + e
    ##
    ## Λ'Λ = I
    ## F'F = V, V (rxr) diagonal
    T, n = size(Z)
    μ = demean ? mean(Z; dims = 1) : zeros(1, n)
    σₓ = scale ? std(Z; dims = 1, corrected = corrected) : ones(1, n)
    X = (Z .- μ)./σₓ
    ev = eigen(Symmetric(X'*X), n-numfactors+1:n)
    neg = findall(x -> x < 0, ev.values)
    if !isempty(neg)
        if any(ev.values[neg] .< -9 * eps(Float64) * first(ev.values))
            error("covariance matrix is not non-negative definite")
        else
            ev.values[neg] = 0.0
        end
    end
    λ  = ev.values[numfactors:-1:1]
    Λ = sqrt(n)*ev.vectors[:, numfactors:-1:1]
    F = X*Λ/n
    ε = (X .- F*Λ')
    (F, Λ, λ, ε, μ, σₓ, Z, X)
end

function Base.view(fm::FactorModel, k::Int)
    @assert k > 0 "Cannot view a FactorModel with 0 factors"
    view(fm, 1:k)
end

function Base.view(fm::FactorModel, rnge::UnitRange)
    T, n = size(fm)
    @assert numfactors(fm) >= maximum(rnge) "Cannot create a FactorModel's view with $(maximum(rnge)) when the parent has $(numfactors(fm)) factors"
    @assert first(rnge) <= maximum(rnge)
    ## To do: check that max(range) < numfactors fm
    FactorModelView(view(factors(fm), :, rnge), view(loadings(fm), :, rnge),
        view(eigvals(fm), rnge), residuals(fm), fm.X̄)
end

## ------------------------------------------------------------
## Methods
## ------------------------------------------------------------
Base.size(fm::AbstractFactorModel) = size(fm.X̄)
numfactors(fm::AbstractFactorModel) = size(loadings(fm), 2)
loadings(fm::AbstractFactorModel) = fm.loadings
factors(fm::AbstractFactorModel) = fm.factors

LinearAlgebra.eigvals(fm::AbstractFactorModel) = fm.eigenvalues

function sdev(fm::AbstractFactorModel)
    T, n = size(fm)
    λ = eigvals(fm)
    sqrt.(λ/(T*n))
end

function explained_variance(fm::FactorModel)
    λ = eigvals(fm)
    λ./sum(λ)
end

function StatsBase.residuals(fm::AbstractFactorModel)
    F = factors(fm)
    Λ = loadings(fm)
    fm.residuals .= fm.X̄ .- F*Λ'
    return fm.residuals
end

X(fm::FactorModel) = fm.X
X̄(fm::AbstractFactorModel) = fm.X̄

## Output
function Base.show(io::IO, fm::AbstractFactorModel)
    printstyled(io, "\nStatic Factor Model\n", color = :green)
    #@printf io "------------------------------------------------------\n"
    @printf io "Dimensions of X..........: %s\n" size(fm)
    @printf io "Number of factors........: %s\n" numfactors(fm)
    @printf io "Factors calculated by....: %s\n" "Principal Component"
    #@printf io "\n"
    #@printf io "------------------------------------------------------\n"
end

describe(fm::FactorModel) = describe(stdout, fm)

function describe(io::IO, fm::FactorModel)
    show(io, fm)
    printstyled(io, "Factors' importance:\n", color = :green)
    factortable(io, fm)
end

function factortable(io::IO, fm::FactorModel)
    explainedvar = explained_variance(fm)
    colnms = "Factor_".*string.(1:numfactors(fm))
    rownms = ["Standard deviation", "Proportion of Variance", "Cumulative Proportion"]
    mat = vcat(sdev(fm)', explainedvar', cumsum(explainedvar)')
    ct = CoefTable(mat, colnms, rownms)
    show(io, ct)
end

## ------------------------------------------------------------
## Information Criteria
## ------------------------------------------------------------
abstract type AbstractInformationCriterion end
struct IC1 <: AbstractInformationCriterion end
struct IC2 <: AbstractInformationCriterion end
struct IC3 <: AbstractInformationCriterion end
struct PCp1 <: AbstractInformationCriterion end
struct PCp2 <: AbstractInformationCriterion end
struct PCp3 <: AbstractInformationCriterion end
struct AIC1 <: AbstractInformationCriterion end
struct AIC2 <: AbstractInformationCriterion end
struct AIC3 <: AbstractInformationCriterion end
struct BIC1 <: AbstractInformationCriterion end
struct BIC2 <: AbstractInformationCriterion end
struct BIC3 <: AbstractInformationCriterion end

struct InformationCriterion{M <: AbstractInformationCriterion, T<:AbstractFloat}
    criterion::M
    crit::Array{T,1}
    rnge::UnitRange{Int64} ## Change name of this field -> rnge
end

## Calculate V(F̂ᵏ) for k ⩽ kₘₐₓ
function V(fmv::FactorModelView)
    ε = residuals(fmv)
    Λ = Factotum.loadings(fmv)
    F = Factotum.factors(fmv)
    mean(ε.^2)
end

V(fm::FactorModel, kₘₐₓ) = [V(view(fm, j)) for j ∈ 1:kₘₐₓ]

variance_factor(::Type{M}, fm, kₘₐₓ) where M <: Union{IC1, IC2, IC3} = 1.0
variance_factor(::Type{M}, fm, kₘₐₓ) where M <: AbstractInformationCriterion = V(view(fm, kₘₐₓ)) ## This is probably transform
transform_V(::Type{M}, V) where M <: Union{IC1, IC2, IC3} = log.(V)
transform_V(::Type{M}, V) where M <: AbstractInformationCriterion = V

function informationcriterion(s::Type{M}, fm::FactorModel, kₘₐₓ::Int64) where M <: AbstractInformationCriterion
    T, n = size(fm)
    rnge = 1:kₘₐₓ
    σ̂²  = variance_factor(s, fm, kₘₐₓ)
    Vₖ = transform_V.(s, V(fm, kₘₐₓ))
    gₜₙ = map(k -> k*penalty(s, T, n, k), rnge)
    InformationCriterion(M(), Vₖ + σ̂².*gₜₙ, rnge)
end

function informationcriteria(criterion::Tuple, fm, kₘₐₓ)
    @assert all(map(x->isa(x(), Factotum.AbstractInformationCriterion), criterion)) "Some of the arguments is not a SelectionCriterion"
    map(x->x(fm, kₘₐₓ), criterion)
end

(for criterion ∈ (:BIC1, :BIC2, :BIC3, :AIC1, :AIC2, :AIC3, :IC1, :IC2, :IC3, :PCp1, :PCp2, :PCp3)
    eval(quote
    ($criterion)(fm, kₘₐₓ::Int64) = Factotum.informationcriterion($criterion, fm, kₘₐₓ)
    Base.string(ic::($criterion)) = string(($criterion))
    ($criterion)(X::Matrix, kₘₐₓ::Int64; kwargs...) = Factotum.informationcriterion($criterion, FactorModel(X, kₘₐₓ; kwargs...), kₘₐₓ)
    end)
end)

function Base.findmin(ic::InformationCriterion)
    fmin = findmin(ic.crit)
    NamedTuple{(Symbol(string(ic.criterion)), :r)}(fmin)
end

function Base.findmin(ic::Tuple{Vararg{InformationCriterion, N}}) where {N}
    fmin = map(x->findmin(x.crit), ic)
    nm = map(x->(Symbol(x.criterion), :r), ic)
    TT = map(x->eltype(x.crit), ic)
    map((nm,x, T) -> NamedTuple{nm}(x), nm, fmin, TT)
end

numfactors(ic::InformationCriterion) = findmin(ic).r
Base.string(ic::InformationCriterion{T, F}) where {T, F}  = string(T)

function Base.show(io::IO, ic::T) where T <: InformationCriterion
    header = ["# of factor"  "Criterion"]
    highlight1 = Highlighter((data, i, j) -> data[i,2] == minimum(data[:,2]), Crayon(background = :light_blue, foreground = :white, bold = true))
    highlight2 = Highlighter((data, i, j) -> (j == 2), Crayon(foreground = :light_blue))
    highlight3 = Highlighter((data, i, j) -> (j == 1), Crayon(foreground = :light_red, bold = true))
    pretty_table([ic.rnge ic.crit], header, header_crayon = crayon"yellow bold", formatters = (ft_printf("%5.0f", 1), ft_printf("%5.3f", 2:2)), highlighters = (highlight1, highlight2, highlight3))
end

function Base.show(io::IO, ic::Tuple{Vararg{InformationCriterion, N}}) where {N}
    header = ["# of factor"  string.(ic)...]
    highlights = map(x->Highlighter((data, i, j) -> data[i,j] == minimum(data[:,x]) && j > 1, Crayon(background = :light_blue, foreground = :white, bold = true)), 2:length(ic)+1)
    tbl = [first(ic).rnge mapreduce(x->x.crit, hcat, ic)]
    pretty_table(tbl, header, header_crayon = crayon"yellow bold", formatters = (ft_printf("%5.0f", 1), ft_printf("%5.3f", 2:length(ic)+1)), highlighters = (highlights...,))
end

function penalty(s::Type{P}, T, N) where P <: Union{IC1, PCp1}
    NtT = N*T
    NpT = N+T
    p1 = NpT/NtT
    p2 = log(NtT/NpT)
    p1*p2
end

function penalty(s::Type{P}, T, N) where P <: Union{IC2, PCp2}
    C2  = min(T, N)
    NtT = N*T
    NpT = N+T
    p1 = NpT/NtT
    p2 = log(C2)
    p1*p2
end

function penalty(s::Type{P}, T, N) where P <: Union{IC3, PCp3}
    C2  = min(T, N)
    log(C2)/C2
end

penalty(s::Type{S}, T, N, k) where S = penalty(s, T, N)

penalty(s::Type{AIC1}, T, N) = 2/T
penalty(s::Type{AIC2}, T, N) = 2/N
penalty(s::Type{AIC3}, T, N, k) = 2*(N+T-k)/(N*T)

penalty(s::Type{BIC1}, T, N) = log(T)/T
penalty(s::Type{BIC2}, T, N) = log(N)/N
penalty(s::Type{BIC3}, T, N, k) = ((N+T-k)*log(N*T))/(N*T)

############################################################
## Wald test
############################################################
# struct WaldTest
#     tbl::NamedTuple
#     rankmin::Int64
#     rankₘₐₓ::Int64
# end

# struct WaldTestFun{F, T, Z}
#     f::F
#     r::Int64
#     vecsigma::T
#     Vhat::Z
# end

# (wf::WaldTestFun)(theta) = wf.f(theta, wf.r, wf.vecsigma, wf.Vhat)

# function waldobjfun(th, r, vecsigma, Vhat)
#     ##r,k = size(theta) ## note that the rank being tested is r0 = r-1
#     theta = reshape(th, r+1, length(th)÷(r+1))
#     sigmamat = diagm(0=>theta[1,:].^2) .+ theta[2:r+1,:]'*theta[2:r+1,:]
#     tempsigma = sigmamat[findall(tril(ones(size(sigmamat))).==1)]
#     (vecsigma -tempsigma)' /Vhat *(vecsigma - tempsigma)
# end

# X = randn(100,10);
# fm = Factotum.FactorModel(X)

# function waldtest(fm::FactorModel, minrank::Int = 0, maxrank::Int = 2)
#     X = copy(fm.X)
#     T, n = size(X)
#     ## Normalize factor
#     Xs = X / diagm(0=>sqrt.(diag(cov(X))))
#     covX = cov(Xs)
#     meanX = mean(Xs, dims=1)
#     vecsigma = Factotum.vech(covX)
#     bigN = length(vecsigma)
#     Vhat = Array{Float64}(undef, bigN, bigN)
#     varvecsig = zeros(n,n,n,n);

#     for i1 in 1:n, i2 = 1:n, i3 = 1:n, i4 = 1:n
#         varvecsig[i1,i2,i3,i4] = sum( (Xs[:,i1] .- meanX[i1]).*(Xs[:,i2] .- meanX[i2]).*(Xs[:,i3] .- meanX[i3]) .*(Xs[:,i4] .- meanX[i4])) / T^2 - covX[i1,i2] *covX[i3,i4] /T
#     end

#     idx = findall(tril(ones(size(covX))).==1)
#     for i=1:bigN, j=1:bigN   ## map elements of varvecsig array into matrix corresponding to
#         Vhat[i,j] = varvecsig[idx[i],idx[j]]
#     end

#     out_table = (rank = -1, waldstat = NaN, df = NaN, pvalue = NaN)

#     ## Initial values
#     for k in minrank:maxrank
#         wf = WaldTestFun(waldobjfun, k, vecsigma, Vhat)
#         df = (n-k)*(n-k+1)/2 - n

#         theta0 = theta_initial_value(n,k)

#         outs = Array{Tuple{Float64,Array{Float64,1},Bool},1}()

#         for j in theta0
#             out = Optim.optimize(wf, j, BFGS(), Optim.Options(allow_f_increases=true); autodiff=:forward)
#             push!(outs, (out.minimum::Float64, out.minimizer::Array{Float64,1}, Optim.converged(out)::Bool))
#         end

#         convouts = outs[map(x->x[3], outs)]
#         out      = convouts[argmin(map(x->x[1], convouts))]

#         dfa = (rank = k, waldstat = out[1], df = df,  pvalue = 1-StatsFuns.chisqcdf(df, out[1]))
#         append!(out_table, dfa)
#     end
#     filter!(raw->raw[:rank]>=0, out_table)
#     WaldTest(out_table, minrank, maxrank)
# end

# function theta_initial_value(n,k)
#     I3 = ones(1,n)/3
#     ek = [diagm(0=>ones(k)) zeros(k, n-k)]
#     t0 = ([I3; zeros(k,n)], [I3; ones(k,n)./(2*k)],[I3; ek./(2*k)], [I3; reverse(ek./(2*k), dims=2)])::NTuple{4,Array{Float64,2}}
#     map(vec, t0)::NTuple{4,Array{Float64,1}}
# end

# function vech(X::Matrix{S}) where S
#     T, n = size(X)
#     r = round(Int64, n*(n+1)/2)
#     x = Array{S, 1}(undef, r)
#     i = 1
#     for j in 1:n, k in j:n
#         x[i] = X[j,k]
#         i += 1
#     end
#     x
# end

export FactorModel, subview, waldtest, describe, PValue, waldstat, Criteria,
       IC1, IC2, IC3, PCp1, PCp2, PCp3,
       AIC1, AIC2, AIC3, BIC1, BIC2, BIC3

end # module"


text:
export BrickletIMUV3Acceleration
struct BrickletIMUV3Acceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickletIMUV3MagneticField
struct BrickletIMUV3MagneticField
    x::Integer
    y::Integer
    z::Integer
end

export BrickletIMUV3AngularVelocity
struct BrickletIMUV3AngularVelocity
    x::Integer
    y::Integer
    z::Integer
end

export BrickletIMUV3Orientation
struct BrickletIMUV3Orientation
    heading::Integer
    roll::Integer
    pitch::Integer
end

export BrickletIMUV3LinearAcceleration
struct BrickletIMUV3LinearAcceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickletIMUV3GravityVector
struct BrickletIMUV3GravityVector
    x::Integer
    y::Integer
    z::Integer
end

export BrickletIMUV3Quaternion
struct BrickletIMUV3Quaternion
    w::Integer
    x::Integer
    y::Integer
    z::Integer
end

export BrickletIMUV3AllData
struct BrickletIMUV3AllData
    acceleration::Vector{Integer}
    magnetic_field::Vector{Integer}
    angular_velocity::Vector{Integer}
    euler_angle::Vector{Integer}
    quaternion::Vector{Integer}
    linear_acceleration::Vector{Integer}
    gravity_vector::Vector{Integer}
    temperature::Integer
    calibration_status::Integer
end

export BrickletIMUV3SensorConfiguration
struct BrickletIMUV3SensorConfiguration
    magnetometer_rate::Integer
    gyroscope_range::Integer
    gyroscope_bandwidth::Integer
    accelerometer_range::Integer
    accelerometer_bandwidth::Integer
end

export BrickletIMUV3AccelerationCallbackConfiguration
struct BrickletIMUV3AccelerationCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3MagneticFieldCallbackConfiguration
struct BrickletIMUV3MagneticFieldCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3AngularVelocityCallbackConfiguration
struct BrickletIMUV3AngularVelocityCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3TemperatureCallbackConfiguration
struct BrickletIMUV3TemperatureCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3OrientationCallbackConfiguration
struct BrickletIMUV3OrientationCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3LinearAccelerationCallbackConfiguration
struct BrickletIMUV3LinearAccelerationCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3GravityVectorCallbackConfiguration
struct BrickletIMUV3GravityVectorCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3QuaternionCallbackConfiguration
struct BrickletIMUV3QuaternionCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3AllDataCallbackConfiguration
struct BrickletIMUV3AllDataCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIMUV3SPITFPErrorCount
struct BrickletIMUV3SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIMUV3Identity
struct BrickletIMUV3Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIMUV3
"""
Full fledged AHRS with 9 degrees of freedom
"""
mutable struct BrickletIMUV3 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIMUV3(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_imu_v3")
        deviceInternal = package.BrickletIMUV3(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_acceleration
"""
    $(SIGNATURES)

Returns the calibrated acceleration from the accelerometer for the
x, y and z axis. The acceleration is in the range configured with
:func:`Set Sensor Configuration`.

If you want to get the acceleration periodically, it is recommended
to use the :cb:`Acceleration` callback and set the period with
:func:`Set Acceleration Callback Configuration`.
"""
function get_acceleration(device::BrickletIMUV3)
    return device.deviceInternal.get_acceleration()
end

export get_magnetic_field
"""
    $(SIGNATURES)

Returns the calibrated magnetic field from the magnetometer for the
x, y and z axis.

If you want to get the magnetic field periodically, it is recommended
to use the :cb:`Magnetic Field` callback and set the period with
:func:`Set Magnetic Field Callback Configuration`.
"""
function get_magnetic_field(device::BrickletIMUV3)
    return device.deviceInternal.get_magnetic_field()
end

export get_angular_velocity
"""
    $(SIGNATURES)

Returns the calibrated angular velocity from the gyroscope for the
x, y and z axis. The angular velocity is in the range configured with
:func:`Set Sensor Configuration`.

If you want to get the angular velocity periodically, it is recommended
to use the :cb:`Angular Velocity` acallback nd set the period with
:func:`Set Angular Velocity Callback Configuration`.
"""
function get_angular_velocity(device::BrickletIMUV3)
    return device.deviceInternal.get_angular_velocity()
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of the IMU Brick.
The temperature is measured in the core of the BNO055 IC, it is not the
ambient temperature
"""
function get_temperature(device::BrickletIMUV3)
    return device.deviceInternal.get_temperature()
end

export get_orientation
"""
    $(SIGNATURES)

Returns the current orientation (heading, roll, pitch) of the IMU Brick as
independent Euler angles. Note that Euler angles always
experience a `gimbal lock <https://en.wikipedia.org/wiki/Gimbal_lock>`__.
We recommend that you use quaternions instead, if you need the absolute
orientation.

If you want to get the orientation periodically, it is recommended
to use the :cb:`Orientation` callback and set the period with
:func:`Set Orientation Callback Configuration`.
"""
function get_orientation(device::BrickletIMUV3)
    return device.deviceInternal.get_orientation()
end

export get_linear_acceleration
"""
    $(SIGNATURES)

Returns the linear acceleration of the IMU Brick for the
x, y and z axis. The acceleration is in the range configured with
:func:`Set Sensor Configuration`.

The linear acceleration is the acceleration in each of the three
axis of the IMU Brick with the influences of gravity removed.

It is also possible to get the gravity vector with the influence of linear
acceleration removed, see :func:`Get Gravity Vector`.

If you want to get the linear acceleration periodically, it is recommended
to use the :cb:`Linear Acceleration` callback and set the period with
:func:`Set Linear Acceleration Callback Configuration`.
"""
function get_linear_acceleration(device::BrickletIMUV3)
    return device.deviceInternal.get_linear_acceleration()
end

export get_gravity_vector
"""
    $(SIGNATURES)

Returns the current gravity vector of the IMU Brick for the
x, y and z axis.

The gravity vector is the acceleration that occurs due to gravity.
Influences of additional linear acceleration are removed.

It is also possible to get the linear acceleration with the influence
of gravity removed, see :func:`Get Linear Acceleration`.

If you want to get the gravity vector periodically, it is recommended
to use the :cb:`Gravity Vector` callback and set the period with
:func:`Set Gravity Vector Callback Configuration`.
"""
function get_gravity_vector(device::BrickletIMUV3)
    return device.deviceInternal.get_gravity_vector()
end

export get_quaternion
"""
    $(SIGNATURES)

Returns the current orientation (w, x, y, z) of the IMU Brick as
`quaternions <https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation>`__.

You have to divide the return values by 16383 (14 bit) to get
the usual range of -1.0 to +1.0 for quaternions.

If you want to get the quaternions periodically, it is recommended
to use the :cb:`Quaternion` callback and set the period with
:func:`Set Quaternion Callback Configuration`.
"""
function get_quaternion(device::BrickletIMUV3)
    return device.deviceInternal.get_quaternion()
end

export get_all_data
"""
    $(SIGNATURES)

Return all of the available data of the IMU Brick.

* acceleration (see :func:`Get Acceleration`)
* magnetic field (see :func:`Get Magnetic Field`)
* angular velocity (see :func:`Get Angular Velocity`)
* Euler angles (see :func:`Get Orientation`)
* quaternion (see :func:`Get Quaternion`)
* linear acceleration (see :func:`Get Linear Acceleration`)
* gravity vector (see :func:`Get Gravity Vector`)
* temperature (see :func:`Get Temperature`)
* calibration status (see below)

The calibration status consists of four pairs of two bits. Each pair
of bits represents the status of the current calibration.

* bit 0-1: Magnetometer
* bit 2-3: Accelerometer
* bit 4-5: Gyroscope
* bit 6-7: System

A value of 0 means for "not calibrated" and a value of 3 means
"fully calibrated". In your program you should always be able to
ignore the calibration status, it is used by the calibration
window of the Brick Viewer and it can be ignored after the first
calibration. See the documentation in the calibration window for
more information regarding the calibration of the IMU Brick.

If you want to get the data periodically, it is recommended
to use the :cb:`All Data` callback and set the period with
:func:`Set All Data Callback Configuration`.
"""
function get_all_data(device::BrickletIMUV3)
    return device.deviceInternal.get_all_data()
end

export save_calibration
"""
    $(SIGNATURES)

A call of this function saves the current calibration to be used
as a starting point for the next restart of continuous calibration
of the IMU Brick.

A return value of *true* means that the calibration could be used and
*false* means that it could not be used (this happens if the calibration
status is not "fully calibrated").

This function is used by the calibration window of the Brick Viewer, you
should not need to call it in your program.
"""
function save_calibration(device::BrickletIMUV3)
    return device.deviceInternal.save_calibration()
end

export set_sensor_configuration
"""
    $(SIGNATURES)

Sets the available sensor configuration for the Magnetometer, Gyroscope and
Accelerometer. The Accelerometer Range is user selectable in all fusion modes,
all other configurations are auto-controlled in fusion mode.
"""
function set_sensor_configuration(device::BrickletIMUV3, magnetometer_rate, gyroscope_range, gyroscope_bandwidth, accelerometer_range, accelerometer_bandwidth)
    device.deviceInternal.set_sensor_configuration(magnetometer_rate, gyroscope_range, gyroscope_bandwidth, accelerometer_range, accelerometer_bandwidth)
end

export get_sensor_configuration
"""
    $(SIGNATURES)

Returns the sensor configuration as set by :func:`Set Sensor Configuration`.
"""
function get_sensor_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_sensor_configuration()
end

export set_sensor_fusion_mode
"""
    $(SIGNATURES)

If the fusion mode is turned off, the functions :func:`Get Acceleration`,
:func:`Get Magnetic Field` and :func:`Get Angular Velocity` return uncalibrated
and uncompensated sensor data. All other sensor data getters return no data.

Since firmware version 2.0.6 you can also use a fusion mode without magnetometer.
In this mode the calculated orientation is relative (with magnetometer it is
absolute with respect to the earth). However, the calculation can't be influenced
by spurious magnetic fields.

Since firmware version 2.0.13 you can also use a fusion mode without fast
magnetometer calibration. This mode is the same as the normal fusion mode,
but the fast magnetometer calibration is turned off. So to find the orientation
the first time will likely take longer, but small magnetic influences might
not affect the automatic calibration as much.
"""
function set_sensor_fusion_mode(device::BrickletIMUV3, mode)
    device.deviceInternal.set_sensor_fusion_mode(mode)
end

export get_sensor_fusion_mode
"""
    $(SIGNATURES)

Returns the sensor fusion mode as set by :func:`Set Sensor Fusion Mode`.
"""
function get_sensor_fusion_mode(device::BrickletIMUV3)
    return device.deviceInternal.get_sensor_fusion_mode()
end

export set_acceleration_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Acceleration` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_acceleration_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_acceleration_callback_configuration(period, value_has_to_change)
end

export get_acceleration_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Acceleration Callback Configuration`.
"""
function get_acceleration_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_acceleration_callback_configuration()
end

export set_magnetic_field_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Magnetic Field` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_magnetic_field_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_magnetic_field_callback_configuration(period, value_has_to_change)
end

export get_magnetic_field_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Magnetic Field Callback Configuration`.
"""
function get_magnetic_field_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_magnetic_field_callback_configuration()
end

export set_angular_velocity_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Angular Velocity` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_angular_velocity_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_angular_velocity_callback_configuration(period, value_has_to_change)
end

export get_angular_velocity_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Angular Velocity Callback Configuration`.
"""
function get_angular_velocity_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_angular_velocity_callback_configuration()
end

export set_temperature_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Temperature` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_temperature_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_temperature_callback_configuration(period, value_has_to_change)
end

export get_temperature_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Temperature Callback Configuration`.
"""
function get_temperature_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_temperature_callback_configuration()
end

export set_orientation_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Orientation` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_orientation_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_orientation_callback_configuration(period, value_has_to_change)
end

export get_orientation_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Orientation Callback Configuration`.
"""
function get_orientation_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_orientation_callback_configuration()
end

export set_linear_acceleration_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Linear Acceleration` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_linear_acceleration_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_linear_acceleration_callback_configuration(period, value_has_to_change)
end

export get_linear_acceleration_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Linear Acceleration Callback Configuration`.
"""
function get_linear_acceleration_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_linear_acceleration_callback_configuration()
end

export set_gravity_vector_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Gravity Vector` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_gravity_vector_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_gravity_vector_callback_configuration(period, value_has_to_change)
end

export get_gravity_vector_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Gravity Vector Callback Configuration`.
"""
function get_gravity_vector_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_gravity_vector_callback_configuration()
end

export set_quaternion_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Quaternion` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_quaternion_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_quaternion_callback_configuration(period, value_has_to_change)
end

export get_quaternion_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Quaternion Callback Configuration`.
"""
function get_quaternion_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_quaternion_callback_configuration()
end

export set_all_data_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`All Data` callback
is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_all_data_callback_configuration(device::BrickletIMUV3, period, value_has_to_change)
    device.deviceInternal.set_all_data_callback_configuration(period, value_has_to_change)
end

export get_all_data_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set All Data Callback Configuration`.
"""
function get_all_data_callback_configuration(device::BrickletIMUV3)
    return device.deviceInternal.get_all_data_callback_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIMUV3)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIMUV3, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIMUV3)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIMUV3, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIMUV3, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIMUV3, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIMUV3)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIMUV3)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIMUV3)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIMUV3, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIMUV3)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIMUV3)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIMUV3, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Soss

mstep = @model pars,state begin
    # Parameters
    α = pars.α        # Transmission rate
    β = pars.β        # Recovery rate
    γ = pars.γ        # Case fatality rate

    # Starting counts
    s0 = state.s      # Susceptible
    i0 = state.i      # Infected
    r0 = state.r      # Recovered
    d0 = state.d      # Deceased
    n0 = state.n       # Population size

    # Transitions between states
    si ~ Binomial(s0, α * i0 / n)
    ir ~ Binomial(i0, β)
    id ~ Binomial(i0-ir, γ)

    # Updated counts
    s = s0 - si
    i = i0 + si - ir - id
    r = r0 + ir
    d = d0 + id
    n = n0 - id

    next = (pars=pars, state=(s=s,i=i,r=r,d=d, n=n))
end;

m = @model s0 begin
    α ~ Uniform()
    β ~ Uniform()
    γ ~ Uniform()
    pars = (α=α, β=β, γ=γ)
    x ~ MarkovChain(pars, mstep(pars=pars, state=s0))
end

s0 = (i = 1, r = 0, d = 0, n = 331000000, s = 330999999);

r = rand(m(s0=s0));

for (n,s) in enumerate(r.x)
    n>100 && break
    n%7==0 && println(s)
end


text:
if haskey(ENV, "CI") && ENV["CI"] == "true"
    ENV["PYTHON"] = ""
    using Pkg
    Pkg.build("PyCall")
end

using CUDA
using PyPlot
using Oceananigans

using ConvergenceTests
using ConvergenceTests.TwoDimensionalDiffusion: run_and_analyze

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

function convergence_test(Nx, Δt, stop_iteration, arch, topo)

    results = [run_and_analyze(architecture=arch, Nx=N, Δt=Δt, stop_iteration=stop_iteration,
                               topo=topo, output=false) for N in Nx]

    L₁ = map(r -> r.L₁, results)
    L∞ = map(r -> r.L∞, results)

    return (L₁=L₁, L∞=L∞)
end

# Setup and run simulations
arch = CUDA.has_cuda() ? GPU() : CPU()
Nx = [8, 16, 32, 64, 128, 256]
stop_time = 1e-4

# Calculate time step based on diffusive time-step constraint for finest mesh
        min_Δx = 2π / maximum(Nx)
   proposal_Δt = 1e-3 * min_Δx^2 # proposal time-step
stop_iteration = round(Int, stop_time / proposal_Δt)
            Δt = stop_time / stop_iteration # ensure time-stepping to exact finish time.

all_topologies = ((Periodic, Periodic, Bounded),
                  (Periodic, Bounded,  Bounded),
                  (Bounded,  Bounded,  Bounded))

topologies = arch isa CPU ? all_topologies : all_topologies[1:2]

errors = [convergence_test(Nx, Δt, stop_iteration, arch, topo) for topo in topologies]

fig, axs = subplots()

for (itopo, topo) in enumerate(topologies)
    L₁ = errors[itopo].L₁
    L∞ = errors[itopo].L∞
    name ="$(topo[1]), $(topo[2])"

    loglog(Nx, L₁, "o", color=defaultcolors[itopo], mfc="None", label="\$L_1\$-norm, $name")
    loglog(Nx, L∞, "^", color=defaultcolors[itopo], mfc="None", label="\$L_\\infty\$-norm, $name")
end

L₁ = errors[1].L₁
L∞ = errors[1].L∞
loglog(Nx, L₁[1] * (Nx[1]./Nx).^2, "k-", alpha=0.6, linewidth=1, label=L"\sim N_x^{-2}")

xscale("log", base=2)
yscale("log", base=10)
xlabel(L"N_x")
ylabel("\$L\$-norms of \$ | c_\\mathrm{sim} - c_\\mathrm{analytical} |\$")
title("Two dimensional diffusion convergence test")
removespines("top", "right")
legend(loc="upper right")

filename = "two_dimensional_diffusion_convergence_$(typeof(arch)).png"
filepath = joinpath(@__DIR__, "figs", filename)
mkpath(dirname(filepath))
savefig(filepath, dpi=480)

for (itopo, topo) in enumerate(topologies)
    L₁_topo = errors[itopo].L₁
    L∞_topo = errors[itopo].L∞
    test_rate_of_convergence(L₁_topo, Nx, expected=-2.0, atol=0.01, name="2D diffusion $topo L₁")
    test_rate_of_convergence(L∞_topo, Nx, expected=-2.0, atol=0.06, name="2D diffusion $topo L∞")
end


text:
"""
    julia2onnx(x)

Convert argument from Julia-friendly to ONNX-friendly format.
For example:

 * tensors are reshaped from column-major to row-major
 * shape tuples are rearranged accordingly
 * attributes in dicts are renamed, their values are converted too

 The reverse operation is available as [onnx2julia](@ref).
"""
julia2onnx(x::AbstractArray) = reshape(x, reverse(size(x)))

"""
    onnx2julia(x)

The reverse of [julia2onnx](@ref).
"""
onnx2julia(x::AbstractArray) = reshape(x, reverse(size(x)))


text:
const FT = Float64
import SciMLBase: step!
using OrdinaryDiffEq
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33, savevalues!, Euler
using LinearAlgebra
using ClimaCore

# pending PR merge:
# import Pkg; Pkg.add(url="https://github.com/CliMA/ClimaCore.jl",rev="sb/online-sphere-remap", subdir = "lib/ClimaCoreTempestRemap")

# using old interface - before the new one is consolidated:
import Pkg;
Pkg.add(url = "https://github.com/CliMA/ClimaAtmos.jl", rev = "jh_ln/external_boundary_conditons");
using ClimaCoreTempestRemap

import Test: @test
import ClimaAtmos.BoundaryConditions: get_boundary_flux, AbstractBoundary
export get_boundary_flux

include("atmos.jl")
include("slab.jl")
include("coupled_bc.jl")
include("conservation.jl")

# initiate spatial and temporal info
tspan = (0.0, 1300.0) # 172800.0)
cpl_Δt = 5.0 # seconds
atmos_dt, slab_dt = (5.0, 5.0) # seconds

slab_sim = slab_init(FT, tspan, dt = atmos_dt, npolynomial = 3)
atmos_sim = atmos_init(FT, tspan, dt = slab_dt, npolynomial = 3)

# parameter exchange
slab_sim.integrator.p.cp_d .= CLIMAParameters.Planet.cp_d(atmos_sim.model.parameters)

# regridding init
atmos_energy_bc = atmos_sim.model.boundary_conditions.ρe_tot.bottom
atmos_momentum_bc = atmos_sim.model.boundary_conditions.uh.bottom

F_a_space = axes(atmos_energy_bc.flux)
F_s_space = axes(slab_sim.integrator.u.T_sfc)

# weightfile = tempname()
# R_atm2slab = ClimaCoreTempestRemap.generate_map( # bring back once CC #614 is resolved/ merged
#     F_s_space, #target
#     F_a_space, #source
#     weightfile = weightfile,
# )
# R_slab2atm = ClimaCoreTempestRemap.generate_map(# bring back once CC #614 is resolved/ merged
#     F_a_space, #target
#     F_s_space, #source
#     weightfile = weightfile,
# )
function dummmy_remap!(target, source)  # delete when can use Tempest again
    parent(target) .= parent(source)
end

# init conservation info collector
CS = ConservationCheck([], [])

# coupling loop
for t in (tspan[1]:cpl_Δt:tspan[end])
    @show t

    ## Atmos
    # calculate surface fluxes for Atmos BCs
    coupler_atmos_boundary_flux(atmos_energy_bc, atmos_sim, slab_sim, F_a_space, F_s_space)
    coupler_atmos_boundary_flux(atmos_momentum_bc, atmos_sim, slab_sim, F_a_space, F_s_space)

    # run
    step!(atmos_sim.integrator, t - atmos_sim.integrator.t, true) # NOTE: use (t - integ_atm.t) here instead of Δt_cpl to avoid accumulating roundoff error in our timestepping.

    ## Slab
    # pre: get accumulated flux from atmos
    F_S = ClimaCore.Fields.zeros(F_s_space)

    # ClimaCoreTempestRemap.remap!(F_S, atmoatmos_energy_bcs_bc.flux, R_atm2slab)
    dummmy_remap!(F_S, atmos_energy_bc.flux)

    # save the accumulated flux
    slab_F_sfc = slab_sim.integrator.p.F_sfc
    slab_F_sfc .= F_S .* cpl_Δt

    # run
    step!(slab_sim.integrator, t - slab_sim.integrator.t, true)

    # conservation info
    check_conservation(CS, atmos_sim, slab_sim)
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
#=
     sqrt(eps(Float64)) == ldexp(0.5, -25) ~ 1.5e-8
     isapprox here uses ldexp(0.5, -30) ~ 1.0e-9 ,   so the system default (sqrt(eps(T))) is ~32x less sensitive 5bits
                        ldexp(0.5, -31) ~ 2.3e-10,                                           ~64x less sensitive 6bits
                                                                                (32.5 sig bits must match)

=#
releps(::Type{Float64}) = ldexp(0.5, -31)
releps(::Type{Float32}) = ldexp(inv(sqrt(2.0f0)), -17)
releps(::Type{Double64}) = ldexp(0.5, -62)
releps(::Type{Double32}) = ldexp(inv(sqrt(2.0f0)), -34)

function Base.isapprox(x::DoubleFloat{T}, y::T; atol::Real=0.0, rtol::Real=atol>0 ? 0 : releps(T), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat}
    return isapprox(x, DoubleFloat{T}(y), atol=atol, rtol=rtol, nans=nans, norm=norm)
end
function Base.isapprox(x::T, y::DoubleFloat{T}; atol::Real=0.0, rtol::Real=atol>0 ? 0 : releps(T), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat}
    return isapprox(DoubleFloat{T}(x), y, atol=atol, rtol=rtol, nans=nans, norm=norm)
end
function Base.isapprox(x::DoubleFloat{T}, y::DoubleFloat{T}; atol::Real=0.0, rtol::Real=atol>0 ? 0 : releps(DoubleFloat{T}), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat}
    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= max(max(1.0e-32, atol), rtol*max(abs(x), abs(y)))) || (nans && isnan(x) && isnan(y))
end

Base.isapprox(x::DoubleFloat{T}, y::F; atol::Real=0.0, rtol::Real=atol>0.0 ? 0.0 : eps(max(abs(x), abs(y)))^(37/64), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat, F<:Real} =
    isapprox(promote(x, y)..., atol=atol, rtol=rtol, nans=nans, norm=norm)

Base.isapprox(x::F, y::DoubleFloat{T}; atol::Real=0.0, rtol::Real=atol>0.0 ? 0.0 : eps(max(abs(x), abs(y)))^(37/64), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat, F<:Real} =
    isapprox(promote(x, y)..., atol=atol, rtol=rtol, nans=nans, norm=norm)

function Base.lerpi(j::Integer, d::Integer, a::DoubleFloat{T}, b::DoubleFloat{T}) where {T}
    t = DoubleFloat{T}(j)/d
    a = fma(-t, a, a)
    return fma(t, b, a)
end

function Base.clamp(x::DoubleFloat{T}, lo::DoubleFloat{T}, hi::DoubleFloat{T}) where {T}
     lo <= x <= hi && return x
     lo <= x && return hi
     return lo
end
function Base.clamp(x::DoubleFloat{T}, lo::T, hi::T) where {T}
     lo <= x <= hi && return x
     lo <= x && return hi
     return lo
end
function Base.clamp(x::T, lo::DoubleFloat{T}, hi::DoubleFloat{T}) where {T}
     lo <= x <= hi && return x
     lo <= x && return hi
     return lo
end

abs2(x::DoubleFloat{T}) where {T} =
    let absx = abs(x)
        absx * absx
    end

# for compatibility with old or unrevised outside linalg functions
function Base.:(+)(v::Vector{DoubleFloat{T}}, x::T) where {T}
    return v .+ x
end
function Base.:(-)(v::Vector{DoubleFloat{T}}, x::T) where {T}
    return v .- x
end
function Base.:(+)(m::Matrix{DoubleFloat{T}}, x::T) where {T}
    return m .+ x
end
function Base.:(-)(m::Matrix{DoubleFloat{T}}, x::T) where {T}
    return m .- x
end


text:
using PyPlot
using CSV
include("structures.jl")
include("utilitaires.jl")
include("parser.jl")
include("concorde.jl")
include("partitionnement_mono.jl")

function mainExp()
    # Tous les types d'instances hormis Poikonen
    run(`clear`)
    println(" Type des instances :")
    println(" --------------------")
    println("  1) Agatz")
    println("  2) Nous")
    println("  3) TSPLIB")
    println("  4) Poikonen")
    typeInstance = parse(Int64, choixListe(4))
    if typeInstance == 4
        expPoiko()
    else
        expAutres(typeInstance)
    end
end

########## Agatz, TSPLIB, Nous ##########

function expAutres(typeInstance)
    f = open("rapport.tex", "w")
    # En-tête
    l = "\\documentclass{article}"
    l *= "\n\\usepackage[utf8]{inputenc}"
    l *= "\n\\usepackage[top=2cm, bottom=3cm, left=2.5cm, right=2.5cm]{geometry}"
    l *= "\n\\usepackage{float}"
    l *= "\n\\usepackage{graphicx}"
    l *= "\n\\usepackage{rotating}"
    l *= "\n\n\\begin{document}"
    l *= "\n\n\\hrule"
    l *= "\n\\begin{center}"
    l *= "\n\\Large{\\textbf{Experimentations}}"
    l *= "\n\\vspace{0.1mm}"
    l *= "\n\\end{center}"
    l *= "\n\\hrule \\vspace{8mm}"
    write(f, l)

    # Table des matières
    #l = "\n\n\\tableofcontents"
    #write(f, l)

    # Logo
    l = "\n\n\\vspace{4mm}\\begin{figure}[H]"
    l *= "\n\\center"
    l *= "\n\\includegraphics[scale=0.22]{logo/logo+.png}"
    l *= "\n\\end{figure}"
    write(f, l)

    # Paramètres
    vitesseCamion = 1
    vitesseDrone = 3
    listeInstances = readdir("experimentations")
    nbInstances = length(listeInstances)
    l = "\n\n\\vspace{5mm}"
    l *= "\n\\section{Parameters}"
    l *= "\n\n\\begin{itemize}"
    l *= "\n    \\item Truck speed: \\textbf{"*string(vitesseCamion)*"}"
    l *= "\n    \\item Drone speed: \\textbf{"*string(vitesseDrone)*"}"
    l *= "\n    \\item Number of instances: \\textbf{"*string(nbInstances)*"}"
    l *= "\n\n\\end{itemize}"
    write(f, l)

    # Environnement
    println("\n Détection du matériel :")
    println(" -----------------------")
    os1 = read(`uname -s`, String) # Read permet de lire le résultat de la commande
    os2 = read(`uname -r`, String)
    # Commande = uname -s -r
    systeme = read(`dmidecode -s system-product-name`, String)
    # Commande = dmidecode -s system-product-name
    processeur = read(`dmidecode -s processor-version`, String)
    # Commande = dmidecode -s processor-version
    tailleMemoire = read(pipeline(`lshw -c memory`, pipeline(`grep taille:`, pipeline(`head -1`,  `awk   '{print $2}'`))), String)
    # Commande = lshw -c memory | grep taille: | head -1 | awk   '{print $2}'
    typeMemoire = read(pipeline(`lshw -c memory`, pipeline(`grep description:`, pipeline(`sed -n 2p`, `cut -c 24-`))), String)
    # Commande = lshw -c memory | grep description: | sed -n 2p | cut -c 24-
    l = "\n\n\\vspace{5mm}"
    l *= "\n\\section{Environment}"
    l *= "\n\n\\begin{itemize}"
    l *= "\n    \\item OS: \\textbf{"*os1*" - "*os2*"}"
    l *= "\n    \\item Manufacturer + Product: \\textbf{"*systeme*"}"
    l *= "\n    \\item CPU: \\textbf{"*processeur*"}"
    l *= "\n    \\item Memory: \\textbf{"*tailleMemoire*" "*typeMemoire*"}"
    l *= "\n    \\item Julia Version: \\textbf{"*string(VERSION)*"}"
    l *= "\n    \\item Code Verbosity: \\textbf{None}"
    l *= "\n\n\\end{itemize}"
    write(f, l)
    println(" --> OK")

    # Résultats en partant de Concorde TSP et avec AEP/EP
    l = "\n\n\\newpage"
    l *= "\n\n\\section{Results}"
    l *= "\n\\subsection{Optimal-TSP}"
    write(f, l)

    deb = 1
    premier = true
    nbTableaux = round(length(listeInstances)/32, RoundUp)
    println("\n Progression : ")
    println(" -------------")
    for i in 1:nbTableaux
        tableau(f, listeInstances, typeInstance, vitesseCamion, vitesseDrone, deb, premier, i, nbTableaux)
        deb += 32
        premier = false
    end

    # Résultats en partant de Concorde TSP vs. LKH TSP
    l = "\n\n\\newpage"
    l *= "\n\\subsection{Heuristic-TSP vs optimal-TSP}"
    write(f,l)

    deb = 1
    nbTableaux = round(length(listeInstances)/20, RoundUp)
    println("\n Progression (bis) : ")
    println(" ---------------------")
    for i in 1:nbTableaux
        tableauBis(f, listeInstances, typeInstance, vitesseCamion, vitesseDrone, deb, i, nbTableaux)
        deb += 20
        premier = false
    end

    ## Pied de page
    write(f, "\n\n\\end{document}")
    close(f)
    run(pipeline(`pdflatex rapport.tex`, stdout=devnull))
    #run(`pdflatex rapport.tex`) # Deuxième compilation pour la table des matières (nécessaire)
    run(`rm rapport.aux rapport.log`)
    run(`mv rapport.tex sorties/`)
    run(`mv rapport.pdf sorties/`)
    println("\n Le rapport est disponible dans le dossier \"sorties/\"")
end

function tableau(f, listeInstances, typeInstance, vitesseCamion, vitesseDrone, deb, premier, i, nbTableaux)
    ## En-tête tableau
    if premier == false # Tableaux suivants
        l = "\n\n\\newpage"
    else
        l = ""
    end
    l *= "\n\n\\renewcommand{\\arraystretch}{1.4}"
    l *= "\n\\begin{table}[H]"
    l *= "\n\\hspace{-10mm}\\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}"
    l *= "\n\\cline{2-10}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{9}{|c|}{\\large{Computational Results - Overview ("*string(trunc(Int, i))*"/"*string(trunc(Int, nbTableaux))*")}} \\\\"
    l *= "\n\\cline{2-10}"
    l *= "\n\\multicolumn{1}{c|}{}& \\multicolumn{3}{|c}{\\textbf{TSP}} & \\multicolumn{3}{|c}{\\textbf{AEP (A1)}} & \\multicolumn{3}{|c|}{\\textbf{EP (A2)}} \\\\"
    l *= "\n\\hline"
    l *= "\n\\textbf{Instance} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{Time (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{Time (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{Time (\$s\$)}  \\\\"
    l *= "\n\\hline"
    l *= "\n\\hline"
    write(f,l)

    # Corps
    compteur = 0 # Maximum 32 instances par page
    for i in deb:length(listeInstances)
        println(" --> "*string(i)*" : OK")
        if typeInstance == 1
            instance = parserAgatzBis(listeInstances[i])
        elseif typeInstance == 2
            instance = parserNousBis(listeInstances[i])
            # On échappe les "_" pour pouvoir les afficher ci-après
            listeInstances[i] = replace(listeInstances[i], "_"=>"\\_")
        elseif typeInstance == 3
            instance = parserTsplibBis(listeInstances[i])
        end
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone = vitesseDrone
        temps = @elapsed concordeFast(instance)
        tspSolution = instance.solution.ordreVisiteCamion
        tspTempsParcours = instance.solution.tempsParcours # Pour ne pas refaire après
        l = "\n\\textit{\""*listeInstances[i]*"\"} & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(temps)
        temps = @elapsed AEP(instance)
        l *= " & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(temps)*" "
        instance.solution = Solution(tspSolution, [], tspTempsParcours, 0) # Reset
        temps = @elapsed EP(instance)
        l *= "& "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(temps)*" \\\\"
        l *= "\n\\hline"
        write(f, l)
        compteur += 1
        if typeInstance == 2
            listeInstances[i] = replace(listeInstances[i], "\\_"=>"_")
        end
        if compteur == 32
            break
        end
    end

    ## Pied tableau
    l = "\n\\end{tabular}"
    l *= "\n\\end{table}"
    write(f, l)
end

function tableauBis(f, listeInstances, typeInstance, vitesseCamion, vitesseDrone, deb, i, nbTableaux)
    ## En-tête tableau
    l = "\n\n\\renewcommand{\\arraystretch}{1.4}"
    l *= "\n\\begin{table}[H]"
    l *= "\n\\rotatebox{90}{\\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|c|c|}"
    l *= "\n\\cline{2-13}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{12}{|c|}{\\large{Computational Results - Overview ("*string(trunc(Int, i))*"/"*string(trunc(Int, nbTableaux))*")}} \\\\"
    l *= "\n\\cline{2-13}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{6}{|c|}{\\textbf{TSP via Concorde}} & \\multicolumn{6}{|c|}{\\textbf{TSP via Linkern}}  \\\\"
    l *= "\n\\cline{2-13}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{3}{|c}{\\textbf{AEP (A1)}} & \\multicolumn{3}{|c|}{\\textbf{EP (A2)}} & \\multicolumn{3}{|c}{\\textbf{AEP (A1)}} & \\multicolumn{3}{|c|}{\\textbf{EP (A2)}} \\\\"
    l *= "\n\\hline"
    l *= "\n\\textbf{Instance} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} \\\\"
    l *= "\n\\hline"
    l *= "\n\\hline"
    write(f,l)

    # Corps
    compteur = 0 # Maximum 20 instances par page (format paysage)
    for i in deb:length(listeInstances)
        println(" --> "*string(i)*" : OK")
        if typeInstance == 1
            instance = parserAgatzBis(listeInstances[i])
        elseif typeInstance == 2
            instance = parserNousBis(listeInstances[i])
            # On échappe les "_" pour pouvoir les afficher ci-après
            listeInstances[i] = replace(listeInstances[i], "_"=>"\\_")
        elseif typeInstance == 3
            instance = parserTsplibBis(listeInstances[i])
        end
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone = vitesseDrone
        l = "\n\\textit{\""*listeInstances[i]*"\"}"
        concordeFast(instance)
        tspSolution = instance.solution.ordreVisiteCamion
        tspTempsParcours = instance.solution.tempsParcours # Pour ne pas refaire après
        temps = @elapsed AEP(instance)
        l *= " & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" "
        instance.solution = Solution(tspSolution, [], tspTempsParcours, 0) # Reset
        temps = @elapsed EP(instance)
        l *= "& "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" "
        # LKH
        instance.solution = Solution([], [], 0, 0)
        lkhFast(instance)
        tspSolution = instance.solution.ordreVisiteCamion
        tspTempsParcours = instance.solution.tempsParcours # Pour ne pas refaire après
        temps = @elapsed AEP(instance)
        l *= " & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" "
        instance.solution = Solution(tspSolution, [], tspTempsParcours, 0) # Reset
        temps = @elapsed EP(instance)
        l *= "& "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" \\\\"
        l *= "\n\\hline"
        write(f, l)
        compteur += 1
        if compteur == 20
            break
        end
    end

    ## Pied tableau
    l = "\n\\end{tabular}}"
    l *= "\n\\end{table}"
    write(f, l)
end

########## Poikonen ##########

function expPoiko()
    f = open("rapport.tex", "w")
    # En-tête
    l = "\\documentclass{article}"
    l *= "\n\\usepackage[utf8]{inputenc}"
    l *= "\n\\usepackage[top=2cm, bottom=3cm, left=2.5cm, right=2.5cm]{geometry}"
    l *= "\n\\usepackage{float}"
    l *= "\n\\usepackage{graphicx}"
    l *= "\n\n\\begin{document}"
    l *= "\n\n\\hrule"
    l *= "\n\\begin{center}"
    l *= "\n\\Large{\\textbf{Experimentations}}"
    l *= "\n\\vspace{0.1mm}"
    l *= "\n\\end{center}"
    l *= "\n\\hrule \\vspace{8mm}"
    write(f, l)

    # Table des matières
    #l = "\n\n\\tableofcontents"
    #write(f, l)

    # Logo
    l = "\n\n\\vspace{4mm}\\begin{figure}[H]"
    l *= "\n\\center"
    l *= "\n\\includegraphics[scale=0.22]{logo/logo+.png}"
    l *= "\n\\end{figure}"
    write(f, l)

    # Lecture instance
    bibliotheque = readdir("experimentations")
    lignes = CSV.read("experimentations/"*bibliotheque[1], header=false)
    println("\n Nombre d'instances : " * string(size(lignes,1)))
    print("\n --> Borne inférieure : ")
    borneInf = parse(Int64, readline())
    print(" --> Borne supérieure : ")
    borneSup = parse(Int64, readline())

    # Paramètres
    vitesseCamion = 1
    vitesseDrone = 3
    l = "\n\n\\vspace{5mm}"
    l *= "\n\\section{Parameters}"
    l *= "\n\n\\begin{itemize}"
    l *= "\n    \\item Truck speed: \\textbf{"*string(vitesseCamion)*"}"
    l *= "\n    \\item Drone speed: \\textbf{"*string(vitesseDrone)*"}"
    l *= "\n    \\item Number of instances: \\textbf{"*string(borneSup - borneInf + 1)*"}"
    l *= "\n\n\\end{itemize}"
    write(f, l)

    # Environnement
    println("\n Détection du matériel :")
    println(" -----------------------")
    os1 = read(`uname -s`, String) # Read permet de lire le résultat de la commande
    os2 = read(`uname -r`, String)
    # Commande = uname -s -r
    systeme = read(`dmidecode -s system-product-name`, String)
    # Commande = dmidecode -s system-product-name
    processeur = read(`dmidecode -s processor-version`, String)
    # Commande = dmidecode -s processor-version
    tailleMemoire = read(pipeline(`lshw -c memory`, pipeline(`grep taille:`, pipeline(`head -1`,  `awk   '{print $2}'`))), String)
    # Commande = lshw -c memory | grep taille: | head -1 | awk   '{print $2}'
    typeMemoire = read(pipeline(`lshw -c memory`, pipeline(`grep description:`, pipeline(`sed -n 2p`, `cut -c 24-`))), String)
    # Commande = lshw -c memory | grep description: | sed -n 2p | cut -c 24-
    l = "\n\n\\vspace{5mm}"
    l *= "\n\\section{Environment}"
    l *= "\n\n\\begin{itemize}"
    l *= "\n    \\item OS: \\textbf{"*os1*" - "*os2*"}"
    l *= "\n    \\item Manufacturer + Product: \\textbf{"*systeme*"}"
    l *= "\n    \\item CPU: \\textbf{"*processeur*"}"
    l *= "\n    \\item Memory: \\textbf{"*tailleMemoire*" "*typeMemoire*"}"
    l *= "\n    \\item Julia Version: \\textbf{"*string(VERSION)*"}"
    l *= "\n    \\item Code Verbosity: \\textbf{None}"
    l *= "\n\n\\end{itemize}"
    write(f, l)
    println(" --> OK")

    # Résultats en partant de Concorde TSP et avec AEP/EP
    l = "\n\n\\newpage"
    l *= "\n\n\\section{Results}"
    l *= "\n\\subsection{Optimal-TSP}"
    write(f, l)

    premier = true
    nbTableaux = round((borneSup - borneInf + 1)/32, RoundUp)
    println("\n Progression : ")
    println(" -------------")
    cpt = borneInf
    for i in 1:nbTableaux
        tableauPoiko(f, vitesseCamion, vitesseDrone, cpt, premier, i, nbTableaux, borneSup, lignes)
        cpt += 32
        premier = false
    end

    # Résultats en partant de Concorde TSP vs. LKH TSP
    l = "\n\n\\newpage"
    l *= "\n\\subsection{Heuristic-TSP vs optimal-TSP}"
    write(f,l)

    nbTableaux = round((borneSup - borneInf + 1)/20, RoundUp)
    println("\n Progression (bis) : ")
    println(" ---------------------")
    cpt = borneInf
    for i in 1:nbTableaux
        tableauPoikoBis(f, vitesseCamion, vitesseDrone, cpt, i, nbTableaux, borneSup, lignes)
        cpt += 20
        premier = false
    end

    ## Pied de page
    write(f, "\n\n\\end{document}")
    close(f)
    run(pipeline(`pdflatex rapport.tex`, stdout=devnull))
    #run(`pdflatex rapport.tex`) # Deuxième compilation pour la table des matières (nécessaire)
    run(`rm rapport.aux rapport.log`)
    run(`mv rapport.tex sorties/`)
    run(`mv rapport.pdf sorties/`)
    println("\n Le rapport est disponible dans le dossier \"sorties/\"")
end

function tableauPoiko(f, vitesseCamion, vitesseDrone, cpt, premier, i, nbTableaux, borneSup, lignes)
    ## En-tête tableau
    if premier == false # Tableaux suivants
        l = "\n\n\\newpage"
    else
        l = ""
    end
    l *= "\n\n\\renewcommand{\\arraystretch}{1.4}"
    l *= "\n\\begin{table}[H]\n\\begin{center}"
    l *= "\n\\hspace{-10mm}\\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}"
    l *= "\n\\cline{2-10}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{9}{|c|}{\\large{Computational Results - Overview ("*string(trunc(Int, i))*"/"*string(trunc(Int, nbTableaux))*")}} \\\\"
    l *= "\n\\cline{2-10}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{3}{|c}{\\textbf{TSP}} & \\multicolumn{3}{|c}{\\textbf{AEP (A1)}} & \\multicolumn{3}{|c|}{\\textbf{EP (A2)}} \\\\"
    l *= "\n\\hline"
    l *= "\n\\textbf{Instance} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{Time (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{Time (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{Time (\$s\$)}  \\\\"
    l *= "\n\\hline"
    l *= "\n\\hline"
    write(f,l)

    # Corps
    compteur = 0 # Maximum 32 instances par page
    for i in cpt:borneSup
        println(" --> "*string(i)*" : OK")
        instance = parserPoikoBis(lignes, i)
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone = vitesseDrone
        temps = @elapsed concordeFast(instance)
        l = "\n\\textit{\"instance\\_"*string(i)*"\"} & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(temps)
        temps = @elapsed AEP(instance)
        l *= " & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(temps)*" "
        instance.solution = Solution([], [], 0, 0) # Reset
        concordeFast(instance)
        temps = @elapsed EP(instance)
        l *= "& "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(temps)*" \\\\"
        l *= "\n\\hline"
        write(f, l)
        compteur += 1
        if compteur == 32
            break
        end
    end

    ## Pied tableau
    l = "\n\\end{tabular}\n\\end{center}"
    l *= "\n\\end{table}"
    write(f, l)
end

function tableauPoikoBis(f, vitesseCamion, vitesseDrone, cpt, i, nbTableaux, borneSup, lignes)
    ## En-tête tableau
    l = "\n\n\\renewcommand{\\arraystretch}{1.4}"
    l *= "\n\\begin{table}[H]"
    l *= "\n\\rotatebox{90}{\\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|c|c|}"
    l *= "\n\\cline{2-13}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{12}{|c|}{\\large{Computational Results - Overview ("*string(trunc(Int, i))*"/"*string(trunc(Int, nbTableaux))*")}} \\\\"
    l *= "\n\\cline{2-13}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{6}{|c|}{\\textbf{TSP via Concorde}} & \\multicolumn{6}{|c|}{\\textbf{TSP via Linkern}}  \\\\"
    l *= "\n\\cline{2-13}"
    l *= "\n\\multicolumn{1}{c|}{} & \\multicolumn{3}{|c}{\\textbf{AEP (A1)}} & \\multicolumn{3}{|c|}{\\textbf{EP (A2)}} & \\multicolumn{3}{|c}{\\textbf{AEP (A1)}} & \\multicolumn{3}{|c|}{\\textbf{EP (A2)}} \\\\"
    l *= "\n\\hline"
    l *= "\n\\textbf{Instance} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} & \\textbf{\$z_1\$} & \\textbf{\$z_2\$} & \\textbf{T (\$s\$)} \\\\"
    l *= "\n\\hline"
    l *= "\n\\hline"
    write(f,l)

    # Corps
    compteur = 0 # Maximum 20 instances par page (format paysage)
    for i in cpt:borneSup
        println(" --> "*string(i)*" : OK")
        instance = parserPoikoBis(lignes, i)
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone = vitesseDrone
        l = "\n\\textit{\"instance\\_"*string(i)*"\"}"
        concordeFast(instance)
        tspSolution = instance.solution.ordreVisiteCamion
        tspTempsParcours = instance.solution.tempsParcours # Pour ne pas refaire après
        temps = @elapsed AEP(instance)
        l *= " & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" "
        instance.solution = Solution(tspSolution, [], tspTempsParcours, 0) # Reset
        temps = @elapsed EP(instance)
        l *= "& "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" "
        # LKH
        instance.solution = Solution([], [], 0, 0)
        lkhFast(instance)
        tspSolution = instance.solution.ordreVisiteCamion
        tspTempsParcours = instance.solution.tempsParcours # Pour ne pas refaire après
        temps = @elapsed AEP(instance)
        l *= " & "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" "
        instance.solution = Solution(tspSolution, [], tspTempsParcours, 0) # Reset
        temps = @elapsed EP(instance)
        l *= "& "*string(round(instance.solution.tempsParcours,digits=2))*" & "*string(round(instance.solution.tempsAttente,digits=2))*" & "*string(round(temps, digits=7))*" \\\\"
        l *= "\n\\hline"
        write(f, l)
        compteur += 1
        if compteur == 20
            break
        end
    end

    ## Pied tableau
    l = "\n\\end{tabular}}"
    l *= "\n\\end{table}"
    write(f, l)
end

mainExp()

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver