text:
"""
```
solve!(m; algorithm = :relaxation, autodiff = :central, verbose = :high, kwargs...)
solve!(m, z0, y0; kwargs...)
solve!(m, z0, y0, Ψ0; kwargs...)
```

computes the risk-adjusted linearization of the dynamic economic model
described by `m` and updates `m` with the solution,
e.g. the coefficients ``(z, y, \\Psi)``.

The three available `solve!` algorithms are slight variations on each other.

- Method 1: uses the `z`, `y`, and `Ψ` fields of `m` as initial guesses
    for ``(z, y, \\Psi)`` and proceeds with the numerical algorithm
    specified by `algorithm`

- Method 2: uses `z0` and `y0` as initial guesses for the deterministic
    steady state, which is then used as the initial guess for ``(z, Y, \\Psi)``
    for the numerical algorithm specified by `algorithm`.

- Method 3: uses `z0`, `y0`, and `Ψ0` as initial guesses for ``(z, Y, \\Psi)``
    and proceeds with the numerical algorithm specified by `algorithm`.

### Inputs
- `m::RiskAdjustedLinearization`: object holding functions needed to calculate
    the risk-adjusted linearization
- `z0::AbstractVector{S1}`: initial guess for ``z``
- `y0::AbstractVector{S1}`: initial guess for ``y``
- `Ψ0::AbstractVector{S1}`: initial guess for ``\\Psi``
- `S1 <: Real`

### Keywords
- `algorithm::Symbol = :relaxation`: which numerical algorithm to use? Can be one of `[:relaxation, :homotopy, :deterministic]`
- `autodiff::Symbol = :central`: use autodiff or not? This keyword is the same as in `nlsolve`
- `use_anderson::Bool = false`: use Anderson acceleration if the relaxation algorithm is applied. Defaults to `false`
- `step::Float64 = .1`: size of step from 0 to 1 if the homotopy algorithm is applied. Defaults to 0.1
- `sparse_jacobian::Bool = false`: if true, exploit sparsity in the Jacobian in calls to `nlsolve` using SparseDiffTools.jl.
    If `jac_cache` and `sparsity` are `nothing`, then `solve!` will attempt to determine the sparsity pattern.
- `sparsity::Union{AbstractArray, Nothing} = nothing`: sparsity pattern for the Jacobian in calls to `nlsolve`
- `colorvec = nothing`: matrix coloring vector for sparse Jacobian in calls to `nlsolve`
- `jac_cache = nothing`: pre-allocated Jacobian cache for calls to `nlsolve` during the numerical algorithms
- `sparsity_detection::Bool = false`: If true, use SparsityDetection.jl to detect sparsity pattern (only relevant if
    both `jac_cache` and `sparsity` are `nothing`). If false,  then the sparsity pattern is
    determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero.
    Currently, SparsityDetection.jl fails to work.

The solution algorithms all use `nlsolve` to calculate the solution to systems of nonlinear
equations. The user can pass in any of the keyword arguments for `nlsolve` to adjust
the settings of the nonlinear solver.

For the keywords relevant to specific methods, see the docstring for the underlying method being called.
Note these methods are not exported.

- `:relaxation` -> `relaxation!`
- `:homotopy` -> `homotopy!`
- `:deterministic` -> `deterministic_steadystate!`
"""
function solve!(m::RiskAdjustedLinearization; algorithm::Symbol = :relaxation,
                autodiff::Symbol = :central, use_anderson::Bool = false,
                step::Float64 = .1, sparse_jacobian::Bool = false,
                sparsity::Union{AbstractArray, Nothing} = nothing, colorvec = nothing,
                jac_cache = nothing, sparsity_detection::Bool = false,
                verbose::Symbol = :high, kwargs...)
    if algorithm == :deterministic
        solve!(m, m.z, m.y; algorithm = algorithm, autodiff = autodiff,
               sparse_jacobian = sparse_jacobian,
               sparsity = sparsity, colorvec = colorvec,
               jac_cache = jac_cache, sparsity_detection = sparsity_detection,
               verbose = verbose, kwargs...)
    else
        solve!(m, m.z, m.y, m.Ψ; algorithm = algorithm, autodiff = autodiff,
               use_anderson = use_anderson, step = step,
               sparse_jacobian = sparse_jacobian, sparsity = sparsity,
               colorvec = colorvec, jac_cache = jac_cache,
               sparsity_detection = sparsity_detection, verbose = verbose, kwargs...)
    end
end

function solve!(m::RiskAdjustedLinearization, z0::AbstractVector{S1}, y0::AbstractVector{S1};
                algorithm::Symbol = :relaxation, autodiff::Symbol = :central,
                use_anderson::Bool = false, step::Float64 = .1,
                sparse_jacobian::Bool = false,
                sparsity::Union{AbstractArray, Nothing} = nothing, colorvec = nothing,
                jac_cache = nothing, sparsity_detection::Bool = false,
                verbose::Symbol = :high, kwargs...) where {S1 <: Real}

    @assert algorithm in [:deterministic, :relaxation, :homotopy] "The algorithm must be :deterministic, :relaxation, or :homotopy"

    # Deterministic steady state
    deterministic_steadystate!(m, vcat(z0, y0); autodiff = autodiff,
                               sparse_jacobian = sparse_jacobian,
                               sparsity = sparsity, colorvec = colorvec,
                               jac_cache = jac_cache,
                               sparsity_detection = sparsity_detection,
                               verbose = verbose, kwargs...)

    # Calculate linearization
    nl = nonlinear_system(m)
    li = linearized_system(m)
    update!(nl, m.z, m.y, m.Ψ; select = Symbol[:μ, :ξ])
    update!(li, m.z, m.y, m.Ψ; select = Symbol[:Γ₁, :Γ₂, :Γ₃, :Γ₄])

    # Back out Ψ
    m.Ψ .= compute_Ψ(m; zero_entropy_jacobian = true)

    # Use deterministic steady state as guess for stochastic steady state?
    if algorithm == :deterministic
        # Zero the entropy and Jacobian terms so they are not undefined or something else
        m.nonlinear[:𝒱_sss]  .= 0.
        m.linearization[:JV] .= 0.

        # Check Blanchard-Kahn
        blanchard_kahn(m; deterministic = true, verbose = verbose)
    else
        solve!(m, m.z, m.y, m.Ψ; algorithm = algorithm,
               use_anderson = use_anderson, step = step,
               sparse_jacobian = sparse_jacobian,
               jac_cache = jac_cache, sparsity = sparsity,
               colorvec = colorvec, sparsity_detection = sparsity_detection,
               verbose = verbose, kwargs...)
    end

    m
end

function solve!(m::RiskAdjustedLinearization, z0::AbstractVector{S1}, y0::AbstractVector{S1}, Ψ0::AbstractMatrix{S1};
                algorithm::Symbol = :relaxation, autodiff::Symbol = :central,
                use_anderson::Bool = false, step::Float64 = .1,
                sparse_jacobian::Bool = false, sparsity::Union{AbstractArray, Nothing} = nothing,
                colorvec = nothing, jac_cache = nothing, sparsity_detection::Bool = false,
                verbose::Symbol = :high, kwargs...) where {S1 <: Number}

    @assert algorithm in [:relaxation, :homotopy] "The algorithm must be :relaxation or :homotopy because this function calculates the stochastic steady state"

    # Stochastic steady state
    if algorithm == :relaxation
        N_zy = m.Nz + m.Ny
        relaxation!(m, vcat(z0, y0), Ψ0; autodiff = autodiff,
                    use_anderson = use_anderson, sparse_jacobian = sparse_jacobian,
                    sparsity = sparsity, colorvec = colorvec,
                    jac_cache = jac_cache, sparsity_detection = sparsity_detection,
                    verbose = verbose, kwargs...)
    elseif algorithm == :homotopy
        homotopy!(m, vcat(z0, y0, vec(Ψ0)); autodiff = autodiff, step = step,
                  sparse_jacobian = sparse_jacobian,
                  sparsity = sparsity, colorvec = colorvec,
                  jac_cache = jac_cache, sparsity_detection = sparsity_detection,
                  verbose = verbose, kwargs...)
    end

    # Check Blanchard-Kahn
    blanchard_kahn(m; deterministic = false, verbose = verbose)

    m
end

"""
```
function deterministic_steadystate!(m, x0; verbose = :none, kwargs...)
```

calculates the deterministic steady state.

### Types:
- `S1 <: Number`
- `S2 <: Real`

### Inputs
- `m::RiskAdjustedLinearization`: object holding functions needed to calculate
    the risk-adjusted linearization
- `x0::AbstractVector{S1}`: initial guess for ``(z, y)``
"""
function deterministic_steadystate!(m::RiskAdjustedLinearization, x0::AbstractVector{S1};
                                    autodiff::Symbol = :central,
                                    sparse_jacobian::Bool = false,
                                    sparsity::Union{AbstractArray, Nothing} = nothing, colorvec = nothing,
                                    jac_cache = nothing, sparsity_detection::Bool = false,
                                    verbose::Symbol = :none, kwargs...) where {S1 <: Real, S2 <: Real}

    # Set up system of equations
    _my_eqn = (F, x) -> _deterministic_equations(F, x, m)

    # Exploit sparsity?
    if sparse_jacobian
        nlsolve_jacobian!, jac =
            construct_sparse_jacobian_function(m, _my_eqn, :deterministic, autodiff;
                                               sparsity = sparsity, colorvec = colorvec,
                                               jac_cache = jac_cache,
                                               sparsity_detection = sparsity_detection)
        out = nlsolve(OnceDifferentiable(_my_eqn, nlsolve_jacobian!, x0, copy(x0), jac), x0; kwargs...)
    else
        out = nlsolve(OnceDifferentiable(_my_eqn, x0, copy(x0), autodiff,
                                         ForwardDiff.Chunk(min(m.Nz, m.Ny))), x0; kwargs...)
    end

    if out.f_converged
        m.z .= out.zero[1:m.Nz]
        m.y .= out.zero[(m.Nz + 1):end]

        if verbose in [:low, :high]
            println("A deterministic steady state has been found")
        end
    else
        error("A deterministic steady state could not be found.")
    end
end

function _deterministic_equations(F::AbstractVector{<: Number}, x::AbstractVector{<: Number},
                                  m::RiskAdjustedLinearization)
    # Unpack input vector
    z = @view x[1:m.Nz]
    y = @view x[(m.Nz + 1):end]

    # Update μ(z, y) and ξ(z, y)
    update!(m.nonlinear, z, y, m.Ψ; select = Symbol[:μ, :ξ])

    # Calculate residuals
    μ_sss             = get_tmp(m.nonlinear.μ.cache, z, y, (1, 1)) # select the first DiffCache b/c that
    ξ_sss             = get_tmp(m.nonlinear.ξ.cache, z, y, (1, 1)) # one corresponds to autodiffing both z and y
    F[1:m.Nz]         = μ_sss - z
    F[(m.Nz + 1):end] = ξ_sss + m.linearization[:Γ₅] * z + m.linearization[:Γ₆] * y
end

"""
```
steady_state_errors(m::RiskAdjustedLinearization, stochastic_steady_state::Bool = true)
```

calculates the errors in the system of equations characterizing the steady state.
When the second input is `true`, the steady state is the stochastic steady state,
and when it is false, the steady state is the deterministic steady state.
"""
function steady_state_errors(m::RiskAdjustedLinearization, stochastic_steady_state::Bool = true)
    if stochastic_steady_state
        return vcat(m[:μ_sss] - m.z, m[:ξ_sss] + m[:Γ₅] * m.z + m[:Γ₆] * m.y + m[:𝒱_sss],
                    vec(m[:Γ₃] + m[:Γ₄] * m.Ψ + (m[:Γ₅] + m[:Γ₆] * m.Ψ) * (m[:Γ₁] + m[:Γ₂] * m.Ψ) + m[:JV]))
    else
        return vcat(m[:μ_sss] - m.z, m[:ξ_sss] + m[:Γ₅] * m.z + m[:Γ₆] * m.y)
    end
end


text:
#=   FlexFloat Typeology    =#

abstract Art

# C is Clay  , the base or substructural concrete type
# Q is Qualia, the manner of enhancement as abstract type
# S is Sculpt, the art unfolding, qualia as dynamism guiding dispatch
#
abstract EnhancedFloat{C, S, Q} <: Real

#=
   These Sculpts are as intervals,
      bounded about lo and hi, covering amidst lo and hi.

   Each boundry either is Closed(Cl) or it is Open(Op).

      A Closed boundry does include its bounding value.
        A Closed boundry is considered complete at its bounding value.
      An Open boundry does not include its bounding value.
        An Open boundry is considered complete before its bounding value
        (above with an open lower bound, below with an open upper bound).

      Open upper bounds are treated as inclusively reaching from
         prevfloat(bounding value) outward, approching arbitrarily
         close to its bounding value.  Open bounds are given logic
         to support participation with the coinciding real value
         span [nextfloat(bounding value) for open lower bounds].
=#

abstract Sculpt <: Art

                         #   explicatively        traditionally
                         #
abstract CLCL <: Sculpt  #    [.......]             [.....]
abstract OPOP <: Sculpt  #     ).....(              (.....)
abstract OPCL <: Sculpt  #     )......]             (.....]
abstract CLOP <: Sculpt  #    [......(              [.....)

#=
   These Qualia are as bistable states,
      each boundry either is exact or it is inexact.
      If both boundries are exact, the interval is exact.
      If either boundry is inexact, the interval is inexact.
=#
abstract Qualia  <: Art

abstract EXACT   <: Qualia
abstract INEXACT <: Qualia

abstract Supple{S<:Sculpt, Q<:Qualia} <: Real

#= Clay is extensible with Union{} =#    typealias Clay AbstractFloat

immutable Flex{S, Q, C<:Clay} <: Supple{S, Q}
    lo::C
    hi::C
end

for (fn,Fn,S) in ((:ClCl,:clcl,:CLCL),(:ClOp,:clop,:CLOP),(:OpCl,:opcl,:OPCL),(:OpOp,:opop,:OPOP))
    @eval begin
        ($Fn){C<:Clay}(fp::C) = Flex{$S,EXACT,C}(fp,fp)
        convert{C<:Clay}(::Type{Flex{$S,EXACT,C}},x::C) = Flex{$S,EXACT,C}(x,x)
        convert{C<:Clay}(::Type{Flex{$S,EXACT,C}},lo::C,hi::C) = Flex{$S,EXACT,C}(minmax(lo,hi)...)
        convert{C<:Clay}(::Type{Flex{$S,EXACT,C}},x::Flex{$S,INEXACT,C}) = Flex{$S,EXACT,C}(x.lo,x.hi)

        ($fn){C<:Clay}(fp::C) = Flex{$S,INEXACT,C}(fp,fp)
        convert{C<:Clay}(::Type{Flex{$S,INEXACT,C}},x::C) = Flex{$S,INEXACT,C}(x,x)
        convert{C<:Clay}(::Type{Flex{$S,INEXACT,C}},lo::C,hi::C) = Flex{$S,INEXACT,C}(minmax(lo,hi)...)
        convert{C<:Clay}(::Type{Flex{$S,INEXACT,C}},x::Flex{$S,EXACT,C}) = Flex{$S,INEXACT,C}(x.lo,x.hi)
   end
end

for (Fn,fn,S) in ((:ClCl,:clcl,:CLCL),(:ClOp,:clop,:CLOP),(:OpCl,:opcl,:OPCL),(:OpOp,:opop,:OPOP))
    @eval begin
        ($fn){C<:Clay}(fp::C) = Flex{$S,EXACT,C}(fp,fp)
        ($Fn){C<:Clay}(fp::C) = Flex{$S,INEXACT,C}(fp,fp)
        function ($fn){C<:Clay}(lo::C,hi::C)
            lo,hi = minmax(lo,hi)
            Flex{$S,EXACT,C}(lo,hi)
        end
        function ($Fn){C<:Clay}(lo::C,hi::C)
            lo,hi = minmax(lo,hi)
            Flex{$S,INEXACT,C}(lo,hi)
        end
        function ($fn)(fp::Real)
            fq = convert(Float64,fp)
            Flex{$S,EXACT,Float64}(fq,fq)
        end
        function ($Fn)(fp::Real)
            fq = convert(Float64,fp)
            Flex{$S,INEXACT,Float64}(fq,fq)
        end
        function ($fn)(lo::Real,hi::Real)
            low = convert(Float64, lo)
            hig = convert(Float64, hi)
            low,hig = minmax(low,hig)
            Flex{$S,EXACT,Float64}(low,hig)
        end
        function ($fn){R<:Real}(lo::R,hi::R)
            low = convert(Float64, lo)
            hig = convert(Float64, hi)
            low,hig = minmax(low,hig)
            Flex{$S,EXACT,Float64}(low,hig)
        end
        function ($Fn){R<:Real}(lo::R,hi::R)
            low = convert(Float64, lo)
            hig = convert(Float64, hi)
            low,hig = minmax(low,hig)
            Flex{$S,INEXACT,Float64}(low,hig)
        end
        function ($fn){R<:Real,S<:Real}(lo::R,hi::S)
            low = convert(Float64, lo)
            hig = convert(Float64, hi)
            low,hig = minmax(low,hig)
            Flex{$S,EXACT,Float64}(low,hig)
        end
        function ($Fn){R<:Real,S<:Real}(lo::R,hi::S)
            low = convert(Float64, lo)
            hig = convert(Float64, hi)
            low,hig = minmax(low,hig)
            Flex{$S,INEXACT,Float64}(low,hig)
        end
    end
end

typealias Flex64 Flex{CLCL,EXACT,Float64}
typealias Flux64 Flex{OPOP,INEXACT,Float64}

for fn in (:clcl, :opcl, :clop, :opop, :ClCl, :OpCl, :ClOp, :OpOp)
    @eval begin
         @vectorize_1arg Real ($fn)
         @vectorize_2arg Real ($fn)
    end
end

ext: jl


text:
"Check assertions on fields of an in-memory MDF by replacing it temporarily."
function test_mdf_replacement(mdf::MDFv2InMemory, getter::Symbol, replacement, message::String)
  # First check if the file is consistent
  @test checkConsistency(mdf)

  # Replace the field
  f = getfield(MPIFiles, getter)
  temp = f(mdf)
  f(mdf, replacement)

  let err = nothing
    try
      checkConsistency(mdf)
    catch err

    end

    @test err isa AssertionError
    @test err.msg == message
  end

  # Change it back
  f(mdf, temp)
end

@testset "MDFInMemory" begin
  fnMeasBruker = joinpath(datadir,"BrukerStore","20150915_102110_Wuerfelphantom_1_1","18")
  fnSMBruker = joinpath(datadir,"BrukerStore","20141121_130749_CalibrationScans_1_1","76")
  fnSM1DBruker = joinpath(datadir,"BrukerStore","20170807_142514_Service_1_1","89")

  fnMeasV1 = joinpath(tmpdir,"mdf","measurement_V1.mdf")
  fnMeasV2 = joinpath(tmpdir,"mdf","measurement_V2.mdf")
  fnSMV1 = joinpath(tmpdir,"mdf","systemMatrix_V1.mdf")
  fnSMV2 = joinpath(tmpdir,"mdf","systemMatrix_V2.mdf")
  fnSMV3 = joinpath(tmpdir,"mdf","systemMatrix_V3.mdf")
  fnSMV4 = joinpath(tmpdir,"mdf","systemMatrix_V4.mdf")
  fnSMV5 = joinpath(tmpdir,"mdf","systemMatrix_V5.mdf")
  fnSMV6 = joinpath(tmpdir,"mdf","systemMatrix_V6.mdf")
  fnSM1DV1 = joinpath(tmpdir,"mdf","systemMatrix1D_V1.mdf")
  fnSM1DV2 = joinpath(tmpdir,"mdf","systemMatrix1D_V2.mdf")

  fnMeasV2_converted = joinpath(tmpdir,"mdf","measurement_V2_converted.mdf")
  fnSMV2_converted = joinpath(tmpdir,"mdf","systemMatrix_V2_converted.mdf")

  @testset "Conversion" begin
    @testset "Fields" begin
      measBruker = MPIFile(fnMeasBruker)
      @test typeof(measBruker) == BrukerFileMeas
      saveasMDF(fnMeasV2, measBruker)#, frames=1:100) <- TODO test this

      mdfv2 = MPIFile(fnMeasV2)
      @test typeof(mdfv2) <: MDFFileV2

      # Back and forth conversion
      mdfv2InMemory = MDFv2InMemory(mdfv2)
      @test typeof(mdfv2InMemory) <: MDFv2InMemory

      message = "There is an inconsistency in the given in-memory MDF. The message is: `Inconsistent dimension J in `gradient` in `acquisition`.`."
      @test_logs (:warn, message) saveasMDF(fnMeasV2_converted, mdfv2InMemory)

      mdf = MPIFile(fnMeasV2_converted)
      @test typeof(mdf) <: MDFFileV2

      @test length(mdf) == 1

      @test studyName(mdf) == "Wuerfelphantom"
      @test studyNumber(mdf) == 1
      @test studyDescription(mdf) == "n.a."
      @test studyTime(mdf) == DateTime( "2015-09-15T10:21:10.992" )
      @test studyUuid(mdf) == studyUuid(measBruker) #UUID("fe9635e0-7cfc-44eb-9f5a-585013a2cb51") see issue #32 on github
      @test uuid_version(studyUuid(mdf)) == 4

      @test experimentName(mdf) == "fuenf (E18)"
      @test experimentNumber(mdf) == 18
      @test experimentDescription(mdf) == "fuenf (E18)"
      @test experimentSubject(mdf) == "Wuerfelphantom"
      @test experimentIsSimulation(mdf) == false
      @test experimentIsCalibration(mdf) == false
      @test experimentUuid(mdf) == experimentUuid(measBruker) # UUID("ef2110b5-cce8-4ca6-b041-38ea01254c47")
      @test uuid_version(experimentUuid(mdf)) == 4

      @test scannerFacility(mdf) == "Universitätsklinikum Hamburg Eppendorf"
      @test scannerOperator(mdf) == "nmrsu"
      @test scannerManufacturer(mdf) == "Bruker/Philips"
      @test scannerName(mdf) == "Preclinical MPI System"
      @test scannerTopology(mdf) == "FFP"

      @test tracerName(mdf) == ["Resovist"]
      @test tracerBatch(mdf) == ["0"]
      @test tracerVendor(mdf) == ["n.a."]
      @test tracerVolume(mdf) == [0.0]
      @test tracerConcentration(mdf) == [0.5]
      @test tracerInjectionTime(mdf) == [DateTime("2015-09-15T11:17:23.011")]

      @test acqStartTime(mdf) == DateTime("2015-09-15T11:17:23.011")
      @test acqGradient(mdf)[:,:,1,1] == [-1.25 0 0; 0 -1.25 0;0 0 2.5]
      @test acqFramePeriod(mdf) == 6.528E-4
      @test acqNumPeriodsPerFrame(mdf) == 1
      @test acqOffsetFieldShift(mdf)[:,1,1] == [0.0; 0.0; -0.0]

      @test dfNumChannels(mdf) == 3
      @test dfWaveform(mdf) == fill("sine", (1, 1))
      @test dfStrength(mdf)[:,:,1] == [0.014 0.014 0.0]
      @test dfPhase(mdf)[:,:,1] == [1.5707963267948966 1.5707963267948966 1.5707963267948966]
      @test dfBaseFrequency(mdf) == 2500000.0
      @test dfDivider(mdf)[:,1] == [102; 96; 99]
      @test dfCycle(mdf) == 6.528E-4

      @test rxNumChannels(mdf) == 3
      @test rxBandwidth(mdf) == 1250000.0
      @test rxNumSamplingPoints(mdf) == 1632
      @test rxDataConversionFactor(mdf) == repeat([1.0, 0.0], outer=(1,rxNumChannels(mdf)))
      @test size(rxTransferFunction(mdf)) == (817, 3)
      @test occursin("../../../transferFunction/tf.h5", rxTransferFunctionFileName(mdf))

      @test acqNumAverages(mdf) == 1
      @test acqNumFrames(mdf) == 500
      @test acqNumPeriodsPerFrame(mdf) == 1
      @test acqNumPeriods(mdf) == 500
      @test acqNumPatches(mdf) == 1
      @test acqNumPeriodsPerPatch(mdf) == 1

      @test size( measData(mdf) ) == (1632,3,1,500)
      @test size( measDataTDPeriods(mdf) ) == (1632,3,500)
      @test size( measDataTDPeriods(mdf, 101:200) ) == (1632,3,100)

      # Utility functions
      @test measIsCalibProcessed(mdf) == false
      @test experimentHasReconstruction(mdf) == false
      @test experimentHasMeasurement(mdf) == true
      @test rxHasTransferFunction(mdf) == true
      @test experimentIsCalibration(mdf) == false

      N = acqNumFrames(mdf)

      @test size(getMeasurements(mdf, numAverages=1,
                  spectralLeakageCorrection=false)) == (1632,3,1,500)

      @test size(getMeasurements(mdf, numAverages=10,
                  spectralLeakageCorrection=false)) == (1632,3,1,50)

      @test size(getMeasurements(mdf, numAverages=10, frames=1:100,
                  spectralLeakageCorrection=true)) == (1632,3,1,10)

      @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100)) == (817,3,1,10)

      @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100, loadasreal=true)) == (1634,3,1,10)

      @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50)
    end

    @testset "Calibration" begin
      smBruker = MPIFile(fnSMBruker)
      @test typeof(smBruker) == BrukerFileCalib
      saveasMDF(fnSMV2, smBruker)

      smv2 = MPIFile(fnSMV2)
      @test typeof(smv2) <: MDFFileV2

      # Back and forth conversion
      smInMemory = MDFv2InMemory(smv2)
      @test typeof(smInMemory) <: MDFv2InMemory

      saveasMDF(fnSMV2_converted, smInMemory)

      sm = MPIFile(fnSMV2_converted)
      @test typeof(sm) <: MDFFileV2

      @test size( systemMatrixWithBG(sm) ) == (1959,817,3,1)
      @test size( systemMatrix(sm,1:10) ) == (1936,10)

      @test size(rxTransferFunction(sm)) == (817, 3)
      @test rxHasTransferFunction(sm) == true
      @test occursin("../../../transferFunction/tf.h5", rxTransferFunctionFileName(sm))

      @test measIsFourierTransformed(sm) == true
      @test measIsTFCorrected(sm) == false
      @test measIsFastFrameAxis(sm) == true
      @test measIsBGCorrected(sm) == false
      @test measFramePermutation(sm) == [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 178, 177, 176, 175, 174, 173, 172, 171, 170, 169, 168, 167, 166, 165, 164, 163, 162, 161, 160, 159, 158, 157, 156, 155, 154, 153, 152, 151, 150, 149, 148, 147, 146, 145, 144, 143, 142, 141, 140, 139, 138, 137, 136, 135, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 267, 266, 265, 264, 263, 262, 261, 260, 259, 258, 257, 256, 255, 254, 253, 252, 251, 250, 249, 248, 247, 246, 245, 244, 243, 242, 241, 240, 239, 238, 237, 236, 235, 234, 233, 232, 231, 230, 229, 228, 227, 226, 225, 224, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 356, 355, 354, 353, 352, 351, 350, 349, 348, 347, 346, 345, 344, 343, 342, 341, 340, 339, 338, 337, 336, 335, 334, 333, 332, 331, 330, 329, 328, 327, 326, 325, 324, 323, 322, 321, 320, 319, 318, 317, 316, 315, 314, 313, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 445, 444, 443, 442, 441, 440, 439, 438, 437, 436, 435, 434, 433, 432, 431, 430, 429, 428, 427, 426, 425, 424, 423, 422, 421, 420, 419, 418, 417, 416, 415, 414, 413, 412, 411, 410, 409, 408, 407, 406, 405, 404, 403, 402, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 534, 533, 532, 531, 530, 529, 528, 527, 526, 525, 524, 523, 522, 521, 520, 519, 518, 517, 516, 515, 514, 513, 512, 511, 510, 509, 508, 507, 506, 505, 504, 503, 502, 501, 500, 499, 498, 497, 496, 495, 494, 493, 492, 491, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 623, 622, 621, 620, 619, 618, 617, 616, 615, 614, 613, 612, 611, 610, 609, 608, 607, 606, 605, 604, 603, 602, 601, 600, 599, 598, 597, 596, 595, 594, 593, 592, 591, 590, 589, 588, 587, 586, 585, 584, 583, 582, 581, 580, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 712, 711, 710, 709, 708, 707, 706, 705, 704, 703, 702, 701, 700, 699, 698, 697, 696, 695, 694, 693, 692, 691, 690, 689, 688, 687, 686, 685, 684, 683, 682, 681, 680, 679, 678, 677, 676, 675, 674, 673, 672, 671, 670, 669, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 801, 800, 799, 798, 797, 796, 795, 794, 793, 792, 791, 790, 789, 788, 787, 786, 785, 784, 783, 782, 781, 780, 779, 778, 777, 776, 775, 774, 773, 772, 771, 770, 769, 768, 767, 766, 765, 764, 763, 762, 761, 760, 759, 758, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 890, 889, 888, 887, 886, 885, 884, 883, 882, 881, 880, 879, 878, 877, 876, 875, 874, 873, 872, 871, 870, 869, 868, 867, 866, 865, 864, 863, 862, 861, 860, 859, 858, 857, 856, 855, 854, 853, 852, 851, 850, 849, 848, 847, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 979, 978, 977, 976, 975, 974, 973, 972, 971, 970, 969, 968, 967, 966, 965, 964, 963, 962, 961, 960, 959, 958, 957, 956, 955, 954, 953, 952, 951, 950, 949, 948, 947, 946, 945, 944, 943, 942, 941, 940, 939, 938, 937, 936, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999, 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023, 1024, 1068, 1067, 1066, 1065, 1064, 1063, 1062, 1061, 1060, 1059, 1058, 1057, 1056, 1055, 1054, 1053, 1052, 1051, 1050, 1049, 1048, 1047, 1046, 1045, 1044, 1043, 1042, 1041, 1040, 1039, 1038, 1037, 1036, 1035, 1034, 1033, 1032, 1031, 1030, 1029, 1028, 1027, 1026, 1025, 1070, 1071, 1072, 1073, 1074, 1075, 1076, 1077, 1078, 1079, 1080, 1081, 1082, 1083, 1084, 1085, 1086, 1087, 1088, 1089, 1090, 1091, 1092, 1093, 1094, 1095, 1096, 1097, 1098, 1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108, 1109, 1110, 1111, 1112, 1113, 1157, 1156, 1155, 1154, 1153, 1152, 1151, 1150, 1149, 1148, 1147, 1146, 1145, 1144, 1143, 1142, 1141, 1140, 1139, 1138, 1137, 1136, 1135, 1134, 1133, 1132, 1131, 1130, 1129, 1128, 1127, 1126, 1125, 1124, 1123, 1122, 1121, 1120, 1119, 1118, 1117, 1116, 1115, 1114, 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170, 1171, 1172, 1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182, 1183, 1184, 1185, 1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194, 1195, 1196, 1197, 1198, 1199, 1200, 1201, 1202, 1246, 1245, 1244, 1243, 1242, 1241, 1240, 1239, 1238, 1237, 1236, 1235, 1234, 1233, 1232, 1231, 1230, 1229, 1228, 1227, 1226, 1225, 1224, 1223, 1222, 1221, 1220, 1219, 1218, 1217, 1216, 1215, 1214, 1213, 1212, 1211, 1210, 1209, 1208, 1207, 1206, 1205, 1204, 1203, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255, 1256, 1257, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266, 1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1276, 1277, 1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, 1288, 1289, 1290, 1291, 1335, 1334, 1333, 1332, 1331, 1330, 1329, 1328, 1327, 1326, 1325, 1324, 1323, 1322, 1321, 1320, 1319, 1318, 1317, 1316, 1315, 1314, 1313, 1312, 1311, 1310, 1309, 1308, 1307, 1306, 1305, 1304, 1303, 1302, 1301, 1300, 1299, 1298, 1297, 1296, 1295, 1294, 1293, 1292, 1337, 1338, 1339, 1340, 1341, 1342, 1343, 1344, 1345, 1346, 1347, 1348, 1349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, 1357, 1358, 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368, 1369, 1370, 1371, 1372, 1373, 1374, 1375, 1376, 1377, 1378, 1379, 1380, 1424, 1423, 1422, 1421, 1420, 1419, 1418, 1417, 1416, 1415, 1414, 1413, 1412, 1411, 1410, 1409, 1408, 1407, 1406, 1405, 1404, 1403, 1402, 1401, 1400, 1399, 1398, 1397, 1396, 1395, 1394, 1393, 1392, 1391, 1390, 1389, 1388, 1387, 1386, 1385, 1384, 1383, 1382, 1381, 1426, 1427, 1428, 1429, 1430, 1431, 1432, 1433, 1434, 1435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444, 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453, 1454, 1455, 1456, 1457, 1458, 1459, 1460, 1461, 1462, 1463, 1464, 1465, 1466, 1467, 1468, 1469, 1513, 1512, 1511, 1510, 1509, 1508, 1507, 1506, 1505, 1504, 1503, 1502, 1501, 1500, 1499, 1498, 1497, 1496, 1495, 1494, 1493, 1492, 1491, 1490, 1489, 1488, 1487, 1486, 1485, 1484, 1483, 1482, 1481, 1480, 1479, 1478, 1477, 1476, 1475, 1474, 1473, 1472, 1471, 1470, 1515, 1516, 1517, 1518, 1519, 1520, 1521, 1522, 1523, 1524, 1525, 1526, 1527, 1528, 1529, 1530, 1531, 1532, 1533, 1534, 1535, 1536, 1537, 1538, 1539, 1540, 1541, 1542, 1543, 1544, 1545, 1546, 1547, 1548, 1549, 1550, 1551, 1552, 1553, 1554, 1555, 1556, 1557, 1558, 1602, 1601, 1600, 1599, 1598, 1597, 1596, 1595, 1594, 1593, 1592, 1591, 1590, 1589, 1588, 1587, 1586, 1585, 1584, 1583, 1582, 1581, 1580, 1579, 1578, 1577, 1576, 1575, 1574, 1573, 1572, 1571, 1570, 1569, 1568, 1567, 1566, 1565, 1564, 1563, 1562, 1561, 1560, 1559, 1604, 1605, 1606, 1607, 1608, 1609, 1610, 1611, 1612, 1613, 1614, 1615, 1616, 1617, 1618, 1619, 1620, 1621, 1622, 1623, 1624, 1625, 1626, 1627, 1628, 1629, 1630, 1631, 1632, 1633, 1634, 1635, 1636, 1637, 1638, 1639, 1640, 1641, 1642, 1643, 1644, 1645, 1646, 1647, 1691, 1690, 1689, 1688, 1687, 1686, 1685, 1684, 1683, 1682, 1681, 1680, 1679, 1678, 1677, 1676, 1675, 1674, 1673, 1672, 1671, 1670, 1669, 1668, 1667, 1666, 1665, 1664, 1663, 1662, 1661, 1660, 1659, 1658, 1657, 1656, 1655, 1654, 1653, 1652, 1651, 1650, 1649, 1648, 1693, 1694, 1695, 1696, 1697, 1698, 1699, 1700, 1701, 1702, 1703, 1704, 1705, 1706, 1707, 1708, 1709, 1710, 1711, 1712, 1713, 1714, 1715, 1716, 1717, 1718, 1719, 1720, 1721, 1722, 1723, 1724, 1725, 1726, 1727, 1728, 1729, 1730, 1731, 1732, 1733, 1734, 1735, 1736, 1780, 1779, 1778, 1777, 1776, 1775, 1774, 1773, 1772, 1771, 1770, 1769, 1768, 1767, 1766, 1765, 1764, 1763, 1762, 1761, 1760, 1759, 1758, 1757, 1756, 1755, 1754, 1753, 1752, 1751, 1750, 1749, 1748, 1747, 1746, 1745, 1744, 1743, 1742, 1741, 1740, 1739, 1738, 1737, 1782, 1783, 1784, 1785, 1786, 1787, 1788, 1789, 1790, 1791, 1792, 1793, 1794, 1795, 1796, 1797, 1798, 1799, 1800, 1801, 1802, 1803, 1804, 1805, 1806, 1807, 1808, 1809, 1810, 1811, 1812, 1813, 1814, 1815, 1816, 1817, 1818, 1819, 1820, 1821, 1822, 1823, 1824, 1825, 1869, 1868, 1867, 1866, 1865, 1864, 1863, 1862, 1861, 1860, 1859, 1858, 1857, 1856, 1855, 1854, 1853, 1852, 1851, 1850, 1849, 1848, 1847, 1846, 1845, 1844, 1843, 1842, 1841, 1840, 1839, 1838, 1837, 1836, 1835, 1834, 1833, 1832, 1831, 1830, 1829, 1828, 1827, 1826, 1871, 1872, 1873, 1874, 1875, 1876, 1877, 1878, 1879, 1880, 1881, 1882, 1883, 1884, 1885, 1886, 1887, 1888, 1889, 1890, 1891, 1892, 1893, 1894, 1895, 1896, 1897, 1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 1914, 1958, 1957, 1956, 1955, 1954, 1953, 1952, 1951, 1950, 1949, 1948, 1947, 1946, 1945, 1944, 1943, 1942, 1941, 1940, 1939, 1938, 1937, 1936, 1935, 1934, 1933, 1932, 1931, 1930, 1929, 1928, 1927, 1926, 1925, 1924, 1923, 1922, 1921, 1920, 1919, 1918, 1917, 1916, 1915, 1, 90, 179, 268, 357, 446, 535, 624, 713, 802, 891, 980, 1069, 1158, 1247, 1336, 1425, 1514, 1603, 1692, 1781, 1870, 1959]

      @test size( calibSNR(sm) ) == (817,3,1)
      @test calibFov(sm) == [0.044; 0.044; 0.001]
      @test calibFovCenter(sm) == [0.0; -0.0; 0.0]
      @test calibSize(sm) == [44; 44; 1]
      @test calibOrder(sm) == "xyz"
      @test calibDeltaSampleSize(sm) == [0.0, 0.0, 0.0] #[0.001; 0.001; 0.001]
      @test calibMethod(sm) == "robot"

      # Utility functions
      @test measIsCalibProcessed(sm) == true
      @test experimentHasReconstruction(sm) == false
      @test experimentHasMeasurement(sm) == true
      @test rxHasTransferFunction(sm) == true
      @test experimentIsCalibration(sm) == true

      @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,)
      #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working

      @test size(getSystemMatrix(sm,1:10)) == (1936,10)
      @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20)
      @test size(getSystemMatrix(sm,1:10,bgCorrection=true)) == (1936,10)
      # test on the data level if the conversion was successfull
      SNRThresh = 2
      freq = filterFrequencies(smBruker,SNRThresh=SNRThresh)
      SBruker = getSystemMatrix(smBruker,frequencies=freq)
      S = getSystemMatrix(sm,frequencies=freq)

      relativeDeviation = zeros(Float32,length(freq))
      for f in 1:length(freq)
        relativeDeviation[f] = norm(SBruker[:,f]-S[:,f])/norm(SBruker[:,f])
      end
      # test if relative deviation for most of the frequency components is below 0.003
      @test quantile(relativeDeviation,0.95)<0.003
    end
  end

  @testset "Defaults" begin
    mdf = defaultMDFv2InMemory()

    @test typeof(time(mdf)) == DateTime
    @test typeof(uuid(mdf)) == UUID
    @test version(mdf) == VersionNumber("2.1.0")
    @test rxUnit(mdf) == "V"
  end

  @testset "Dict" begin
    @testset "String keys" begin
      dict = Dict{String, Any}()
      dict["dfBaseFrequency"] = 1600/40e3
      dict["experimentUuid"] = UUID("946a039e-48de-47ee-957e-a15af437e0be")
      dict["acqGradient"] = fill(0.0, (1, 1, 3, 3))
      dict["calibFov"] = fill(0.2, 3)

      mdf = inMemoryMDFFromDict(dict)

      @test dfBaseFrequency(mdf) == 1600/40e3
      @test experimentUuid(mdf) == UUID("946a039e-48de-47ee-957e-a15af437e0be")
      @test acqGradient(mdf) == fill(0.0, (1, 1, 3, 3))
      @test calibFov(mdf) == fill(0.2, 3)
      @test calibFieldOfView(mdf) == fill(0.2, 3)

      mdf = MDFv2InMemory(dict)

      @test dfBaseFrequency(mdf) == 1600/40e3
      @test experimentUuid(mdf) == UUID("946a039e-48de-47ee-957e-a15af437e0be")
      @test acqGradient(mdf) == fill(0.0, (1, 1, 3, 3))
      @test calibFov(mdf) == fill(0.2, 3)
      @test calibFieldOfView(mdf) == fill(0.2, 3)
    end

    @testset "Symbol keys" begin
      dict = Dict{Symbol, Any}()
      dict[:dfBaseFrequency] = 1600/40e3
      dict[:experimentUuid] = UUID("946a039e-48de-47ee-957e-a15af437e0be")
      dict[:acqGradient] = fill(0.0, (1, 1, 3, 3))
      dict[:calibFov] = fill(0.2, 3)

      mdf = inMemoryMDFFromDict(dict)

      @test dfBaseFrequency(mdf) == 1600/40e3
      @test experimentUuid(mdf) == UUID("946a039e-48de-47ee-957e-a15af437e0be")
      @test acqGradient(mdf) == fill(0.0, (1, 1, 3, 3))
      @test calibFov(mdf) == fill(0.2, 3)
      @test calibFieldOfView(mdf) == fill(0.2, 3)

      mdf = MDFv2InMemory(dict)

      @test dfBaseFrequency(mdf) == 1600/40e3
      @test experimentUuid(mdf) == UUID("946a039e-48de-47ee-957e-a15af437e0be")
      @test acqGradient(mdf) == fill(0.0, (1, 1, 3, 3))
      @test calibFov(mdf) == fill(0.2, 3)
      @test calibFieldOfView(mdf) == fill(0.2, 3)
    end

    @testset "MDF to Dict" begin
      mdf = defaultMDFv2InMemory()
      calibFov(mdf, fill(0.2, 3)) # One alias should be tested
      dict = inMemoryMDFToDict(mdf)

      @test typeof(dict["time"]) == DateTime
      @test typeof(dict["uuid"]) == UUID
      @test dict["version"] == VersionNumber("2.1.0")
      @test dict["rxUnit"] == "V"
      @test dict["calibFov"] == fill(0.2, 3)
    end
  end

  @testset "Consistency" begin
    A = 1 # tracer materials/injections for multi-color MPI
    N = 2 # acquired frames (N = O + E), same as a spatial position for calibration
    E = 1 # acquired background frames (E = N − O)
    O = 1 # acquired foreground frames (O = N − E)
    B = 1 # coefficients stored after sparsity transformation (B \\le O)
    J = 1 # periods within one frame
    Y = 1 # partitions of each patch position
    C = 1 # receive channels
    D = 1 # drive-field channels
    F = 1 # frequencies describing the drive-field waveform
    V = 800 # points sampled at receiver during one drive-field period
    W = V # sampling points containing processed data (W = V if no frequency selection or bandwidth reduction has been applied)
    K = Int(V/2 + 1)# frequencies describing the processed data (K = V/2 + 1 if no frequency selection or bandwidth reduction has been applied)
    Q = 1 # frames in the reconstructed MPI data set
    P = 1 # voxels in the reconstructed MPI data set
    S = 1 # channels in the reconstructed MPI data set

    mdf = MDFv2InMemory()
    mdf.root = MDFv2Root(
      time=DateTime("2021-04-26T17:12:21.686"),
      uuid=UUID("946a039e-48de-47ee-957e-a15af437e0be"),
      version=VersionNumber("2.1.0")
    )
    mdf.study = MDFv2Study(
      description = "n.a.",
      name = "n.a.",
      number = 1,
      time = DateTime("2021-04-26T17:12:21.686"),
      uuid = UUID("f2f3ae66-2fc3-49f7-91f9-71c1c2dc15e7")
    )
    mdf.experiment = MDFv2Experiment(
      description = "n.a.",
      isSimulation = true,
      name = "n.a.",
      number = 1,
      subject = "n.a.",
      uuid = UUID("3076d4bc-a2ef-46ef-8a99-dcd047379b8a")
    )
    mdf.tracer = MDFv2Tracer(
      batch = fill("n.a.", A),
      concentration = fill(1.0, A),
      injectionTime = fill(DateTime("2021-04-26T17:12:21.686"), A),
      name = fill("MyAwesomeTracer", A),
      solute = fill("n.a.", A),
      vendor = fill("Me", A),
      volume = fill(1.0, A),
    )
    mdf.scanner = MDFv2Scanner(
      boreSize = 0.3,
      facility = "MyAwesomeInstitute",
      manufacturer = "MeMyselfAndI",
      name = "MyAwesomeScanner",
      operator = "JustMe",
      topology = "FFL" # Of course ;)
    )
    drivefield = MDFv2Drivefield(
      baseFrequency = 40e3,
      cycle = 1600/40e3,
      divider = fill(1600, (D, F)),
      numChannels = D,
      phase = fill(0.0, (J, D, F)),
      strength = fill(1.0, (J, D, F)),
      waveform = fill("sine", (D, F)),
    )
    receiver = MDFv2Receiver(
      bandwidth = 20e3,
      dataConversionFactor = repeat([1/2^16 0], C),
      inductionFactor = fill(1.0, C),
      numChannels = C,
      numSamplingPoints = V,
      transferFunction = fill(1+0.5im, (C, K)),
      unit = "V"
    )
    mdf.acquisition = MDFv2Acquisition(
      gradient = fill(0.0, (J, Y, 3, 3)),
      numAverages = 1,
      numFrames = N,
      numPeriodsPerFrame = J,
      offsetField = fill(0.0, (J, Y, 3)),
      startTime = DateTime("2021-04-26T17:12:21.686"),
      drivefield = drivefield,
      receiver = receiver
    )
    mdf.measurement = MDFv2Measurement(;
      data = fill(0, (N, J, C, K)),
      framePermutation = fill(0, N),
      frequencySelection = collect(1:K),
      isBackgroundCorrected = true,
      isBackgroundFrame = vcat(fill(true, E), fill(false, O)),
      isFastFrameAxis = true,
      isFourierTransformed = true,
      isFramePermutation = true,
      isFrequencySelection = true,
      isSparsityTransformed = true,
      isSpectralLeakageCorrected = true,
      isTransferFunctionCorrected = true,
      sparsityTransformation = "DCT-I",
      subsamplingIndices = fill(0, (J, C, K, B))
    )
    mdf.calibration = MDFv2Calibration(
      deltaSampleSize = fill(0.001, 3),
      fieldOfView = fill(0.2, 3),
      fieldOfViewCenter = fill(0.0, 3),
      method = "robot",
      offsetFields = fill(0.0, (O, 3)),
      order = "xyz",
      positions = fill(0.0, (O, 3)),
      size = [1, 1, 1],
      snr = fill(0.0, (J, C, K))
    )
    mdf.reconstruction = MDFv2Reconstruction(
      data = fill(0, (Q, P, S)),
      fieldOfView = fill(0.2, 3),
      fieldOfViewCenter = fill(0.0, 3),
      isOverscanRegion = fill(false, P),
      order = "xyz",
      positions = fill(0.0, (P, 3)),
      size = [1, 1, 1]
    )

    # TODO: As soon as there is a @test_nothrow available, change this.
    # See also: https://github.com/JuliaLang/julia/issues/18780
    @test checkConsistency(mdf) # This call has a side-effect on the variables!

    # Check if non-required fields make consistency checking fail
    tracer = mdf.tracer
    mdf.tracer = nothing
    @test checkConsistency(mdf)
    mdf.tracer = tracer

    @testset "Derived variables" begin
      @test mdf.variables.A == A
      @test mdf.variables.N == N
      @test mdf.variables.E == E
      @test mdf.variables.O == O
      @test mdf.variables.B == B
      @test mdf.variables.J == J
      @test mdf.variables.Y == Y
      @test mdf.variables.C == C
      @test mdf.variables.D == D
      @test mdf.variables.F == F
      @test mdf.variables.V == V
      @test mdf.variables.W == W
      @test mdf.variables.K == K
      @test mdf.variables.Q == Q
      @test mdf.variables.P == P
      @test mdf.variables.S == S
    end

    @testset "Missing" begin
      # Structs
      test_mdf_replacement(mdf, :root, missing, "The field `root` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :study, missing, "The field `study` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :experiment, missing, "The field `experiment` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :scanner, missing, "The field `scanner` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :acquisition, missing, "The field `acquisition` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :drivefield, missing, "The field `drivefield` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :receiver, missing, "The field `receiver` is missing in the given in-memory MDF.")

      # Root
      test_mdf_replacement(mdf, :time, missing, "The field `time` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :uuid, missing, "The field `uuid` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :version, missing, "The field `version` is missing in the given in-memory MDF.")

      # Study
      test_mdf_replacement(mdf, :studyDescription, missing, "The field `description` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :studyName, missing, "The field `name` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :studyNumber, missing, "The field `number` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :studyUuid, missing, "The field `uuid` is missing in the given in-memory MDF.")

      # Experiment
      test_mdf_replacement(mdf, :experimentDescription, missing, "The field `description` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :experimentIsSimulation, missing, "The field `isSimulation` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :experimentName, missing, "The field `name` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :experimentNumber, missing, "The field `number` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :experimentSubject, missing, "The field `subject` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :experimentUuid, missing, "The field `uuid` is missing in the given in-memory MDF.")

      # Tracer
      test_mdf_replacement(mdf, :tracerBatch, missing, "The field `batch` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :tracerConcentration, missing, "The field `concentration` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :tracerName, missing, "The field `name` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :tracerSolute, missing, "The field `solute` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :tracerVendor, missing, "The field `vendor` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :tracerVolume, missing, "The field `volume` is missing in the given in-memory MDF.")

      # Scanner
      test_mdf_replacement(mdf, :scannerFacility, missing, "The field `facility` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :scannerManufacturer, missing, "The field `manufacturer` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :scannerName, missing, "The field `name` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :scannerOperator, missing, "The field `operator` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :scannerTopology, missing, "The field `topology` is missing in the given in-memory MDF.")

      # Drivefield
      test_mdf_replacement(mdf, :dfBaseFrequency, missing, "The field `baseFrequency` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :dfCycle, missing, "The field `cycle` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :dfDivider, missing, "The field `divider` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :dfNumChannels, missing, "The field `numChannels` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :dfPhase, missing, "The field `phase` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :dfStrength, missing, "The field `strength` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :dfWaveform, missing, "The field `waveform` is missing in the given in-memory MDF.")

      # Receiver
      test_mdf_replacement(mdf, :rxBandwidth, missing, "The field `bandwidth` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :rxNumChannels, missing, "The field `numChannels` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :rxNumSamplingPoints, missing, "The field `numSamplingPoints` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :rxUnit, missing, "The field `unit` is missing in the given in-memory MDF.")

      # Acquisition
      test_mdf_replacement(mdf, :acqNumAverages, missing, "The field `numAverages` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :acqNumFrames, missing, "The field `numFrames` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :acqNumPeriodsPerFrame, missing, "The field `numPeriodsPerFrame` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :acqStartTime, missing, "The field `startTime` is missing in the given in-memory MDF.")

      # Measurement
      #test_mdf_replacement(mdf, :measData, missing, "The field `data` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsBackgroundCorrected, missing, "The field `isBackgroundCorrected` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsBackgroundFrame, missing, "The field `isBackgroundFrame` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsFastFrameAxis, missing, "The field `isFastFrameAxis` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsFourierTransformed, missing, "The field `isFourierTransformed` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsFramePermutation, missing, "The field `isFramePermutation` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsFrequencySelection, missing, "The field `isFrequencySelection` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsSparsityTransformed, missing, "The field `isSparsityTransformed` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsSpectralLeakageCorrected, missing, "The field `isSpectralLeakageCorrected` is missing in the given in-memory MDF.")
      test_mdf_replacement(mdf, :measIsTransferFunctionCorrected, missing, "The field `isTransferFunctionCorrected` is missing in the given in-memory MDF.")

      # Calibration
      test_mdf_replacement(mdf, :calibMethod, missing, "The field `method` is missing in the given in-memory MDF.")

      # Reconstruction
      test_mdf_replacement(mdf, :recoData, missing, "The field `data` is missing in the given in-memory MDF.")
    end

    @testset "Sizes" begin
      # Tracer
      test_mdf_replacement(mdf, :tracerBatch, fill("n.a.", A+1), "Inconsistent dimensions in `batch` in `tracer`.")
      test_mdf_replacement(mdf, :tracerConcentration, fill(1.0, A+1), "Inconsistent dimensions in `concentration` in `tracer`.")
      test_mdf_replacement(mdf, :tracerInjectionTime, fill(DateTime("2021-04-26T17:12:21.686"), A+1), "Inconsistent dimensions in `injectionTime` in `tracer`.")
      test_mdf_replacement(mdf, :tracerName, fill("MyAwesomeTracer", A+1), "Inconsistent dimensions in `name` in `tracer`.")
      test_mdf_replacement(mdf, :tracerSolute, fill("n.a.", A+1), "Inconsistent dimensions in `solute` in `tracer`.")
      test_mdf_replacement(mdf, :tracerVendor, fill("Me", A+1), "Inconsistent dimensions in `vendor` in `tracer`.") # Works because A has been determined earlier
      test_mdf_replacement(mdf, :tracerVolume, fill(1.0, A+1), "Inconsistent dimensions in `volume` in `tracer`.")

      # Drivefield
      test_mdf_replacement(mdf, :dfDivider, fill(1600, (D+1, F)), "Inconsistent dimension D in `divider` in `drivefield`.")
      test_mdf_replacement(mdf, :dfDivider, fill(1600, (D, F+1)), "Inconsistent dimension F in `divider` in `drivefield`.") # Works because F has been determined earlier
      test_mdf_replacement(mdf, :dfPhase, fill(0.0, (J+1, D, F)), "Inconsistent dimension J in `phase` in `drivefield`.") # Works because J has been determined earlier
      test_mdf_replacement(mdf, :dfPhase, fill(0.0, (J, D+1, F)), "Inconsistent dimension D in `phase` in `drivefield`.")
      test_mdf_replacement(mdf, :dfPhase, fill(0.0, (J, D, F+1)), "Inconsistent dimension F in `phase` in `drivefield`.")
      test_mdf_replacement(mdf, :dfStrength, fill(1.0, (J+1, D, F)), "Inconsistent dimension J in `strength` in `drivefield`.")
      test_mdf_replacement(mdf, :dfStrength, fill(1.0, (J, D+1, F)), "Inconsistent dimension D in `strength` in `drivefield`.")
      test_mdf_replacement(mdf, :dfStrength, fill(1.0, (J, D, F+1)), "Inconsistent dimension F in `strength` in `drivefield`.")
      test_mdf_replacement(mdf, :dfWaveform, fill("sine", (D+1, F)), "Inconsistent dimension D in `waveform` in `drivefield`.")
      test_mdf_replacement(mdf, :dfWaveform, fill("sine", (D, F+1)), "Inconsistent dimension F in `waveform` in `drivefield`.")

      # Receiver
      test_mdf_replacement(mdf, :rxDataConversionFactor, repeat([1/2^16 0], C+1), "Inconsistent dimension C in `dataConversionFactor` in `receiver`.")
      test_mdf_replacement(mdf, :rxDataConversionFactor, repeat([1/2^16 0 0], C), "Inconsistent second dimension in `dataConversionFactor` in `receiver`.")
      test_mdf_replacement(mdf, :rxInductionFactor, fill(1.0, C+1), "Inconsistent dimension C in `inductionFactor` in `receiver`.")
      test_mdf_replacement(mdf, :rxTransferFunction, fill(1+0.5im, (C+1, K)), "Inconsistent dimension C in `transferFunction` in `receiver`.")
      test_mdf_replacement(mdf, :rxTransferFunction, fill(1+0.5im, (C, K+1)), "Inconsistent dimension K in `transferFunction` in `receiver`.")

      # Measurement
      test_mdf_replacement(mdf, :measFramePermutation, fill(0, N+1), "Inconsistent dimension N in `framePermutation` in `measurement`.")
      test_mdf_replacement(mdf, :measFrequencySelection, collect(1:K+1), "Inconsistent dimension K in `frequencySelection` in `measurement`.")
      test_mdf_replacement(mdf, :measIsBackgroundFrame, vcat(fill(true, E+1), fill(false, O)), "Inconsistent dimension N in `isBackgroundFrame` in `measurement`.")
      test_mdf_replacement(mdf, :measSparsityTransformation, nothing, "Field `sparsityTransformation` must be set when `isSparsityTransformed` is set in in `measurement`.")

      # Calibration
      test_mdf_replacement(mdf, :calibDeltaSampleSize, fill(0.001, 4), "Inconsistent length in `deltaSampleSize` in `calibration`.")
      test_mdf_replacement(mdf, :calibFieldOfView, fill(0.2, 4), "Inconsistent length in `fieldOfView` in `calibration`.")
      test_mdf_replacement(mdf, :calibFieldOfViewCenter, fill(0.0, 4), "Inconsistent length in `fieldOfViewCenter` in `calibration`.")
      test_mdf_replacement(mdf, :calibOffsetFields, fill(0.0, (O+1, 3)), "Inconsistent dimension O in `offsetFields` in `calibration`.")
      test_mdf_replacement(mdf, :calibOffsetFields, fill(0.0, (O, 4)), "Inconsistent second dimension in `offsetFields` in `calibration`.")
      test_mdf_replacement(mdf, :calibOrder, "bla", "Wrong `order` of `bla` in `calibration`.")
      test_mdf_replacement(mdf, :calibPositions, fill(0.0, (O+1, 3)), "Inconsistent dimension O in `positions` in `calibration`.")
      test_mdf_replacement(mdf, :calibPositions, fill(0.0, (O, 4)), "Inconsistent second dimension in `positions` in `calibration`.")
      test_mdf_replacement(mdf, :calibSize, [1, 1, 1, 1], "Inconsistent length in `size` in `calibration`.")
      test_mdf_replacement(mdf, :calibSize, [1, 1, 2], "The product of `size` with `[1, 1, 2]` must equal O.")
      test_mdf_replacement(mdf, :calibSnr, fill(0.0, (J+1, C, K)), "Inconsistent dimension J in `snr` in `calibration`.")
      test_mdf_replacement(mdf, :calibSnr, fill(0.0, (J, C+1, K)), "Inconsistent dimension C in `snr` in `calibration`.")
      test_mdf_replacement(mdf, :calibSnr, fill(0.0, (J, C, K+1)), "Inconsistent dimension K in `snr` in `calibration`.")

      # Reconstruction
      test_mdf_replacement(mdf, :recoFieldOfView, fill(0.2, 4), "Inconsistent length in `fieldOfView` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoFieldOfViewCenter, fill(0.0, 4), "Inconsistent length in `fieldOfViewCenter` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoIsOverscanRegion, fill(false, P+1), "Inconsistent length in `isOverscanRegion` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoOrder, "bla", "Wrong `order` of `bla` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoPositions, fill(0.0, (P+1, 3)), "Inconsistent dimension P in `positions` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoPositions, fill(0.0, (P, 4)), "Inconsistent second dimension in `positions` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoSize, [1, 1, 1, 1], "Inconsistent length in `size` in `reconstruction`.")
      test_mdf_replacement(mdf, :recoSize, [1, 1, 2], "The product of `size` with `[1, 1, 2]` must equal P.")
    end
  end
end


text:
#
# This file is the beginning a compilation of useful aerodynamic calculations
#
# August 20, 2021 		Initial code
# September 11, 2021		Added horizontal and vertical stabilizer volume ratios
# September 12, 2021		Restructured code to add the "Extract_data function"
# October 1, 2021		Data Conversions, AR, Sw and MAC functions added
# October 7, 2021		AC calculations
# October 22, 2021		Added conversion to Newtons and stall speed calculation
# October 23, 2021		Added Cg, wing loading and wing load
#
# This code is licensed under MIT license (see LICENSE.txt for details)

	#-------------------------- Data Conversion Dictionaries ---------------------------#

# convert arbitrary units to meters
To_meters 	 = Dict("in" => 0.0254, "ft" => 0.3048, "yd" => 0.9144, "cm" => 0.01)
To_sq_meters 	 = Dict("in^2" => 0.00064516, "ft^2" => 0.092903, "yd^2" => 0.836127, "cm^2" => 0.0001)
To_N		 = Dict("kg" => 9.80665, "lbs" => 4.448221615)
To_kg		 = Dict("lb" => 0.453592, "N" => 0.1019716213)
To_kgm		 = Dict("in lb" => 0.01152125)

	#--------------------------------- Extract Data ------------------------------------#
	# This is a helper function to pull data from the input matrix, load it into the
	# working structure and check the data for presnece.
	# This function also attempts to convert all units to metric

function Extract_data(df, var_array)
	# load the data into the structure
	df_rows = size(df)

	for i = 1:df_rows[1]												# read all lines up to the maximum row count
		for j in eachindex(var_array)
			if var_array[j].name == df[i,1]									# look for first value
				var_array[j].value = df[i,2]								# load data from spreadsheet into structure

				if df[i,3] != var_array[j].units							# check the units to make sure they match

					# this series of if statements could be collapsed, but the result was difficult to read
					if (var_array[j].units == "m")
						if (get(To_meters, df[i,3], 0) != 0)					# check to see if scale factor exists
							corrected_value = df[i,2] * To_meters[df[i,3]]
							@printf("\n>>> Converted %2.2f %s to %2.2f m", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "m"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned
						else
							print("\nError: Could not convert unit type = \"", df[i,3], "\"\n")
							exit(86)							# variable units are incorrect, stop the program
						end

					elseif (var_array[j].units == "m^2")
						if (get(To_sq_meters, df[i,3], 0) != 0)					# check to see if scale factor exists
							corrected_value = df[i,2] * To_sq_meters[df[i,3]]
							@printf("\n>>> Converted %2.2f %s to %2.2f m^2", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "m^2"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned
						else
							print("\nError: Could not convert unit type = \"", df[i,3], "\"\n")
							exit(86)							# variable units are incorrect, stop the program
						end

					elseif (var_array[j].units == "N")
						if (get(To_N, df[i,3], 0) != 0)						# check to see if scale factor exists
							corrected_value = df[i,2] * To_N[df[i,3]]
							@printf("\n>>> Converted %2.2f %s to %2.2f N", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "N"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned
						else
							print("\nError: Could not convert unit type = \"", df[i,3], "\"\n")
							exit(86)							# variable units are incorrect, stop the program
						end

					elseif (var_array[j].units == "kg")
						if (get(To_kg, df[i,3], 0) != 0)					# check to see if scale factor exists
							corrected_value = df[i,2] * To_kg[df[i,3]]
							@printf("\n>>> Converted %2.2f %s to %2.2f kg", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "kg"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned
						else
							print("\nError: Could not convert unit type = \"", df[i,3], "\"\n")
							exit(86)							# variable units are incorrect, stop the program
						end

					elseif (var_array[j].units == "kgm")
						if (get(To_kgm, df[i,3], 0) != 0)					# check to see if scale factor exists
							corrected_value = df[i,2] * To_kgm[df[i,3]]
							@printf("\n>>> Converted %2.2f %s to %2.2f kgm", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "kgm"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned
						else
							print("\nError: Could not convert unit type = \"", df[i,3], "\"\n")
							exit(86)							# variable units are incorrect, stop the program
						end

					# note this is a special case since it is more than just multiplying by a coefficient
					elseif (var_array[j].units == "C")

						if (df[i,3] == "F")							# check to see if scale factor exists
							corrected_value = (df[i,2] - 32) * (5/9)
							@printf("\n>>> Converted %2.2f %s to %2.2f C", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "C"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned

						elseif (df[i,3] == "K")
							corrected_value = df[i,2] - 273.15
							@printf("\n>>> Converted %2.2f %s to %2.2f C", df[i,2], df[i,3], corrected_value)
							df[i,2] = corrected_value					# update the data structure with the updated value
							df[i,3] = "C"							# update the data structure with the updated units
							var_array[j].value = df[i,2]					# update the data being returned

#						elseif (df[i,3] == "C")
#							print("\n Celsius\n")

						else
							print("\nError: Could not convert unit type = \"", df[i,3], "\"\n")
							exit(86)							# variable units are incorrect, stop the program
						end

					else
						print("\nInternal Error: ", var_array[j].name, " unknown units!!\n")	# report unit mismatch error

						exit(86)								# variable units are incorrect, stop the program
					end
				end

			var_array[j].stored = 1										# data was loaded into this variable
			end
		end
	end

	# Check for missing variables in the return data
	for j in eachindex(var_array)
		if var_array[j].stored != 1
			print("\nError: ", var_array[j].name, " value is missing!\nExiting\n")
			exit(86)
		end
	end

	return  				# note the array was passed in so no additional return pointer is needed
end

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

	#---------------------------------- Wing Loading (Wload) -------------------------------------#

function Wload(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	W  		= req_vars("W",  	0, "kg", 	0)	# W = weight - a vector defined as the product of mass and acceleration
	Sw	   	= req_vars("Sw",	0, "m^2",	0)	# Surface area reference wing

	var_array = [W, Sw]						# Array of structures to enable easy iteration

	Extract_data(df, var_array)					# Loads the array of structures with data, checks units and presence

	return (W.value / Sw.value)
end

	#----------------------------------- Wing Load  (WL) ----------------------------------------#

function WL(df)

	# define structures to hold definitions for the call: (name, value, units, stored)
	m  	   = req_vars("m",  	0, "kg", 	0)		# m = mass
	Cg	   = req_vars("Cg",	0, "m",		0)		# Cg location from nose
	lwdle  	   = req_vars("lwdle",  0, "m", 	0)		# Distance datum (nose) to wing LE
	lhdle  	   = req_vars("lhdle",  0, "m", 	0)		# Distance datum (nose) to horizontal stab LE
	AC	   = req_vars("AC",	0, "m",		0)		# Distance from wing LE to AC
	ACh	   = req_vars("ACh",	0, "m",		0)		# Distance from horizontal stab LE to AC

	var_array = [m, Cg, lwdle, lhdle, AC, ACh]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)					# Loads the array of structures with data, checks units and presence

	moment_1 = ((lwdle.value + AC.value) - Cg.value) * m.value
	tail_down_force = moment_1 / ((lhdle.value + ACh.value) - (lwdle.value + AC.value))

	return ((m.value + tail_down_force), tail_down_force)
end

	#-------------------------------- Center of Gravity (Cg) ------------------------------------#

function Cg(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	m  	= req_vars("m",  	0, "kg", 	0)		# m = mass
	M  	= req_vars("M",  	0, "kgm", 	0)		# M = moment kg meters

	var_array = [m, M]						# Array of structures to enable easy iteration

	Extract_data(df, var_array)					# Loads the array of structures with data, checks units and presence

	return (M.value / m.value)
end

	#--------------------- Horizontal Stabilizer Stall Speed (Vhstall) ------------------------#

function Vhstall(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Wh  	= req_vars("Wh",  	0,  "N", 	0)		# W = weight - a vector defined as the product of mass and acceleration (newtons)
	Clh  	= req_vars("Clh",  	0,  "kg/m^3",	0)		# Coefficient of lift
	Sht	= req_vars("Sht",	0,  "m^2",	0)		# Stabilizer surface area

	var_array = [Wh, Clh, Sht]					# Array of structures to enable easy iteration

	Extract_data(df, var_array)					# Loads the array of structures with data, checks units and presence

	ρ = 1.225							# freestream density of air - mass/unit volume of air upstream of a body

	intermediate_value = (2 * Wh.value) /  (ρ * Sht.value * Clh.value)

	return (sqrt(intermediate_value))
end

	#----------------------------- Wing Stall Speed (Vstall) ---------------------------------#

function Vstall(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	W  	= req_vars("W",  	0, "N", 	0)	# W = weight - a vector defined as the product of mass and acceleration (newtons)
	Cl  	= req_vars("Cl",  	0, "kg/m^3", 	0)	# Coefficient of lift
	Sw   	= req_vars("Sw",	0, "m^2",	0)	# Surface area reference wing

	var_array = [W, Cl, Sw]					# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and presence

	ρ = 1.225						# freestream density of air - mass/unit volume of air upstream of a body

	intermediate_value = (2 * W.value) /  (ρ * Sw.value * Cl.value)

	return (sqrt(intermediate_value))
end

	#--------------- Distance between AC of Wing and Vertical Stab (Lvt) -------------------#

function Lvt(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	lwdle	= req_vars("lwdle",	0, "m", 0)		# Distance datum (nose) to wing LE
	lvdle	= req_vars("lvdle",	0, "m", 0)		# Distance datum (nose) to vertical stab LE
	AC	= req_vars("AC",	0, "m",	0)		# Distance from wing LE to AC
	ACv	= req_vars("ACv",	0, "m",	0)		# Distance from horizontal stab LE to AC

	var_array = [lwdle, lvdle, AC, ACv]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and presence

	return (abs((lvdle.value + ACv.value) - (lwdle.value + AC.value)))
end

	#--------------- Distance between AC of Wing and Horizontal Stab (Lht) -------------------#

function Lht(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	lwdle  	= req_vars("lwdle",  	0, "m", 0)		# Distance datum (nose) to wing LE
	lhdle  	= req_vars("lhdle",  	0, "m", 0)		# Distance datum (nose) to horizontal stab LE
	AC	= req_vars("AC",	0, "m",	0)		# Distance from wing LE to AC
	ACh	= req_vars("ACh",	0, "m",	0)		# Distance from horizontal stab LE to AC

	var_array = [lwdle, lhdle, AC, ACh]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	return (abs((lhdle.value + ACh.value) - (lwdle.value + AC.value)))
end

	#----------------- Vertical Stab Aerodynamic Center Calculation (ACv) --------------------#

function ACv(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	λv	  = req_vars("λv",    0, "degrees", 	0)	# Vertical stab sweep angle
	bv	  = req_vars("bv",    0, "m",   	0)	# Vertical stab span
	cMACv	  = req_vars("cMACv", 0, "m",   	0)	# Vertical stab chord
	Crv	  = req_vars("Crv",   0, "m",   	0)	# Vertical stab root chord
	Ctv	  = req_vars("Ctv",   0, "m",   	0)	# Vertical stab tip chord

	var_array = [λv, bv, cMACv, Crv, Ctv]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	S = tan(deg2rad(λv.value))*bv.value/2						# Wing sweep distance from 0 sweep LE
	C = (S * (Crv.value + (2 * Ctv.value))) / (3 * (Crv.value + Ctv.value))		# Chord setback from 0 sweep LE

	return (C + (cMACv.value * 0.25))
end

	#------------ Vertical Stab Mean Aerodynamic Chord Calculation (cMACv) ----------------#

function MACv(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Ctv  = req_vars("Ctv",  0, "m", 0)			# Vertical stab tip chord
	Crv  = req_vars("Crv",  0, "m",   0)			# Vertical stab root chord

	var_array = [Ctv, Crv]					# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	taper_ratio = Ctv.value / Crv.value			# Taper ratio
	return (Crv.value*(2/3)*((1 + taper_ratio + taper_ratio^2)/(1+taper_ratio)))
end

	#------------ Horizontal Stab Mean Aerodynamic Chord Calculation (cMACh) ----------------#

function MACh(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Cth  = req_vars("Cth",  0, "m",   0)			# Horizontal stab tip chord
	Crh  = req_vars("Crh",  0, "m",   0)			# Horizontal stab root chord

	var_array = [Cth, Crh]					# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	taper_ratio = Cth.value / Crh.value			# Taper ratio
	return (Crh.value*(2/3)*((1 + taper_ratio + taper_ratio^2)/(1+taper_ratio)))
end

	#------------- Horizontal Stabilizer Aerodynamic Center Calculation (ACh) ---------------#

function ACh(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	λh	  = req_vars("λh",    0, "degrees", 	0)		# Horizontal Stab sweep angle
	bh	  = req_vars("bh",    0, "m",		0)		# Horizontal Stab span
	cMACh	  = req_vars("cMACh", 0, "m",   	0)		# Mean horizontal stab chord
	Crh	  = req_vars("Crh",   0, "m",		0)		# Horizontal stab root chord
	Cth	  = req_vars("Cth",   0, "m",   	0)		# Horizontal stab tip chord

	var_array = [λh, bh, cMACh, Crh, Cth]				# Array of structures to enable easy iteration

	Extract_data(df, var_array)					# Loads the array of structures with data, checks units and prenence

	S = tan(deg2rad(λh.value))*bh.value/2						# Wing sweep distance from 0 sweep LE
	C = (S * (Crh.value + (2 * Cth.value))) / (3 * (Crh.value + Cth.value))		# Chord setback from 0 sweep LE

	return (C + (cMACh.value * 0.25))
end

	#------------------------ Wing MAC distance from Chord (d) ----------------------------#

function d(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	b	 = req_vars("b",   0, "m",   0)			# Wing span
	Cr	 = req_vars("Cr",  0, "m",   0)			# Wing root chord
	Ct	 = req_vars("Ct",  0, "m",   0)			# Wing tip chord

	var_array = [b, Cr, Ct]					# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	return ((b.value*((0.5*Cr.value)+Ct.value)) / (3 * (Cr.value+Ct.value)))
end

	#---------------------- Wing Aerodynamic Center Calculation (AC) ----------------------#

function AC(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	λ    	 = req_vars("λ",    0, "degrees", 0)		# Wing sweep angle
	b	 = req_vars("b",    0, "m",   	0)		# Wing span
	cMAC 	 = req_vars("cMAC", 0, "m",   	0)		# Mean wing chord
	Cr	 = req_vars("Cr",   0, "m",   	0)		# Wing root chord
	Ct	 = req_vars("Ct",   0, "m",   	0)		# Wing tip chord

	var_array = [λ, b, cMAC, Cr, Ct]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	S = tan(deg2rad(λ.value))*b.value/2						# Wing sweep distance from 0 sweep LE
	C = (S * (Cr.value + (2 * Ct.value))) / (3 * (Cr.value + Ct.value))		# Chord setback from 0 sweep LE

	return (C + (cMAC.value * 0.25))
end

	#------------------- Vertical Stab Aspect Ratio Calculation (AR) ---------------------#

function ARv(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Svt 	= req_vars("Svt",  0, "m^2", 0)			# Vertical stab area
	bv	= req_vars("bv",   0, "m",   0)			# Vertical stab span
	Nv	= req_vars("Nv",   0, "count",  0)		# Number of vertical stabs

	var_array = [Nv, Svt, bv]				# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	return (((bv.value^2)*Nv.value)/Svt.value)
end

	#----------------------- Vertical Stab Area Calculation (Svt) ----------------------#

function Svt(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	bv	= req_vars("bv",   0, "m",   	0)		# Wing span
	Crv	= req_vars("Crv",  0, "m",   	0)		# Wing root chord
	Ctv	= req_vars("Ctv",  0, "m",  	0)		# Wing tip chord
	Nv	= req_vars("Nv",   0, "count",  0)		# Number of vertical stabs

	var_array = [Nv, bv, Crv, Ctv]				# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	average_chord = (Crv.value + Ctv.value)/2

	return (Nv.value*bv.value*average_chord)
end

	#------------------- Horizontal Stab Aspect Ratio Calculation (AR) ---------------------#

function ARh(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Sht = req_vars("Sht",  0, "m^2", 0)		# Wing area
	bh	= req_vars("bh",   0, "m",   0)		# Wing span

	var_array = [Sht, bh]				# Array of structures to enable easy iteration

	Extract_data(df, var_array)			# Loads the array of structures with data, checks units and prenence

	return ((bh.value^2)/Sht.value)
end

	#----------------------- Horizontal Stab Area Calculation (Swh) ----------------------#

function Sht(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	bh	= req_vars("bh",   0, "m",   0)		# Wing span
	Crh	= req_vars("Crh",  0, "m",   0)		# Wing root chord
	Cth	= req_vars("Cth",  0, "m",   0)		# Wing tip chord

	var_array = [bh, Crh, Cth]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)			# Loads the array of structures with data, checks units and prenence

	average_chord = (Crh.value + Cth.value)/2

	return (bh.value*average_chord)
end

	#----------------------------- Wing Area Calculation (Sw) ----------------------------#

function Sw(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	b	= req_vars("b",   0, "m",   0)		# Wing span
	Cr	= req_vars("Cr",  0, "m",   0)		# Wing root chord
	Ct	= req_vars("Ct",  0, "m",   0)		# Wing tip chord

	var_array = [b, Cr, Ct]				# Array of structures to enable easy iteration

	Extract_data(df, var_array)			# Loads the array of structures with data, checks units and prenence

	average_chord = (Cr.value + Ct.value)/2

	return (b.value*average_chord)
end

	#--------------------------- Aspect Ratio Calculation (AR) --------------------------#

function AR(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Sw  = req_vars("Sw",  0, "m^2", 0)		# Wing area
	b	= req_vars("b",   0, "m",   0)		# Wing span

	var_array = [Sw, b]				# Array of structures to enable easy iteration

	Extract_data(df, var_array)			# Loads the array of structures with data, checks units and prenence

	return ((b.value^2)/Sw.value)
end

	#------------------------ Mean Aerodynamic Chord Calculation -----------------------#

function MAC(df)
	# define structures to hold definitions for the call: (name, value, units, stored)
	Ct  = req_vars("Ct",  0, "m",   0)			# Wing tip chord
	Cr  = req_vars("Cr",  0, "m",   0)			# Wing root chord

	var_array = [Ct, Cr]					# Array of structures to enable easy iteration

	Extract_data(df, var_array)				# Loads the array of structures with data, checks units and prenence

	taper_ratio = Ct.value / Cr.value						# Wing taper ratio
	return (Cr.value*(2/3)*((1 + taper_ratio + taper_ratio^2)/(1+taper_ratio)))
end

	#--------------------------- Reynolds Number Calculation ---------------------------#
	# air density and dynamic_viscosity should be changed to calculations for improved accuracy
	# The numbers included below will get the calculations in the right order of magnitude for STP

function Re(velocity, chord)		# velocity = m/sec, chord = m

	air_density = 1.138				# kg/m^3
	dynamic_viscosity = 1.73E-05	# m/sec^2

	return air_density*velocity*chord/dynamic_viscosity
end

	#--------------------------- Horizontal Tail Volume Ratio ---------------------------#
function Cht(df)			# df = data matrix from excel file

	# define structures to hold definitions for the call: (name, value, units, stored)
	Sht  = req_vars("Sht",  0, "m^2", 0)		# Horizontal tail area
	Lht  = req_vars("Lht",  0, "m",   0)		# Length between the aerodynamic centers of the wing and horizontal tailplane
	Sw   = req_vars("Sw",   0, "m^2", 0)		# Wing area
	cMAC = req_vars("cMAC", 0, "m",   0)		# Mean aerodynamic chord

	var_array = [Sht, Lht, Sw, cMAC]		# Array of structures to enable easy iteration

	Extract_data(df, var_array)			# Loads the array of structures with data, checks units and prenence

	return ((Sht.value*Lht.value)/(Sw.value*cMAC.value))
end

	#--------------------------- Vertical Tail Volume Ratio ---------------------------#
function Cvt(df)			# df = data matrix from excel file

	# define structures to hold definitions for the call: (name, value, units, stored)
	Svt  = req_vars("Svt",  0, "m^2", 0)		# Vertical tail area
	Lvt  = req_vars("Lvt",  0, "m",   0)		# Length between the aerodynamic centers of the wing and vertical tailplane
	Sw   = req_vars("Sw",   0, "m^2", 0)		# Wing area
	b    = req_vars("b",    0, "m",   0)		# Wing span

	var_array = [Svt, Lvt, Sw, b]			# Array of structures to enable easy iteration

	Extract_data(df, var_array)			# Loads the array of structures with data, checks units and prenence

	return ((Svt.value*Lvt.value)/(Sw.value*b.value))
end


text:
export project, circle_position

"""
    project(x::Number, x1::Number, y1::Number, x2::Number, y2::Number)

Calculate for a `x` from the range `[x1,y1]` its value in the range `[x2,y2]`.

This includes `x`s outside of the input range and reverse.
See the Examples section below.

# Examples

```jldoctest
julia> project(5.0,4.0,6.0,3.0,8.0)
5.5
julia> project(2.0,0.0,1.0,0.0,1.0)
2.0
julia> project(0.3,1.0,0.0,0.0,1.0)
0.7
```
"""
@inline function project(x::Number, x1::Number, y1::Number, x2::Number, y2::Number)
    return (x-x1)/(y1-x1)*(y2-x2)+x2
end

"""
    circle_position(t::Int; radius::Real=1.0, center::Vector{Real}=[0.0, 0.0], number_of_points::Int=20)

Devides a circle in `t` equidistant sections and gives for every t between 0 and `number_of_points` - 1
the beginning postion of the secion on the circle beginning on the right most point on the circle.
With increasing `t` the position moves counter clockwise.

# Examples

```jldoctest
julia> circle_position(10)
[-1.0, 0.0]
```
"""
function circle_position(t::I;
                         radius::Number=1.0, center::AbstractVector{<:Number}=[0.0, 0.0], number_of_points::I=20) where I<:Integer
    x = t/number_of_points * 2π
    radius.*[cos(x), sin(x)] .+ center
end

"""
    function swap!(vec::AbstractVector, i::I, j::I) where I<:Integer

Fast swap of the elements `i` with the element `j` in a Vector `vec`.
WARNING: No bounds check will be applied for performance reasons.
"""
@inline function swap!(vec::AbstractVector, i::I, j::I) where I<:Integer
    @inbounds vec[i], vec[j] = vec[j], vec[i]
    return
end


text:
include("model/get_entities_by_color.jl")
include("model/safe_add_physical_group.jl")
include("model/safe_fragment.jl")
include("model/add_cad_names_to_physical_groups.jl")
include("model/add_materials_to_physical_groups_by_color.jl")
include("model/color_material_physical_group_entities.jl")
include("model/import_model.jl")
include("model/overlay_mpact_grid_hierarchy.jl")

include("mesh/set_mesh_field_using_materials.jl")
include("mesh/generate_mesh.jl")
include("mesh/get_cad_to_mesh_error.jl")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
######################################################################
# fov3.jl
# Like FOV, but has 3 groups instead of 2
######################################################################

type FOV3 <: Sensor
	region_probs::Vector{NTuple{2,Float64}}

	# for discretized measurements
	num_bins::Int
	bin_range::UnitRange{Int64}

	#FOV3(region_probs::Vector{NTuple{2,Float64}}) = new(region_probs, 2, 0:1)

    FOV3() = new([(60.0,.9), (120.0,0.9), (180.0,0.9)], 3, 1:3)
end

# 1 means it is in the field of view
# 0 means it is not
function observe(m::SearchDomain, s::FOV3, p::Pose)

	# ensure bearing is reflected across
	rel_bearing = fit_180(p[3] - true_bearing(p, m.theta))
	if rel_bearing < 0.0
		rel_bearing = -1.0 * rel_bearing
	end

    # find the region in which this rel bearing lies
	prob_in_view = 0.0
	n = length(s.region_probs)
    expected_obs = 1
	for i = 1:n
		temp_angle, temp_prob = s.region_probs[i]
		if rel_bearing  <= temp_angle
            expected_obs = i
			prob_in_view = temp_prob
			break
		end
	end

    # TODO: make this less inefficient
    obs_vals = [1,2,3]
    o_probs  = 0.5*(1.0-0.9) * ones(3)

    o_probs[expected_obs] = 0.9 #prob_in_view

    r = rand()
    o = 1
    if r < o_probs[1]
        o = obs_vals[1]
    elseif r < o_probs[2]
        o = obs_vals[2]
    else
        o = obs_vals[3]
    end

    return o * 1.0
end

# o technically has to be a float, but we only have two obs:
#  1 = in field of view
#  0 = not in field of view
function O(s::FOV3, theta::LocTuple, p::Pose, o::Float64)

	# determine relative bearing and fix it in 0 to 180
	rel_bearing = fit_180(p[3] - true_bearing(p, theta))
	if rel_bearing < 0.0
		rel_bearing = -1.0 * rel_bearing
	end

	prob_in_view = 0.0
	n = length(s.region_probs)
    expected_obs = 1
	for i = 1:n
		temp_angle, temp_prob = s.region_probs[i]
		if rel_bearing  <= temp_angle
			prob_in_view = temp_prob
            expected_obs = i * 1.0
			break
		end
	end

    ret_val = (o == expected_obs) ? 0.9 : 0.05

	return ret_val
end

# I don't really get why this is here...
#function O(x::Vehicle,s::FOV3,xp::Pose, theta::LocTuple, o::ObsBin, df::DF)
function O(s::FOV3, theta::LocTuple, xp::Pose, o::ObsBin)
	# determine relative bearing and fix it in 0 to 180
	rel_bearing = fit_180(xp[3] - true_bearing(xp, theta))
	if rel_bearing < 0.0
		rel_bearing = -1.0 * rel_bearing
	end

	prob_in_view = 0.0
	n = length(s.region_probs)
    expected_obs = 1
	for i = 1:n
		temp_angle, temp_prob = s.region_probs[i]
		if rel_bearing  <= temp_angle
			prob_in_view = temp_prob
            expected_obs = i * 1.0
			break
		end
	end

    ret_val = (o == expected_obs) ? 0.9 : 0.05
	return ret_val
end

obs2bin(o::Float64, s::FOV3) = round(Int, o)


text:
function horner(coefs, x)
    s = coefs[end]
    for k in length(coefs)-1:-1:1
        s = coefs[k] + x * s
    end
    return s
end

@show horner([-19, 7, -4, 6], 3)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    OILMM(fs, H)
An Orthogonal Instantaneous Linear Mixing Model (OILMM) -- a distribution over vector-
valued functions. Let `p` be the number of observed outputs, and `m` the number of latent
processes, then `H`, also known as the mixing matrix, is a `p x m` orthogonal matrix whose
column space spans the output space. The latent processes are represented an `Independent_MOGP`
as the latent processes remain decoupled.

# Arguments:
- f: a length-`m` vector of Gaussian process objects as an IndependentMOGP.
- H: a `p x m` orthogonal matrix representing a fixed basis of our p-dim target: h_1,...,h_m
"""
const OILMM = ILMM{<:IndependentMOGP, <:Orthogonal}

"""
    project(H, σ²)

Computes the projection `T` and `ΣT` given the mixing matrix and noise.
"""
function project(H::Orthogonal{Z}, σ²::Z) where {Z<:Real}
    U, S = H.U, H.S

    # Compute the projection of the data.
    T = sqrt(S) \ U'

    # Compute the projected noise, which is a matrix of size `size(Yproj)`.
    ΣT = diag(σ² * inv(S))

    return T, ΣT
end

"""
    rand(rng::AbstractRNG, fx::FiniteGP{<:OILMM})
Sample from the OILMM, including the observation noise.
Follows generative structure of model 2 from [1].
Follows the AbstractGPs.jl API.
See also `rand_latent`.
[1] - Bruinsma et al 2020.
"""
function AbstractGPs.rand(rng::AbstractRNG, fx::FiniteGP{<:OILMM})
    fs, H, σ², x = unpack(fx)

    # Obtain U and S
    U, S = H.U, H.S

    # Generate from the latent processes.
    X = reshape(reduce(vcat, map(f -> rand(rng, f(x)), fs.fs)), length(x), :)

    # Transform latents into observed space.
    F = vec((U * sqrt(S) * X')')

    # Generate iid noise and add to each output.
    return F .+ sqrt(noise_var(fx.Σy)) .* randn(rng, size(F))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.mean_and_var(fx::FiniteGP{<:OILMM})
    fs, H, σ², x = unpack(fx)

    # Compute the marginals over the independent latents.
    fs_marginals = reduce(hcat, map(f -> AbstractGPs.marginals(f(x)), fs.fs))
    M_latent = mean.(fs_marginals)'
    V_latent = var.(fs_marginals)'

    # Obtain U and S
    U, S = H.U, H.S

    # Compute the means.
    M = U * sqrt(S) * M_latent

    # Compute the variances.
    V = abs2.(U * sqrt(S)) * (V_latent) .+ σ²

    # Package everything into independent Normal distributions.
    return collect(vec(M')), collect(vec(V'))
end

# See AbstractGPs.jl API docs.
function AbstractGPs.logpdf(fx::FiniteGP{<:OILMM}, y::AbstractVector{<:Real})
    fs, H, σ², x = unpack(fx)

    # Projection step.
    Y = reshape_y(y, length(x))
    T, ΣT = project(H, σ²)
    Ty = T*Y

    # Latent process log marginal likelihood calculation.
    y_rows = collect(eachrow(Ty))
    ΣT_rows = collect(eachrow(repeat(ΣT, 1, size(Ty, 2))))
    lmls_latents = map((f, s, y) -> logpdf(f(x, collect(s)), collect(y)), fs.fs, ΣT_rows, y_rows)

    return sum(lmls_latents) + regulariser(H, σ², Y)
end

"""
    regulariser(fx, y)

Computes the regularisation term of the logpdf.
See e.g. appendix A.4 of [1] - Bruinsma et al 2020.
"""
function regulariser(
    H::Orthogonal{T},
    σ²::T,
    Y::AbstractMatrix{T},
) where {T<:Real}
    U, S = H.U, H.S

    n = size(Y, 2)
    p, m = size(U)

    return -(n * (logdet(S) + (p - m) * log(2π * σ²)) +
        sum(abs2, (I - U * U') * Y) / σ²) / 2
end

# See AbstractGPs.jl API docs.
function AbstractGPs.posterior(fx::FiniteGP{<:OILMM}, y::AbstractVector{<:Real})
    fs, H, σ², x = unpack(fx)

    # Projection step.
    Y = reshape_y(y, length(x))
    T, ΣT = project(H, σ²)
    Ty = T*Y

    # Condition each latent process on the projected observations.
    y_rows = collect(eachrow(Ty))
    ΣT_rows = collect(eachrow(repeat(ΣT, 1, size(Ty, 2))))
    fs_posterior = map(
        (f, s, y) -> AbstractGPs.posterior(f(x, collect(s)), collect(y)),
        fs.fs,
        ΣT_rows,
        y_rows,
    )
    return ILMM(independent_mogp(fs_posterior), H)
end


text:
# Dogleg is appropriate here because Jx'*Jx is positive definite, and in that
# case we only need to calculate one newton step. In the secular equation version
# we need repeaed factorizations, and that is not as easy to exploit (but maybe
# there's a good shifted one out there?)
struct NormedResiduals{Tx, Tfx, Tf}
    x::Tx
    Fx::Tfx
    F::Tf
end
has_batched_f(::NormedResiduals) = false
has_param(::NormedResiduals) = false
function value(nr::NormedResiduals, x)
    Fx = nr.F.F(nr.Fx, x)

    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
    nr.x .= x
    nr.Fx .= Fx

    f = (norm(Fx)^2)/2
    return f
end

function upto_gradient(nr::NormedResiduals, Fx, x)
    Fx, Jx = nr.F.FJ(nr.Fx, nr.Fx*x', x)

    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
    nr.x .= x
    nr.Fx .= Fx

    f = (norm(Fx)^2)/2
    Fx .= Jx'*Fx
   return f, Fx
end
function upto_hessian(nr::NormedResiduals, Fx, Jx, x)  #Fx is the gradient and Jx is the Hessian
    Fx, Jx = nr.F.FJ(nr.Fx, Jx, x)

    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
    nr.x .= x
    nr.Fx .= Fx

    f = (norm(Fx)^2)/2
    # this is the gradient
    Fx .= Jx'*Fx
    # As you may notice, this can be expensive... Because the gradient
    # is going to be very simple. May want to create a
    # special type or way to hook into trust regions here. We can exploit
    # that we only need the cauchy and the newton steps, not any shifted
    # systems. There is no need to get the full hessian. because these two
    # steps are don't need these multiplies
    # This is the Hessian
    Jx .= Jx'*Jx
    return f, Fx, Jx
end

function solve(prob::NEqProblem, x, approach::TrustRegion{<:Union{SR1, DBFGS, BFGS, Newton}, <:Any, <:Any}, options::NEqOptions)
    if !(mstyle(prob) === InPlace())
        throw(ErrorException("solve() not defined for OutOfPlace() with Trustregion for NEqProblem"))
    end
    F = prob.R
    # should we wrap a Fx here so we can log F0 info here?
    # and so we can extract it at the end as well?
    # xcache = copy(x).-1
    Fx_outer = copy(x)
    x_outer = copy(x)

    normed_residual = NormedResiduals(x_outer, Fx_outer, F)
    ρ2F0 = sqrt(value(normed_residual, x_outer)*2)
    ρF0 = norm(normed_residual.Fx, Inf)
    td = OptimizationProblem(normed_residual)
    res = solve(td, x, approach, OptimizationOptions(maxiter=options.maxiter))
    # normed_residual to get Fx
    # f0*2 is wrong because it's 2norm
    newinfo = (solution=solution(res), best_residual=Fx_outer, ρF0=ρF0, ρ2F0=ρ2F0, time=res.info.time, iter=res.info.iter)
    return ConvergenceInfo(approach, newinfo, options)
end


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct VoltageModeControl <: InnerControl
        kpv::Float64
        kiv::Float64
        kffv::Float64
        rv::Float64
        lv::Float64
        kpc::Float64
        kic::Float64
        kffi::Float64
        ωad::Float64
        kad::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of an inner loop current control PID using virtual impedance based on D'Arco, Suul and Fosso.
"A Virtual Synchronous Machine implementation for distributed control of power converters in SmartGrids."
Electric Power Systems Research 122 (2015) 180–197.

# Arguments
- `kpv::Float64`: voltage controller proportional gain, validation range: `(0, nothing)`
- `kiv::Float64`: voltage controller integral gain, validation range: `(0, nothing)`
- `kffv::Float64`: Binary variable to enable feed-forward gain of voltage., validation range: `(0, nothing)`
- `rv::Float64`: virtual resistance, validation range: `(0, nothing)`
- `lv::Float64`: virtual inductance, validation range: `(0, nothing)`
- `kpc::Float64`: current controller proportional gain, validation range: `(0, nothing)`
- `kic::Float64`: current controller integral gain, validation range: `(0, nothing)`
- `kffi::Float64`: Binary variable to enable feed-forward gain of current, validation range: `(0, nothing)`
- `ωad::Float64`: active damping filter cutoff frequency (rad/sec), validation range: `(0, nothing)`
- `kad::Float64`: active damping gain, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the VoltageModeControl model are:
	ξd_ic: d-axis integrator state of the PI voltage controller,
	ξq_ic: q-axis integrator state of the PI voltage controller,
	γd_ic: d-axis integrator state of the PI current controller,
	γq_ic: q-axis integrator state of the PI current controller,
	ϕd_ic: d-axis low-pass filter of active damping,
	ϕq_ic: q-axis low-pass filter of active damping
- `n_states::Int`: VoltageModeControl has 6 states
"""
mutable struct VoltageModeControl <: InnerControl
    "voltage controller proportional gain"
    kpv::Float64
    "voltage controller integral gain"
    kiv::Float64
    "Binary variable to enable feed-forward gain of voltage."
    kffv::Float64
    "virtual resistance"
    rv::Float64
    "virtual inductance"
    lv::Float64
    "current controller proportional gain"
    kpc::Float64
    "current controller integral gain"
    kic::Float64
    "Binary variable to enable feed-forward gain of current"
    kffi::Float64
    "active damping filter cutoff frequency (rad/sec)"
    ωad::Float64
    "active damping gain"
    kad::Float64
    ext::Dict{String, Any}
    "The states of the VoltageModeControl model are:
	ξd_ic: d-axis integrator state of the PI voltage controller,
	ξq_ic: q-axis integrator state of the PI voltage controller,
	γd_ic: d-axis integrator state of the PI current controller,
	γq_ic: q-axis integrator state of the PI current controller,
	ϕd_ic: d-axis low-pass filter of active damping,
	ϕq_ic: q-axis low-pass filter of active damping"
    states::Vector{Symbol}
    "VoltageModeControl has 6 states"
    n_states::Int
end

function VoltageModeControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext=Dict{String, Any}(), )
    VoltageModeControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext, [:ξd_ic, :ξq_ic, :γd_ic, :γq_ic, :ϕd_ic, :ϕq_ic], 6, )
end

function VoltageModeControl(; kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext=Dict{String, Any}(), states=[:ξd_ic, :ξq_ic, :γd_ic, :γq_ic, :ϕd_ic, :ϕq_ic], n_states=6, )
    VoltageModeControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function VoltageModeControl(::Nothing)
    VoltageModeControl(;
        kpv=0,
        kiv=0,
        kffv=0,
        rv=0,
        lv=0,
        kpc=0,
        kic=0,
        kffi=0,
        ωad=0,
        kad=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`VoltageModeControl`](@ref) `kpv`."""
get_kpv(value::VoltageModeControl) = value.kpv
"""Get [`VoltageModeControl`](@ref) `kiv`."""
get_kiv(value::VoltageModeControl) = value.kiv
"""Get [`VoltageModeControl`](@ref) `kffv`."""
get_kffv(value::VoltageModeControl) = value.kffv
"""Get [`VoltageModeControl`](@ref) `rv`."""
get_rv(value::VoltageModeControl) = value.rv
"""Get [`VoltageModeControl`](@ref) `lv`."""
get_lv(value::VoltageModeControl) = value.lv
"""Get [`VoltageModeControl`](@ref) `kpc`."""
get_kpc(value::VoltageModeControl) = value.kpc
"""Get [`VoltageModeControl`](@ref) `kic`."""
get_kic(value::VoltageModeControl) = value.kic
"""Get [`VoltageModeControl`](@ref) `kffi`."""
get_kffi(value::VoltageModeControl) = value.kffi
"""Get [`VoltageModeControl`](@ref) `ωad`."""
get_ωad(value::VoltageModeControl) = value.ωad
"""Get [`VoltageModeControl`](@ref) `kad`."""
get_kad(value::VoltageModeControl) = value.kad
"""Get [`VoltageModeControl`](@ref) `ext`."""
get_ext(value::VoltageModeControl) = value.ext
"""Get [`VoltageModeControl`](@ref) `states`."""
get_states(value::VoltageModeControl) = value.states
"""Get [`VoltageModeControl`](@ref) `n_states`."""
get_n_states(value::VoltageModeControl) = value.n_states

"""Set [`VoltageModeControl`](@ref) `kpv`."""
set_kpv!(value::VoltageModeControl, val) = value.kpv = val
"""Set [`VoltageModeControl`](@ref) `kiv`."""
set_kiv!(value::VoltageModeControl, val) = value.kiv = val
"""Set [`VoltageModeControl`](@ref) `kffv`."""
set_kffv!(value::VoltageModeControl, val) = value.kffv = val
"""Set [`VoltageModeControl`](@ref) `rv`."""
set_rv!(value::VoltageModeControl, val) = value.rv = val
"""Set [`VoltageModeControl`](@ref) `lv`."""
set_lv!(value::VoltageModeControl, val) = value.lv = val
"""Set [`VoltageModeControl`](@ref) `kpc`."""
set_kpc!(value::VoltageModeControl, val) = value.kpc = val
"""Set [`VoltageModeControl`](@ref) `kic`."""
set_kic!(value::VoltageModeControl, val) = value.kic = val
"""Set [`VoltageModeControl`](@ref) `kffi`."""
set_kffi!(value::VoltageModeControl, val) = value.kffi = val
"""Set [`VoltageModeControl`](@ref) `ωad`."""
set_ωad!(value::VoltageModeControl, val) = value.ωad = val
"""Set [`VoltageModeControl`](@ref) `kad`."""
set_kad!(value::VoltageModeControl, val) = value.kad = val
"""Set [`VoltageModeControl`](@ref) `ext`."""
set_ext!(value::VoltageModeControl, val) = value.ext = val