text:
# # Langmuir turbulence example
#
# This example implements the Langmuir turbulence simulation reported in section
# 4 of McWilliams, J. C. et al., "Langmuir Turbulence in the ocean," Journal of
# Fluid Mechanics (1997). This example demonstrates:
#
#   * how to run a simulation with surface wave effects via the Craik-Leibovich
#     approximation
#
# In addition to `Oceananigans.jl` we need `PyPlot` for plotting, `Random` for
# generating random initial conditions, and `Printf` for printing progress messages.

using Oceananigans, Oceananigans.Diagnostics, Oceananigans.OutputWriters, Random, Printf
using Oceananigans: g_Earth
using Oceananigans.SurfaceWaves: UniformStokesDrift

# ## Model parameters
#
# Here we use an anisotropic grid with `Nh` horizontal grid points, `Nz` vertical
# grid points, `Δh` horizontal grid spacing, and `Δz` vertical grid spacing.
# We specify fluxes of buouyancy and momentum that match the specifications in
# McWilliams and Sullivan (1997) via
#
#   1. A buoyancy flux `Qᵇ` at the top of the domain.
#
#   2. A velocity flux `Qᵘ` at the top of the domain, which is related
#       to the `x` momentum flux `τˣ` via `τˣ = ρ₀ * Qᵘ`, where `ρ₀` is a reference density.
#
# We also choose the wavenumber and amplitude of the surface wave field, and use an
# initial buoyancy gradient of `N² = 1.936e-5`.

      Nh = 32       # Number of grid points in x, y, z
      Nz = 48       # Number of grid points in x, y, z
      Δh = 4.0      # [m] Grid spacing in x, y, z (meters)
      Δz = 2.0      # [m] Grid spacing in x, y, z (meters)
      Qᵘ = -3.72e-5 # [m² s⁻²] Velocity flux / stress at surface
      Qᵇ = 2.307e-9 # [m³ s⁻²] Buoyancy flux at surface
      N² = 1.936e-5 # [s⁻²] Initial buoyancy gradient
       f = 1e-4     # [s⁻¹] Coriolis parameter
end_time = 2hour    # [s] End time for the simulation

# Surface wave stokes drift profile:
#
# $ uˢ(z) = Uˢ exp(2kˢʷ z) $
#
# and
#
# $ ∂_z uˢ(z) = 2 kˢʷ Uˢ exp(2kˢʷ z) $
#
# Notes:
#
#  * We use `const`s so Stokes drift functions compile on GPU.
#
#  * The `Oceananigans.jl` implementation of Craik-Leibovich uses the Lagrangian-mean
#    formulation, which means that the vertical derivative or 'Stokes shear' is specified.
#    If the surface wave field Stokes drift changes in time, the time-derivative of the
#    Stokes drift profile must be provided as well.

const kˢʷ = 0.105    # [m⁻¹] Surface wave wavenumber
const aˢʷ = 0.8      # [m] Surface wave amplitude
const Uˢ = aˢʷ^2 * kˢʷ * sqrt(g_Earth * kˢʷ)

## Oceananigans
∂z_uˢ(z, t) = 2kˢʷ * Uˢ * exp(2kˢʷ * z)

# We use the Stokes drift profile for the initial condition
uˢ(z) = Uˢ * exp(2kˢʷ * z)

# ## Boundary conditions
#
# Here we define `Flux` boundary conditions at the surface for `u`, `T`, and `S`,
# and a `Gradient` boundary condition on `T` that maintains a constant stratification
# at the bottom. Our flux boundary condition for salinity uses a function that calculates
# the salinity flux in terms of the evaporation rate.

grid = RegularCartesianGrid(size=(Nh, Nh, Nz), extent=(Δh*Nh, Δh*Nh, Δz*Nz))

u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))

b_bcs = TracerBoundaryConditions(grid,    top = BoundaryCondition(Flux, Qᵇ),
                                       bottom = BoundaryCondition(Gradient, N²))

# ## Model instantiation
#
# We instantiate a horizontally-periodic `Model` on the CPU with on a `RegularCartesianGrid`,
# using a `FPlane` model for rotation (constant rotation rate), the Anisotropic Minimum
# Dissipation closure to model the effects of unresolved turbulence, and the previously defined
# boundary conditions for `u` and `b`.

model = Model(
         architecture = CPU(),
                 grid = grid,
             buoyancy = BuoyancyTracer(), tracers = :b,
              closure = AnisotropicMinimumDissipation(),
        surface_waves = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
  boundary_conditions = SolutionBoundaryConditions(grid, u=u_bcs, b=b_bcs),
)

# ## Initial conditions
#
# Out initial condition for buoyancy consists of a linear stratification superposed with
# random noise damped at the walls.

## Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

## Temperature initial condition: a stable density tradient with random noise superposed.
b₀(x, y, z) = N² * z + 1e-1 * Ξ(z) * N² * model.grid.Lz

# Velocity initial condition: note that McWilliams and Sullivan (1997) use a model
# that is formulated in terms of the Eulerian-mean velocity field. Due to this,
# their 'resting' initial condition actually consists of a sheared and
# unbalanced Lagrangian-mean velocity field uˢ(z). To this we add noise.
u₀(x, y, z) = uˢ(z) + sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

# Finally, we introduce noise into the vertical velocity:
w₀(x, y, z) = sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

set!(model, u=u₀, w=w₀, b=b₀)

# ## Set up output
#
# We set up an output writer that saves all velocity fields, tracer fields, and the subgrid
# turbulent diffusivity associated with `model.closure`. The `prefix` keyword argument
# to `JLD2OutputWriter` indicates that output will be saved in
# `ocean_wind_mixing_and_convection.jld2`.

## Create a NamedTuple containing all the fields to be outputted.
fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,))

## Instantiate a JLD2OutputWriter to write fields.
field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4,
                                prefix="langmuir_turbulence", force=true)

## Add the output writer to the models `output_writers`.
model.output_writers[:fields] = field_writer;

# ## Running the simulation
#
# To run the simulation, we instantiate a `TimeStepWizard` to ensure stable time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 0.2.

wizard = TimeStepWizard(cfl=0.2, Δt=5.0, max_change=1.1, max_Δt=10.0)

# Diagnostic that returns the maximum absolute value of `u, v, w` by calling
# `umax(), vmax(), `wmax()`:

umax = FieldMaximum(abs, model.velocities.u)
vmax = FieldMaximum(abs, model.velocities.v)
wmax = FieldMaximum(abs, model.velocities.w)

# We also create a figure and define a plotting function for live plotting of results.
makeplot = false

using PyPlot

fig, axs = subplots(ncols=3, figsize=(12, 5))

"""
    makeplot!(axs, model)

Make a triptych of x-z slices of vertical velocity, temperature, and salinity
associated with `model` in `axs`.
"""
function makeplot!(axs, model)
    jhalf = floor(Int, model.grid.Ny/2)

    ## Coordinate arrays for plotting
    xCˣᶻ = repeat(model.grid.xC, 1, model.grid.Nz)
    xFˣᶻ = repeat(model.grid.xF[1:end-1], 1, model.grid.Nz)
    zCˣᶻ = repeat(reshape(model.grid.zC, 1, model.grid.Nz), model.grid.Nx, 1)
    zFˣᶻ = repeat(reshape(model.grid.zF[1:end-1], 1, model.grid.Nz), model.grid.Nx, 1)

    xCˣʸ = repeat(model.grid.xC, 1, model.grid.Ny)
    yCˣʸ = repeat(reshape(model.grid.yC, 1, model.grid.Ny), model.grid.Nx, 1)

    sca(axs[1]); cla()
    title("Horizontal velocity")
    pcolormesh(xFˣᶻ, zCˣᶻ, Array(interior(model.velocities.u))[:, jhalf, :])
    xlabel("\$ x \$ (m)"); ylabel("\$ z \$ (m)")

    sca(axs[2]); cla()
    title("Vertical velocity")
    pcolormesh(xCˣᶻ, zFˣᶻ, Array(interior(model.velocities.w))[:, jhalf, :])
    xlabel("\$ x \$ (m)")

    sca(axs[3]); cla()
    title("Vertical velocity at \$ z = $(model.grid.zF[Nz-2]) \$ meters")
    pcolormesh(xCˣʸ, yCˣʸ, Array(interior(model.velocities.w))[:, :, Nz-2])
    xlabel("\$ x \$ (m)"); ylabel("\$ y \$ (m)")

    axs[3].yaxis.set_label_position("right")
    axs[3].tick_params(right=true, labelright=true, left=false, labelleft=false)

    [ax.set_aspect(1) for ax in axs]

    tight_layout()

    pause(0.01)

    return nothing
end

# Finally, we run the the model in a `while` loop.

while model.clock.time < end_time

    ## Update the time step associated with `wizard`.
    update_Δt!(wizard, model)

    ## Time step the model forward
    walltime = @elapsed time_step!(model, 100, wizard.Δt)

    ## Print a progress message
    @printf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
            model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
            umax(), vmax(), wmax(), prettytime(walltime))

    makeplot && makeplot!(axs, model)
end

# Show the reults in a plot

makeplot!(axs, model)
gcf()


text:
# For testing this config re-dispatches Xrule_via_ad to Xrule without config argument
struct ADviaRuleConfig <: RuleConfig{Union{HasReverseMode, HasForwardsMode}} end

function ChainRulesCore.frule_via_ad(config::ADviaRuleConfig, ȧrgs, f, args...; kws...)
    ret = frule(config, ȧrgs, f, args...; kws...)
    # we don't support actually doing AD: the rule has to exist. lets give helpfulish error
    ret === nothing && throw(MethodError(frule, (ȧrgs, f, args...)))
    return ret
end

function ChainRulesCore.rrule_via_ad(config::ADviaRuleConfig, f, args...; kws...)
    ret = rrule(config, f, args...; kws...)
    # we don't support actually doing AD: the rule has to exist. lets give helpfulish error
    ret === nothing && throw(MethodError(rrule, (f, args...)))
    return ret
end

ext: jl


text:
module SavitzkyGolay

export savitzky_golay, SGolay, SGolayResults

using LinearAlgebra

struct SGolay{T1 <: Signed, T2 <: Real}
    w::T1        # Window size
    order::T1    # Polynomial order
    deriv::T1    # Derivative order
    rate::T2      # Rate
    function SGolay(w::T1, order::T1, deriv::T1, rate::T2) where {T1 <: Signed, T2 <: Real}
        isodd(w) || throw(ArgumentError("w must be an odd number."))
        w ≥ 1 || throw(ArgumentError("w must greater than or equal to 1."))
        w ≥ order + 2 || throw(ArgumentError("w too small for the polynomial order chosen (w ≥ order + 2)."))
        return new{T1, T2}(w, order, deriv, rate)
    end
end

SGolay(w, order) = SGolay(w, order, 0, 1.0)
SGolay(w, order, deriv) = SGolay(w, order, deriv, 1.0)

struct SGolayResults{T <: Float64}
    y::Vector{T}
    params::SGolay
    coeff::Vector{T}
    Vdm::Matrix{T}
end

function (p::SGolay)(y::AbstractVector)
    return _savitzky_golay(y, p)
end

function savitzky_golay(
    y::AbstractVector, window_size::T0, order::T0;
    deriv::T0=0, rate::T1=1.0,
    ) where {T0 <: Signed, T1 <: Real}
    y_, p = _check_input_sg(y, window_size, order, deriv, rate)
    return _savitzky_golay(y_, p)
end

function _savitzky_golay(y::AbstractVector, p::SGolay)
    order_range = 0 : p.order
    hw = Int64((p.w - 1) / 2) # half-window size
    V = zeros(2*hw + 1, length(order_range))
    _vandermonde!(V, hw, order_range)
    c = _coefficients(V, order_range, p)
    y_ = _padding_signal(y, hw)
    y_conv = _convolve_1d(y_, c)
    return SGolayResults(y_conv, p, c, V)
end

function _check_input_sg(y::Vector, w, order, deriv, rate)
    isodd(w) || throw(ArgumentError("w must be an even number."))
    w ≥ 1 || throw(ArgumentError("w must greater than or equal to 1."))
    w ≥ order + 2 || throw(ArgumentError("w too small for the polynomial order chosen (w ≥ order + 2)."))
    length(y) > 1 || throw(ArgumentError("vector x must have more than one element."))
    return Float64.(y), SGolay(w, order, deriv, rate)
end

function _convolve_1d(u::AbstractVector, v::Vector)
    m = length(u)
    n = length(v)
    m > n || throw(ArgumentError("length of signal u must be greater than length of kernel v."))
    w = zeros(m + n - 1)
    @inbounds for j in 1:m, k in 1:n
        w[j+k-1] += u[j]*v[k]
    end
    return w[n:end-n+1]
end

function _vandermonde!(V::Matrix{T1}, hw::T2, order_range::UnitRange{T2}) where {T1 <: Float64, T2 <: Int64}
    @inbounds for i in -hw:hw, j in order_range
        V[i+hw+1, j+1] = i^j
    end
end

function _coefficients(V::Matrix{T1}, order_range::UnitRange{T2}, p::SGolay) where {T1 <: Float64, T2 <: Int64}
    Vqr = qr(V')
    c = Vqr.R \ (Vqr.Q' * _onehot(p.deriv + 1, length(order_range)))
    c .*= (p.rate)^(p.deriv) * factorial(p.deriv)
    return reverse(c)
end

function _onehot(i::T, m::T) where T <: Int64
    m > i || throw(ArgumentError("length of vector must be greater than the position"))
    oh = zeros(m)
    oh[i] = 1.0
    return oh
end

function _padding_signal(y::AbstractVector, hw::Int64)
    initvals = y[1] .- abs.(reverse(y[2:hw+1]) .- y[1])
    endvals = y[end] .+ abs.(reverse(y[end-hw:end-1] .- y[end]))
    return vcat(initvals, y, endvals)
end

end  # module SavitzkyGolay


text:
```@meta
DocTestSetup = quote
    using DataFrames, MixedModelsSim, Tables
end
```

"""
    factorproduct(facs...)

Return a `Vector{NamedTuple}` obtained by crossing `facs...`.

The arguments should be coerceable to a `Tables.RowTable` with `rowtable`.

The value is a `Tables.RowTable` and hence can be converted to a `DataFrame`.

# Example
```julia-repl
julia> DataFrame(factorproduct((item=nlevels(3,'I'),), (subj=nlevels(5), age=["Y","Y","Y","O","O"])))
15×3 DataFrame
│ Row │ item   │ subj   │ age    │
│     │ String │ String │ String │
├─────┼────────┼────────┼────────┤
│ 1   │ I1     │ S1     │ Y      │
│ 2   │ I2     │ S1     │ Y      │
│ 3   │ I3     │ S1     │ Y      │
│ 4   │ I1     │ S2     │ Y      │
│ 5   │ I2     │ S2     │ Y      │
│ 6   │ I3     │ S2     │ Y      │
│ 7   │ I1     │ S3     │ Y      │
│ 8   │ I2     │ S3     │ Y      │
│ 9   │ I3     │ S3     │ Y      │
│ 10  │ I1     │ S4     │ O      │
│ 11  │ I2     │ S4     │ O      │
│ 12  │ I3     │ S4     │ O      │
│ 13  │ I1     │ S5     │ O      │
│ 14  │ I2     │ S5     │ O      │
│ 15  │ I3     │ S5     │ O      │
```
"""
factorproduct(facs...) = factorproduct(Val(:rowtable), facs...)

function factorproduct(::Val{:rowtable}, facs...)
    vec([merge(s...) for s in Iterators.product(rowtable.(facs)...)])
end

# currently not more efficient, but reserving this use for later efficiency improvements
factorproduct(::Val{:columntable}, facs...) = columntable(factorproduct(facs...))

"""
    cyclicshift(v::AbstractVector, nrow)

Return an `eltype(v)` matrix of size `nrow` by `length(v)` with each column consisting
of `v` followed by a cyclic shift of `v` followed by ...

```@example
cyclicshift('a':'d', 8)
```
"""
function cyclicshift(v::AbstractVector, nrow)
    # The return value i used to counterbalance levels of conditions in [`withinitem`](@ref).
    vlen = length(v)
    [v[1 + (i + j) % vlen] for i in 0:(nrow - 1), j in 0:(vlen - 1)]
end

"""
    nlevels(nlev, tag='S')

Return a `Vector{String}` of `tag` followed by `1:nlev` left-padded with zeros

# Examples
```julia-repl
julia> show(nlevels(10))
["S01", "S02", "S03", "S04", "S05", "S06", "S07", "S08", "S09", "S10"]
```
"""
function nlevels(nlev, tag='S')
    string.(tag, lpad.(1:nlev, ndigits(nlev), '0'))
end

"""
    nlevstbl(nm::Symbol, n, vars::Pair{Symbol, Vector{String}}...)

Return a `Tables.columntable` with a `nm` column as a `PooledArray` with `n` levels.

If any `vars` pairs are given they are expanded to columns representing characteristics
of the `nm` column.  In experimental design terminology, if say `nm` is `:item` then these
represent between-item experimental factors.

The `nm` column is generated as `nlevels(n, uppercase(first(string(nm))))`

# Examples
```julia-repl
julia> nlevstbl(:item, 10)
(item = ["I01", "I02", "I03", "I04", "I05", "I06", "I07", "I08", "I09", "I10"],)

julia> nlevstbl(:item, 9, :level => ["low", "medium", "high"])
(item = ["I1", "I2", "I3", "I4", "I5", "I6", "I7", "I8", "I9"], level = ["low", "medium", "high", "low", "medium", "high", "low", "medium", "high"])
```
"""
function nlevstbl(nm::Symbol, n::Integer, vars::Pair{Symbol, Vector{String}}...)
    nms = [nm]
    # need to specify PooledArray[] so that the array doesn't specialize on
    # PooledArrays's type parameters
    vals = PooledArray[PooledArray(nlevels(n, uppercase(first(string(nm)))); signed=true, compress=true)]
    inner = 1
    for var in vars
        levs = last(var)
        nlev = length(levs)
        rept = inner * nlev
        q, r = divrem(n, rept)
        iszero(r) || throw(ArgumentError("n = $n is not a multiple of repetition block size $rept"))
        push!(vals, PooledArray(repeat(levs, inner=inner, outer=q); signed=true, compress=true))
        push!(nms, first(var))
        inner = rept
    end
    NamedTuple{(nms...,)}((vals...,))
end

"""
    pooled!(df, cols::Type=Union{AbstractString,Missing})

Like `DataFrames.categorical!` but converting columns to `PooledArray`s

!!! warning
    This method is not type-specific in the first argument, in order to eliminate
    a dependency on `DataFrames.jl`. It nonetheless expects a `DataFrame` as
    its first argument
"""
function pooled!(df, cols::Type=Union{AbstractString,Missing})
    try
        for n in propertynames(df)
            v = df[!, n]
            if eltype(v) <: cols
                df[!, n] = PooledArray(v)
            end
        end
    catch e
        if e isa MethodError
            throw(ArgumentError("Something went wrong. Make sure you're using a DataFrame."))
        else
            rethrow(e)
        end
    end

    return df
end

"""
    flatlowertri(::LowerTriangular)

Returns the lower triangular flattened into 1D array in column-major order.
"""
function flatlowertri(l::LowerTriangular)
    rr, cc = size(l)
    return [l[i, j] for j = 1:cc for i = j:rr]
end

"""
    _update!(m::MixedModel, θ)

Update the mixed model to use θ as its new parameter vector.

!!! note
    This is a convenience function for installing a particular parameter vector
    and the resulting model fit. It does not actually perform any type of
    optimization.

!!! note
    For GLMMs, this only sets θ and not β, even for `fast=false` fits.
"""
_update!(m::LinearMixedModel, θ) = updateL!(MixedModels.setθ!(m, θ))
_update!(m::GeneralizedLinearMixedModel, θ) = pirls!(MixedModels.setθ!(m, θ), false)
# arguably type piracy, but we're all the same developers....

"""
    update!(m::MixedModel, re...)

Update the mixed model to use the random-effects covariance matrices.

The `re` can be created using [`create_re`](@ref).

They should be specified in the order specified in `VarCorr(m)`.

!!! warning
    We recommend against calling this method directly. Instead, use the method
    with keyword arguments to specify the  different `re` by name.

!!! note
    This is a convenience function for installing a particular parameter vector
    and the resulting model fit. It does not actually perform any type of
    optimization.

Details
========
The `re` used as the λ fields of the model's `ReTerm`s and should be specified
as the lower Cholesky factor of covariance matrices.
"""
function update!(m::MixedModel, re...)
    Base.depwarn("Specifying the random effects by position instead of name is deprecated",
                 :update!vararg)

    θ = vcat((flatlowertri(rr) for rr in re)...)
    return _update!(m, θ)
end

"""
    update!(m::MixedModel; namedre...)
    update!(m::MixedModel; θ)

Update the mixed model to use the random-effects covariance matrices.

The `namedre` can be created using [`create_re`](@ref). The `namedre` are specified
by the name of the blocking variable, e.g. `subj=create_re(...)`.

!!! warning
    Setting θ directly as a keyword-argument is deprecated.

!!! note
    This is a convenience function for installing a particular parameter vector
    and the resulting model fit. It does not actually perform any type of
    optimization.

Details
========
The `re` is used as the λ fields of the model's `ReTerm`s and should be specified
as the lower Cholesky factor of covariance matrices.
"""
function update!(m::MixedModel; θ=nothing, namedre...)
    if θ === nothing
        length(namedre) == 0 && throw(ArgumentError("no random effects specified"))
        θ = createθ(m; namedre...)
    elseif θ !== nothing && length(namedre) != 0
        throw(ArgumentError("θ must not be specified when named random effects are"))
    end
    return _update!(m, θ)
end

"""
    createθ(m::MixedModel; named_re)
    create_theta(m::MixedModel; named_re)

Create the parameter vector θ corresponding to the random effects.

The `named_re` can be created using [`create_re`](@ref). The `named_re` are specified
by the name of the blocking variable, e.g. `subj=create_re(...)`.

The model must be specified because the parameters are sorted internally for computational
efficiency.
"""
function createθ(m::MixedModel; named_re...)
    newre = Dict(named_re...)
    fns = fnames(m)

    if Set(fns) != keys(newre)
        throw(ArgumentError("Specified blocking variables do not match model blocking variables."))
    end

    return mapfoldl(vcat, fns) do fname
        flatlowertri(newre[fname])
    end
end

const create_theta = createθ

"""
    create_re(sigmas...; corrmat=Matrix{Float64}(I, length(sigmas), length(sigmas))

Create the covariance factor for a random effect from the standard deviations and correlation matrix.

The `sigmas` should be specified in the same order as the random slopes in the
output of `VarCorr(m)`.

The correlation matrix defaults to the identity matrix, i.e. no correlation between
random effects.

!!! note
    The return value is the lower Cholesky factor of the covariance matrix, which is what
    [`update!`](@ref) requires.
"""
function create_re(sigmas...; corrmat=nothing)
    ss = Diagonal([sigmas...])

    isnothing(corrmat) ? LowerTriangular(ss) : ss * cholesky(Symmetric(corrmat, :L)).L
end

# """
#     withinitem(nitem, df)

# Return a `DataFrame` of `item` with `nitem` levels and balanced columns of conditions from `df`
# """
# function withinitem(nitem, df; tag = 'I')
#     nrow = size(df, 1)
#     q, r = divrem(nitem, nrow)
#     iszero(r) || throw(ArgumentError("nitem = $nitem is not a multiple of nrow = $nrow"))
#     value = df[vec(cyclicshift(1:nrow, nitem)), :]
#     value.item = repeat(PooledArray(nlevels(nitem, tag=tag)), outer=nrow)
#     value
# end


text:
module Day16

export day16

struct Ticket
    fields::Vector{Int}
end

struct Rule
    name::String
    range1::UnitRange
    range2::UnitRange
end

using CombinedParsers
using CombinedParsers.Regexp: words, whitespace_maybe

function parserule(str)
    line = Sequence(:name => CombinedParsers.Regexp.words, ": ", :first => Numeric(Int), '-', :second => Numeric(Int), " or ", :third => Numeric(Int), '-', :fourth =>Numeric(Int))
    t = str |> line
    Rule(t[:name], UnitRange(t[:first],t[:second]), UnitRange(t[:third],t[:fourth]))
end

function parseticket(str)
    fields = join(Repeat(Numeric(Int)), ',')
    str |> fields
end

function parseinput(filename)
    rules = Dict{String,Tuple{UnitRange,UnitRange}}()
    myticket = Ticket(Int[])
    tickets = Vector{Ticket}()
    mode = :rules
    for line in readlines(filename)
        if line == "your ticket:"
            mode = :myticket
            continue
        elseif line == "nearby tickets:"
            mode = :tickets
            continue
        elseif line |> strip == ""
            continue
        end
        if mode == :rules
            rule = parserule(line)
            rules[rule.name] = (rule.range1, rule.range2)
        elseif mode == :myticket
            myticket = Ticket(parseticket(line))
        elseif mode == :tickets
            push!(tickets, Ticket(parseticket(line)))
        end
    end

    (rules, myticket, tickets)
end

function iscompletelyinvalid(rules, n)
    for t in values(rules)
        if in(n, t[1]) || in(n, t[2])
            return 0
        end
    end
    n
end

function day16(::Val{:a}, filename)
    (rules, myticket, tickets) = parseinput(filename)

    suma = 0
    for ticket in tickets
        for n in ticket.fields
            suma += iscompletelyinvalid(rules, n)
        end
    end
    suma
end

isvalid(rules, ticket) = Base.all(Base.any(in(n, t[1]) || in(n, t[2]) for t in values(rules)) for n in ticket.fields)

function correct!(labels, rules, ticket)
    for (i,n) in enumerate(ticket.fields)
        for (name, t) in rules
            if !in(n, t[1]) && !in(n, t[2])
                delete!(labels[i], name)
            end
        end
    end
end

function day16(::Val{:b}, filename)
    (rules, myticket, tickets) = parseinput(filename)
    # Find the valid tickets
    valid_tickets = Ticket[]
    for ticket in tickets
        if isvalid(rules, ticket)
            push!(valid_tickets, ticket)
        end
    end

    # Initialize labels: every field may have any label
    labels = Dict{Int, Set{String}}()
    for field in 1:length(myticket.fields)
        labels[field] = Set(keys(rules))
    end

    # Remove labels from fields that cannot have the
    # label due to range constraints
    for ticket in valid_tickets
        correct!(labels, rules, ticket)
    end

    # Start with fields with only one possibility.
    # Remove this possibility from the other field's
    # label set. Continue until all fields have a
    # unique label
    label = Dict{Int, String}()
    todo = Set(1:length(myticket.fields))
    while !isempty(todo)
        for (n,s) in labels
            if length(s) == 1
                l = pop!(s)
                label[n] = l
                for (m,t) in labels
                    delete!(t, l)
                end
                delete!(todo, n)
            end
        end
    end

    pr = 1
    for (i,n) in enumerate(myticket.fields)
        if length(label[i]) < 9
            continue
        elseif label[i][1:9] == "departure"
            pr *= n
        end
    end
    pr
end

end


text:
const tbbox = (min_x = 440720.0, min_y = 3.74532e6, max_x = 446720.0, max_y = 3.75132e6)

@testset "GeoUtils" begin
    @testset "bbox" begin
        ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
        @test GeoArrays.bbox(ga) == tbbox
    end

    @testset "bbox!" begin
        ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
        @test GeoArrays.bbox(ga) == tbbox
        GeoArrays.bbox!(ga, GeoArrays.bbox(ga))
        @test GeoArrays.bbox(ga) == tbbox
    end

    @testset "flip upside down" begin
        ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
        ga_identical = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
        @test ga.A == ga_identical.A
        @test ga.f == ga_identical.f
        @test ga.crs == ga_identical.crs
        GeoArrays.flipud!(ga)
        @test ga.f != ga_identical.f
        @test GeoArrays.bbox(ga) == tbbox
        GeoArrays.flipud!(ga)
        @test ga.f == ga_identical.f
        @test ga.A == ga_identical.A
    end

    @testset "geotransform should remain the same" begin
        gt = [1000.,1.,0.2,2000.,0.,1.]
        f = GeoArrays.geotransform_to_affine(gt)
        gtnew = GeoArrays.affine_to_geotransform(f)
        @test gt == gtnew
    end

    @testset "bbox overlapping" begin
        @test GeoArrays.bbox_overlap((min_x = 10, min_y = 10, max_x = 1000, max_y = 1000), (min_x = 20, min_y = 20, max_x = 900, max_y = 1100)) == true
        @test GeoArrays.bbox_overlap((min_x = 10, min_y = 10, max_x = 100, max_y = 100), (min_x = 200, min_y = 200, max_x = 400, max_y = 400)) == false
        @test GeoArrays.bbox_overlap((min_x = 440720.0, min_y = 3.74532e6, max_x = 446720.0, max_y = 3.75132e6), (min_x = 442720.0, min_y = 3.74532e6, max_x = 456720.0, max_y = 3.75032e6)) == true
        @test GeoArrays.bbox_overlap((min_x = 440720.0, min_y = 3.74532e6, max_x = 446720.0, max_y = 3.75132e6), (min_x = 342720.0, min_y = 3.94532e6, max_x = 356720.0, max_y = 3.95032e6)) == false
    end

    @testset "crop by bbox or ga" begin
        ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
        cbox = (min_x = 441260.0, min_y = 3.75012e6, max_x = 442520.0, max_y = 3.75078e6)
        ga_s = ga[10:30,20:50,begin:end]

        ga_c = GeoArrays.crop(ga, cbox)
        ga_s_c = GeoArrays.crop(ga, ga_s)

        @test size(ga_c) == (21, 11, 1)
        @test GeoArrays.bbox(ga_c) == cbox

        @test size(ga_s_c) == size(ga_s)
        @test GeoArrays.bbox(ga_s_c) == GeoArrays.bbox(ga_s)
    end

end


text:
#=
estimate a covariance matrix that satisfies some given prior information and
minimizes a given convex spectral function

P ∈ 𝕊ᵈ is the covariance variable
minimize    f(P)                    (note: enforces P ⪰ 0)
subject to  tr(P) = 1               (normalize)
            B vec(P) = b            (prior info as equalities)
            C vec(P) ≤ c            (prior info as inequalities)
where f is a convex spectral function
=#

struct CovarianceEstJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    d::Int
    ext::MatSpecExt # formulation specifier
end

function build(inst::CovarianceEstJuMP{T}) where {T <: Float64}
    d = inst.d
    @assert d >= 1
    @assert is_domain_pos(inst.ext)

    P0 = randn(T, d, d)
    P0 = Hermitian(P0 * P0' + 0.5 * I, :U)
    P0 *= d / tr(P0)
    vec_dim = Cones.svec_length(d)
    P0_vec = zeros(T, vec_dim)
    Cones.smat_to_svec!(P0_vec, P0, sqrt(T(2)))

    model = JuMP.Model()
    JuMP.@variable(model, P[1:d, 1:d], Symmetric)
    JuMP.@constraint(model, tr(P) == d)
    P_vec = zeros(JuMP.AffExpr, vec_dim)
    Cones.smat_to_svec!(P_vec, one(T) * P, sqrt(T(2)))

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(inst.ext, d, vcat(1.0 * epi, P_vec), model)

    # linear prior constraints
    B = randn(T, div(d, 3), vec_dim)
    b = B * P0_vec
    JuMP.@constraint(model, B * P_vec .== b)
    C = randn(T, round(Int, log(d)), vec_dim)
    c = C * P0_vec
    JuMP.@constraint(model, C * P_vec .<= c)

    # save for use in tests
    model.ext[:P_var] = P

    return model
end

function test_extra(inst::CovarianceEstJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective
    tol = eps(T)^0.2
    P_opt = JuMP.value.(model.ext[:P_var])
    λ = eigvals(Symmetric(P_opt, :U))
    @test minimum(λ) >= -tol
    obj_result = get_val(pos_only(λ), inst.ext)
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

#=
  File i/o for nii/mgh files based on MATLAB code by Doug Greve, Bruce Fischl:
    MRIfspec.m
    MRIread.m
    MRIwrite.m
    load_mgh.m
    load_nifti.m
    load_nifti_hdr.m
    save_nifti.m
    fsgettmppath.m
    vox2ras_0to1.m
    vox2ras_tkreg.m
    vox2rasToQform.m
=#

using LinearAlgebra, Printf, DelimitedFiles

export MRI, NIfTIheader, get_tmp_path, mri_filename, mri_read, mri_write,
       mri_read_bfiles, mri_read_bfiles!

"Container for header of a volume stored in NIfTI format"
struct NIfTIheader
  # NIfTI standard header fields
  sizeof_hdr::Int32
  data_type::NTuple{10, UInt8}
  db_name::NTuple{18, UInt8}
  extents::Int32
  session_error::Int16
  regular::UInt8
  dim_info::UInt8
  dim::NTuple{8,Int16}
  intent_p1::Float32
  intent_p2::Float32
  intent_p3::Float32
  intent_code::Int16
  datatype::Int16
  bitpix::Int16
  slice_start::Int16
  pixdim::NTuple{8,Float32}
  vox_offset::Float32
  scl_slope::Float32
  scl_inter::Float32
  slice_end::Int16
  slice_code::Int8
  xyzt_units::Int8
  cal_max::Float32
  cal_min::Float32
  slice_duration::Float32
  toffset::Float32
  glmax::Int32
  glmin::Int32
  descrip::NTuple{80,UInt8}
  aux_file::NTuple{24,UInt8}
  qform_code::Int16
  sform_code::Int16
  quatern_b::Float32
  quatern_c::Float32
  quatern_d::Float32
  quatern_x::Float32
  quatern_y::Float32
  quatern_z::Float32
  srow_x::NTuple{4,Float32}
  srow_y::NTuple{4,Float32}
  srow_z::NTuple{4,Float32}
  intent_name::NTuple{16,UInt8}
  magic::NTuple{4,UInt8}

  # Additional fields
  do_bswap::Bool
  sform::Matrix{Float32}
  qform::Matrix{Float32}
  vox2ras::Matrix{Float32}
end

"Container for header and image data of an MRI volume or volume series"
mutable struct MRI{A<:AbstractArray}
  vol::A
  ispermuted::Bool
  niftihdr::NIfTIheader

  fspec::String
  pwd::String

  flip_angle::Float32
  tr::Float32
  te::Float32
  ti::Float32

  vox2ras0::Matrix{Float32}
  volsize::Vector{Int32}
  height::Int32
  width::Int32
  depth::Int32
  nframes::Int32

  vox2ras::Matrix{Float32}
  nvoxels::Int32
  xsize::Float32
  ysize::Float32
  zsize::Float32

  x_r::Float32
  x_a::Float32
  x_s::Float32

  y_r::Float32
  y_a::Float32
  y_s::Float32

  z_r::Float32
  z_a::Float32
  z_s::Float32

  c_r::Float32
  c_a::Float32
  c_s::Float32

  vox2ras1::Matrix{Float32}
  Mdc::Matrix{Float32}
  volres::Vector{Float32}
  tkrvox2ras::Matrix{Float32}

  bval::Vector{Float32}
  bvec::Matrix{Float32}
end

"""
    MRI(vol::Array{T,N}) where T<:Number

Return an empty `MRI` structure
"""
MRI(vol::Array{T}) where T<:Number = MRI(
  vol,
  false,
  NIfTIheader(Int32(0),
              Tuple(zeros(UInt8, 10)),
              Tuple(zeros(UInt8, 18)),
              Int32(0),
              Int16(0),
              UInt8(0),
              UInt8(0),
              Tuple(zeros(Int16, 8)),
              Float32(0),
              Float32(0),
              Float32(0),
              Int16(0),
              Int16(0),
              Int16(0),
              Int16(0),
              Tuple(zeros(Float32, 8)),
              Float32(0),
              Float32(0),
              Float32(0),
              Int16(0),
              Int8(0),
              Int8(0),
              Float32(0),
              Float32(0),
              Float32(0),
              Float32(0),
              Int32(0),
              Int32(0),
              Tuple(zeros(UInt8, 80)),
              Tuple(zeros(UInt8, 24)),
              Int16(0),
              Int16(0),
              Float32(0),
              Float32(0),
              Float32(0),
              Float32(0),
              Float32(0),
              Float32(0),
              Tuple(zeros(Float32, 4)),
              Tuple(zeros(Float32, 4)),
              Tuple(zeros(Float32, 4)),
              Tuple(zeros(UInt8, 16)),
              Tuple(zeros(UInt8, 4)),

              UInt8(0),
              Matrix{Float32}(undef, 0, 0),
              Matrix{Float32}(undef, 0, 0),
              Matrix{Float32}(undef, 0, 0)
             ),

  "",
  "",

  Float32(0),
  Float32(0),
  Float32(0),
  Float32(0),

  Matrix{Float32}(undef, 0, 0),
  Vector{Int32}(undef, 0),
  Int32(0),
  Int32(0),
  Int32(0),
  Int32(0),

  Matrix{Float32}(undef, 0, 0),
  Int32(0),
  Float32(0),
  Float32(0),
  Float32(0),

  Float32(0),
  Float32(0),
  Float32(0),

  Float32(0),
  Float32(0),
  Float32(0),

  Float32(0),
  Float32(0),
  Float32(0),

  Float32(0),
  Float32(0),
  Float32(0),

  Matrix{Float32}(undef, 0, 0),
  Matrix{Float32}(undef, 0, 0),
  Vector{Float32}(undef, 0),
  Matrix{Float32}(undef, 0, 0),

  Vector{Float32}(undef, 0),
  Matrix{Float32}(undef, 0, 0)
)

"""
    MRI(ref::MRI{A}, nframes::Integer=ref.nframes, datatype::DataType=eltype(A)) where A<:AbstractArray

Return an `MRI` structure whose header fields are populated based on a
reference `MRI` structure `ref`, and whose image array are populated with
zeros.

Optionally, the new `MRI` structure can be created with a different number of
frames (`nframes`) than the reference MRI structure.
"""
function MRI(ref::MRI{A}, nframes::Integer=ref.nframes, datatype::DataType=eltype(A)) where A<:AbstractArray

  if nframes == 1
    mri = MRI(zeros(datatype, ref.volsize...))
  else
    mri = MRI(zeros(datatype, ref.volsize..., nframes))
  end

  for var in fieldnames(MRI)
    any(var .== (:vol, :fspec, :bval, :bvec)) && continue
    setfield!(mri, var, getfield(ref, var))
  end

  mri.nframes = nframes

  return mri
end

"""
    get_tmp_path(tmpdir::String="")

Return path to a directory where temporary files can be stored.

Search for candidate directories in the following order:
 1. `\$``TMPDIR`:  Check if environment variable is defined and directory exists
 2. `\$``TEMPDIR`: Check if environment variable is defined and directory exists
 3. `/scratch`:    Check if directory exists
 4. `/tmp`:        Check if directory exists
 5. `tmpdir`:      Check if `tmpdir` argument was passed and directory exists

If none of the above exist, use current directory (`./`) and print warning.
"""
function get_tmp_path(tmpdir::String="")

  if haskey(ENV, "TMPDIR")
    tmppath = ENV["TMPDIR"]
    if isdir(tmppath)
      return tmppath
    end
  end

  if haskey(ENV, "TEMPDIR")
    tmppath = ENV["TEMPDIR"]
    if isdir(tmppath)
      return tmppath
    end
  end

  tmppath = "/scratch"
  if isdir(tmppath)
    return tmppath
  end

  tmppath = "/tmp"
  if isdir(tmppath)
    return tmppath
  end

  tmppath = tmpdir
  if isdir(tmppath)
    return tmppath
  end

  tmppath = "./"
  println("WARNING: get_tmp_path could not find a temporary folder, " *
	  "using current folder")
  return tmppath
end

"""
    vox2ras_0to1(M0::Matrix)

Convert a 0-based vox2ras matrix `M0` to a 1-based vox2ras matrix such that:

Pxyz = M_0 * [c r s 1]' = M_1 * [c+1 r+1 s+1 1]'
"""
function vox2ras_0to1(M0::Matrix)

  if size(M0) != (4,4)
    error("Input must be a 4x4 matrix")
  end

  Q = zeros(4, 4)
  Q[1:3, 4] = ones(3, 1)

  M1 = inv(inv(M0)+Q)

  return M1
end

"""
    vox2ras_tkreg(voldim::Vector, voxres::Vector)

Return a 0-based vox2ras transform of a volume that is compatible with the
registration matrix produced by tkregister. May work with MINC xfm.

# Arguments
- voldim = [ncols;  nrows;  nslices ...]
- volres = [colres; rowres; sliceres ...]
"""
function vox2ras_tkreg(voldim::Vector, voxres::Vector)

  if length(voldim) < 3 | length(voxres) < 3
    error("Input vectors must have at least 3 elements")
  end

  T = zeros(4,4)
  T[4,4] = 1

  T[1,1] = -voxres[1]
  T[1,4] =  voxres[1] * voldim[1]/2

  T[2,3] =  voxres[3]
  T[2,4] = -voxres[3] * voldim[3]/2

  T[3,2] = -voxres[2]
  T[3,4] =  voxres[2] * voldim[2]/2

  return T
end

"""
    vox2ras_to_qform(vox2ras::Matrix)

Convert a vox2ras matrix to NIfTI qform parameters. The vox2ras should be 6 DOF.

Return the following NIfTI header fields:
- hdr.quatern_b
- hdr.quatern_c
- hdr.quatern_d
- hdr.qoffset_x
- hdr.qoffset_y
- hdr.qoffset_z
- hdr.pixdim[1]

From DG's vox2rasToQform.m:
  This code mostly just follows CH's mriToNiftiQform() in mriio.c
"""
function vox2ras_to_qform(vox2ras::Matrix)

  if size(vox2ras) != (4, 4)
    error("vox2ras size=" * string(size(vox2ras)) * ", must be (4, 4)")
  end

  x = vox2ras[1,4]
  y = vox2ras[2,4]
  z = vox2ras[3,4]

  d = sqrt.(sum(vox2ras[:, 1:3].^2; dims=1))
  Mdc = vox2ras[1:3, 1:3] ./ repeat(d, 3)
  if det(Mdc) == 0
    error("vox2ras determinant is 0")
  end

  r11 = Mdc[1,1]
  r21 = Mdc[2,1]
  r31 = Mdc[3,1]
  r12 = Mdc[1,2]
  r22 = Mdc[2,2]
  r32 = Mdc[3,2]
  r13 = Mdc[1,3]
  r23 = Mdc[2,3]
  r33 = Mdc[3,3]

  if det(Mdc) > 0
    qfac = 1.0
  else
    r13 = -r13
    r23 = -r23
    r33 = -r33
    qfac = -1.0
  end

  # From DG's vox2rasToQform.m: "following mat44_to_quatern()"
  a = r11 + r22 + r33 + 1.0
  if a > 0.5
    a = 0.5 * sqrt(a)
    b = 0.25 * (r32-r23) / a
    c = 0.25 * (r13-r31) / a
    d = 0.25 * (r21-r12) / a
  else
    xd = 1.0 + r11 - (r22+r33)
    yd = 1.0 + r22 - (r11+r33)
    zd = 1.0 + r33 - (r11+r22)
    if xd > 1
      b = 0.5 * sqrt(xd)
      c = 0.25 * (r12+r21) / b
      d = 0.25 * (r13+r31) / b
      a = 0.25 * (r32-r23) / b
    elseif yd > 1
      c = 0.5 * sqrt(yd)
      b = 0.25 * (r12+r21) / c
      d = 0.25 * (r23+r32) / c
      a = 0.25 * (r13-r31) / c
    else
      d = 0.5 * sqrt(zd)
      b = 0.25 * (r13+r31) / d
      c = 0.25 * (r23+r32) / d
      a = 0.25 * (r21-r12) / d
    end
    if a < 0
      a = -a
      b = -b
      c = -c
      d = -d
    end
  end

  return [b, c, d, x, y, z, qfac]
end

"""
    vox2ras_to_orient(vox2ras::Matrix)

Convert a vox2ras matrix to a 3-character array indicating image orientation,
e.g., RAS, LIA, etc.
"""
function vox2ras_to_orient(vox2ras::Matrix)

  orient = Vector{Char}(undef, 3)

  for idim in 1:3
    (amax, imax) = findmax(abs.(vox2ras[1:3, idim]))
    if imax == 1
      if vox2ras[imax, idim] > 0
        orient[idim] = 'R'
      else
        orient[idim] = 'L'
      end
    elseif imax == 2
      if vox2ras[imax, idim] > 0
        orient[idim] = 'A'
      else
        orient[idim] = 'P'
      end
    else
      if vox2ras[imax, idim] > 0
        orient[idim] = 'S'
      else
        orient[idim] = 'I'
      end
    end
  end

  return orient
end

"""
    mri_filename(fstring::String, checkdisk::Bool=true)

Return a valid file name, file stem, and file extension, given a string
`fstring` that can be either a file name or a file stem.

Valid extensions are: mgh, mgz, nii, nii.gz. Thus a file name is expected to
have the form stem.{mgh, mgz, nii, nii.gz}.

If `fstring` is a file name, then the stem and extension are determined from
`fstring`.

If `fstring` is a file stem and `checkdisk` is true (default), then the
file name and extension are determined by searching for a file on disk
named `fstring`.{mgh, mgz, nii, nii.gz}, where the possible extensions are
searched for in this order. If no such file is found, then empty strings are
returned.
"""
function mri_filename(fstring::String, checkdisk::Bool=true)

  fname = ""
  fstem = ""
  fext   = ""

  # List of supported file extensions
  extlist = ["mgh", "mgz", "nii", "nii.gz"]

  idot = findlast(isequal('.'), fstring)

  if isnothing(idot) && checkdisk
    # Filename has no extension, check if file exists with a supported extension
    for ext in extlist
      name = fstring * '.' * ext

      if isfile(name)
        fname = name
        fstem = fstring
        fext = ext
      end
    end
  else
    # Filename has an extension, check if it is one of the supported ones
    ext = lowercase(fstring[(idot+1):end]);

    if cmp(ext, "gz") == 0
      idot = findprev(isequal('.'), fstring, idot-1)

      if !isnothing(idot)
        ext = lowercase(fstring[(idot+1):end])
      end
    end

    if any(cmp.(ext, extlist) .== 0)
      fname = fstring
      fstem = fstring[1:idot-1]
      fext = ext
    end
  end

  return fname, fstem, fext
end

"""
    mri_read(infile::String, headeronly::Bool=false, permuteflag::Bool=false)

Read an image volume from disk and return an `MRI` structure similar to the
FreeSurfer MRI struct defined in mri.h.

# Arguments
- `infile::String`: Path to the input, which can be:
  1. An MGH file, e.g., f.mgh or f.mgz
  2. A NIfTI file, e.g., f.nii or f.nii.gz
  3. A file stem, in which case the format and full file name are determined
     by finding a file on disk named `infile`.{mgh, mgz, nii, nii.gz}
  4. A Bruker scan directory, which is expected to contain the following files:
     method, acqp, pdata/1/reco, pdata/1/2dseq

- `headeronly::Bool=false`: If true, the pixel data are not read in.

- `permuteflag::Bool==false`: If true, the first two dimensions of the image
  volume are permuted in the .vol, .volsize, and .volres fields of the output
  structure (this is the default behavior of the MATLAB version). The
  `permuteflag` will not affect the vox2ras matrices, which always map indices
  in the (column, row, slice) convention.

# Output
In the `MRI` structure:
- Times are in ms and angles are in radians.

- `vox2ras0`: This field contains the vox2ras matrix that converts 0-based
  (column, row, slice) voxel indices to (x, y, z) coordinates. This is the also
  the matrix that mri_write() uses to derive all geometry information
  (e.g., direction cosines, voxel resolution, P0). Thus if any other geometry
  fields of the structure are modified, the change will not be reflected in
  the output volume.

- `vox2ras1`: This field contains the vox2ras matrix that converts 1-based
  (column, row, slice) voxel indices to (x, y, z) coordinates.

- `niftihdr`: If the input file is NIfTI, this field contains a `NIfTIheader`
  structure.
"""
function mri_read(infile::String, headeronly::Bool=false, permuteflag::Bool=false)

  if isdir(infile)					#------ Bruker -------#
    mri = load_bruker(infile)
  else
    (fname, fstem, fext) = mri_filename(infile)
    if isempty(fname)
      error("Cannot determine format of " * infile)
    end

    if any(fext .== ["mgh", "mgz"])			#-------- MGH --------#
      (vol, M, mr_parms, volsz) =
       load_mgh(fname, nothing, nothing, headeronly)

      mri = MRI(vol)

      mri.fspec = fname
      mri.pwd = pwd()

      if !isempty(mr_parms)
        (mri.tr, mri.flip_angle, mri.te, mri.ti) = mr_parms
      end

      if isempty(M)
        error("Loading " * fname * " as MGH")
      else
        mri.vox2ras0 = M
      end

      mri.volsize = volsz[1:3]
      mri.nframes = length(volsz) < 4 ? 1 : volsz[4]
    elseif any(fext .== ["nii", "nii.gz"])		#------- NIfTI -------#
      if headeronly
        hdr = load_nifti(fname, headeronly)
      else
        hdr, vol = load_nifti(fname, headeronly)

        isempty(vol) && error("Loading " * fname * " as NIfTI")
      end

      # Compatibility with MRIread.m:
      # When data have > 4 dims, put all data into dim 4.
      volsz = Int.(hdr.dim[2:end])
      volsz = volsz[findall(volsz.>0)]
      if ~headeronly
        if length(volsz) < 5
          mri = MRI(vol)
        else
          mri = MRI(reshape(vol, volsz[1], volsz[2], volsz[3], :))
        end
      end

      mri.fspec = fname
      mri.pwd = pwd()

      mri.niftihdr = hdr

      mri.tr = hdr.pixdim[5]	# Already in msec

      mri.flip_angle = mri.te = mri.ti = 0

      mri.vox2ras0 = hdr.vox2ras

      mri.volsize = collect(volsz[1:3])
      mri.nframes = length(volsz) < 4 ? 1 : volsz[4]
    else
      error("File extension " * fext * " not supported")
    end

    # Optional DWI tables -----------------------------------------------#

    bfile = fstem * ".bvals"

    if ~isfile(bfile)
      bfile = fstem * ".bval"

      if ~isfile(bfile)
        bfile = ""
      end
    end

    gfile = fstem * ".bvecs"

    if ~isfile(gfile)
      gfile = fstem * ".bvec"

      if ~isfile(gfile)
        gfile = ""
      end
    end

    if ~isempty(bfile) && ~isempty(gfile)
      (b, g) = mri_read_bfiles(bfile, gfile)

      if length(b) == mri.nframes
        mri.bval = b
        mri.bvec = g

        # Normalize gradient vectors
        mri.bvec = mri.bvec ./ sqrt.(sum(mri.bvec.^2, dims=2))
        mri.bvec[isnan.(mri.bvec)] .= Float32(0)
      end
    end
  end

  # Dimensions not redundant when using header only
  (mri.width, mri.height, mri.depth) = collect(mri.volsize)

  # Set additional header fields related to volume geometry
  mri_set_geometry!(mri)

  # Permute volume from col-row-slice to row-col-slice, if desired
  # (This is the default behavior in the MATLAB version, but not here)
  if permuteflag
    mri.vol = permutedims(mri.vol, [2; 1; 3:ndims(mri.vol)])
    mri.volsize = mri.volsize[[2,1,3]]
    mri.volres = mri.volres[[2,1,3]]
    mri.ispermuted = true
  end

  return mri
end

"""
    mri_set_geometry!(mri::MRI)

Set header fields related to volume geometry in an `MRI` structure.

These are redundant fields that can be derived from the `vox2ras0`,
`width`, `height`, `depth` fields. They are in the MRI struct defined
in mri.h, as well as the MATLAB version of this structure, so they have
been added here for completeness.

Note: mri_write() derives all geometry information (i.e., direction cosines,
voxel resolution, and P0 from vox2ras0. If any of the redundant geometry
elements below are modified, the change will not be reflected in the output
volume.
"""
function mri_set_geometry!(mri::MRI)

  mri.vox2ras = mri.vox2ras0

  mri.nvoxels = mri.width * mri.height * mri.depth	# Number of voxels
  mri.xsize = sqrt(sum(mri.vox2ras[:,1].^2))	# Col
  mri.ysize = sqrt(sum(mri.vox2ras[:,2].^2))	# Row
  mri.zsize = sqrt(sum(mri.vox2ras[:,3].^2))	# Slice

  mri.x_r = mri.vox2ras[1,1]/mri.xsize	# Col
  mri.x_a = mri.vox2ras[2,1]/mri.xsize
  mri.x_s = mri.vox2ras[3,1]/mri.xsize

  mri.y_r = mri.vox2ras[1,2]/mri.ysize	# Row
  mri.y_a = mri.vox2ras[2,2]/mri.ysize
  mri.y_s = mri.vox2ras[3,2]/mri.ysize

  mri.z_r = mri.vox2ras[1,3]/mri.zsize	# Slice
  mri.z_a = mri.vox2ras[2,3]/mri.zsize
  mri.z_s = mri.vox2ras[3,3]/mri.zsize

  ic = [mri.width/2; mri.height/2; mri.depth/2; 1]
  c = mri.vox2ras * ic
  mri.c_r = c[1]
  mri.c_a = c[2]
  mri.c_s = c[3]

  #--------------------------------------------------------------------#
  # These are in the MATLAB version of this structure for convenience

  # 1-based vox2ras: Good for doing transforms on julia indices
  mri.vox2ras1 = vox2ras_0to1(mri.vox2ras)

  # Matrix of direction cosines
  mri.Mdc = mri.vox2ras[1:3, 1:3] *
            Diagonal(1 ./ [mri.xsize, mri.ysize, mri.zsize])

  # Vector of voxel resolutions
  mri.volres = [mri.xsize, mri.ysize, mri.zsize]

  mri.tkrvox2ras = vox2ras_tkreg(mri.volsize, mri.volres)
end

"""
    load_bruker(indir::String, headeronly::Bool=false)

Read Bruker image data from disk and return an `MRI` structure similar to the
FreeSurfer MRI struct defined in mri.h.

# Arguments
- `indir::String`: Path to a Bruker scan directory (files called method, acqp,
pdata/1/reco, and pdata/1/2dseq are expected to be found in it)

- `headeronly::Bool=false`: If true, the pixel data are not read in.
"""
function load_bruker(indir::String, headeronly::Bool=false)

  dname = abspath(indir)
  methfile = dname * "/method"
  acqpfile = dname * "/acqp"
  recofile = dname * "/pdata/1/reco"
  imgfile  = dname * "/pdata/1/2dseq"

  if any(.!isfile.([methfile, acqpfile, recofile, imgfile]))
    error("Input directory must contain the files: " *
          "method, acqp, pdata/1/reco, pdata/1/2dseq")
  end

  mri = MRI(Array{Float32, 4}(undef, (0, 0, 0, 0)))

  mri.fspec = imgfile
  mri.pwd = pwd()

  nb0 = 0

  # Read information for the image header from the Bruker method file
  io = open(methfile, "r")

  while !eof(io)
    ln = readline(io)

    if startswith(ln, "##\$PVM_SpatResol=")		# Voxel size
      ln = readline(io)
      ln = split(ln)
      mri.volres = parse.(Float32, ln)
    elseif startswith(ln, "##\$PVM_Matrix=")		# Matrix size
      ln = readline(io)
      ln = split(ln)
      mri.volsize = parse.(Float32, ln)
    elseif startswith(ln, "##\$EchoTime=")		# TE
      ln = split(ln, "=")[2]
      mri.te = parse(Float32, ln)
    elseif startswith(ln, "##\$PVM_RepetitionTime=")	# TR
      ln = split(ln, "=")[2]
      mri.tr = parse(Float32, ln)
    elseif startswith(ln, "##\$PVM_DwAoImages=")	# Number of b=0 volumes
      ln = split(ln, "=")[2]
      nb0 = parse(Int64, ln)
    elseif startswith(ln, "##\$PVM_DwDir=")		# Diffusion gradients
      nval = split(ln, "(")[2]
      nval = split(nval, ")")[1]
      nval = split(nval, ",")
      nval = prod(parse.(Int64, nval))

      nread = 0
      bvec = Vector{Float32}(undef, 0)

      while nread < nval
        ln = readline(io)
        ln = split(ln)

        nread += length(ln)
	push!(bvec, parse.(Float32, ln)...)
      end

      mri.bvec = permutedims(reshape(bvec, 3, :), [2 1])

      # Normalize gradient vectors
      mri.bvec = mri.bvec ./ sqrt.(sum(mri.bvec.^2, dims=2))
      mri.bvec[isnan.(mri.bvec)] .= Float32(0)
    elseif startswith(ln, "##\$PVM_DwEffBval=")		# b-values
      nval = split(ln, "(")[2]
      nval = split(nval, ")")[1]
      nval = parse(Int64, nval)

      nread = 0
      bval = Vector{Float32}(undef, 0)

      while nread < nval
        ln = readline(io)
        ln = split(ln)

        nread += length(ln)
	push!(bval, parse.(Float32, ln)...)
      end

      mri.bval = bval
    end
  end

  close(io)

  # Add b=0 volumes to the gradient table
  # (The method file includes them in the list of b-values but not vectors)
  for ib0 in 1:nb0
    mri.bvec = vcat([0 0 0], mri.bvec)
  end

  # Read receiver gain from Bruker acqp file
  io = open(acqpfile, "r")

  gain = Float32(1)

  while !eof(io)
    ln = readline(io)

    if startswith(ln, "##\$RG=")			# Receiver gain
      ln = split(ln, "=")[2]
      gain = parse(Float32, ln)
    end
  end

  gain /= 64

  close(io)

  # Read information about image binary data from Bruker reco file
  io = open(recofile, "r")

  data_type = Int32
  int_offset = Vector{Float32}(undef, 0)
  int_slope = Vector{Float32}(undef, 0)
  byte_order = ""

  while !eof(io)
    ln = readline(io)

    if startswith(ln, "##\$RECO_wordtype=")		# Bytes per voxel
      ln = split(ln, "=")[2]

      if ln == "_32BIT_SGN_INT"
        data_type = Int32
      elseif ln == "_16BIT_SGN_INT"
        data_type = Int16
      elseif ln == "_8BIT_SGN_INT"
        data_type = Int8
      end
    elseif startswith(ln, "##\$RECO_map_offset=")	# Intensity offset
      nval = split(ln, "(")[2]
      nval = split(nval, ")")[1]
      nval = parse(Int64, nval)

      nread = 0

      while nread < nval
        ln = readline(io)
        ln = split(ln)

        nread += length(ln)
        push!(int_offset, parse.(Float32, ln)...)
      end
    elseif startswith(ln, "##\$RECO_map_slope")		# Intensity slope
      nval = split(ln, "(")[2]
      nval = split(nval, ")")[1]
      nval = parse(Int64, nval)

      nread = 0

      while nread < nval
        ln = readline(io)
        ln = split(ln)

        nread += length(ln)
        push!(int_slope, parse.(Float32, ln)...)
      end
    elseif startswith(ln, "##\$RECO_byte_order=")	# Byte order
      byte_order = split(ln, "=")[2]
    end
  end

  close(io)

  mri.nframes = length(int_slope)

  mri.vox2ras0 = Matrix(Diagonal([mri.volres; 1]))

  if headeronly
    return mri
  end

  # Read image data
  io = open(imgfile, "r")

  vol = read!(io, Array{data_type}(undef, (mri.volsize..., mri.nframes)))

  close(io)

  if byte_order == "littleEndian"
    vol = ltoh.(vol)
  else
    vol = ntoh.(vol)
  end

  # Apply intensity offset and slope
  vol = Int32.(vol)

  if data_type == Float32
    mri.vol = Float32.(vol)
  else
    mri.vol = Array{Float32}(undef, size(vol))
    for iframe in 1:mri.nframes
      mri.vol[:,:,:,iframe] = vol[:,:,:,iframe] ./ int_slope[iframe] .+
                                                   int_offset[iframe]
    end
  end

  # Normalize by receiver gain
  mri.vol = mri.vol / gain

  return mri
end

"""
    load_mgh(fname::String, slices::Union{Vector{Unsigned}, Nothing}=nothing, frames::Union{Vector{Unsigned}, Nothing}=nothing, headeronly::Bool=false)

Load a .mgh or .mgz file from disk.

Return:
- The image data as an array

- The 4x4 vox2ras matrix that transforms 0-based voxel indices to coordinates
  such that: [x y z]' = M * [i j k 1]'

- MRI parameters as a vector: [tr, flip_angle, te, ti]

- The image volume dimensions as a vector (useful when `headeronly` is true,
  in which case the image array will be empty)

# Arguments
- `fname::String`: Path to the .mgh/.mgz file

- `slices::Vector`: 1-based slice numbers to load (default: read all slices)

- `frames::Vector`: 1-based frame numbers to load (default: read all frames)

- `headeronly::Bool=false`: If true, the pixel data are not read in.
"""
function load_mgh(fname::String, slices::Union{Vector{Unsigned}, Nothing}=nothing, frames::Union{Vector{Unsigned}, Nothing}=nothing, headeronly::Bool=false)

  vol = M = mr_parms = volsz = []

  # Unzip if it is compressed
  ext = lowercase(fname[(end-1):end])

  if cmp(ext, "gz") == 0
    # Create unique temporary file name
    tmpfile = tempname(get_tmp_path()) * ".load_mgh.mgh"
    gzipped = true

    if Sys.isapple()
      cmd = `gunzip -c $fname`
    else
      cmd = `zcat $fname`
    end
    run(pipeline(cmd, stdout=tmpfile))
  else
    tmpfile = fname
    gzipped = false
  end

  io = open(tmpfile, "r")

  v       = ntoh(read(io, Int32))
  ndim1   = ntoh(read(io, Int32))
  ndim2   = ntoh(read(io, Int32))
  ndim3   = ntoh(read(io, Int32))
  nframes = ntoh(read(io, Int32))
  type    = ntoh(read(io, Int32))
  dof     = ntoh(read(io, Int32))

  if !isnothing(slices) && any(slices .> ndim3)
    error("Some slices=" * string(slices) * " exceed nslices=" * string(dim3))
  end

  if !isnothing(frames) && any(frames .> nframes)
    error("Some frames=" * string(frames) * " exceed nframes=" * string(nframes))
  end

  UNUSED_SPACE_SIZE = 256
  USED_SPACE_SIZE = 3*4+4*3*4		# Space for RAS transform

  unused_space_size = UNUSED_SPACE_SIZE-2
  ras_good_flag = ntoh(read(io, Int16))

  if ras_good_flag > 0
    delta  = ntoh.(read!(io, Vector{Float32}(undef, 3)))
    Mdc    = ntoh.(read!(io, Vector{Float32}(undef, 9)))
    Mdc    = reshape(Mdc, (3,3))
    Pxyz_c = ntoh.(read!(io, Vector{Float32}(undef, 3)))

    D = Diagonal(delta)

    Pcrs_c = [ndim1; ndim2; ndim3]/2	# Should this be kept?

    Pxyz_0 = Pxyz_c - Mdc*D*Pcrs_c

    M = [Mdc*D Pxyz_0; 0 0 0 1]
    ras_xform = [Mdc Pxyz_c; 0 0 0 1]
    unused_space_size = unused_space_size - USED_SPACE_SIZE
  end

  skip(io, unused_space_size)
  nv = ndim1 * ndim2 * ndim3 * nframes
  volsz = [ndim1 ndim2 ndim3 nframes]

  MRI_UCHAR  =  0
  MRI_INT    =  1
  MRI_LONG   =  2
  MRI_FLOAT  =  3
  MRI_SHORT  =  4
  MRI_BITMAP =  5
  MRI_USHRT  = 10

  # Determine number of bytes per voxel
  if type == MRI_FLOAT
    nbytespervox = 4
    dtype = Float32
  elseif type == MRI_UCHAR
    nbytespervox = 1
    dtype = UInt8
  elseif type == MRI_SHORT
    nbytespervox = 2
    dtype = Int16
  elseif type == MRI_USHRT
    nbytespervox = 2
    dtype = UInt16
  elseif type == MRI_INT
    nbytespervox = 4
    dtype = Int32
  end

  if headeronly
    skip(io, nv*nbytespervox)
    if !eof(io)
      mr_parms = ntoh.(read!(io, Vector{Float32}(undef, 4)))
    end
    close(io)

    if gzipped		# Clean up
      cmd = `rm -f $tmpfile`
      run(cmd)
    end

    return vol, M, mr_parms, volsz
  end

  if isnothing(slices) && isnothing(frames)		# Read the entire volume
    vol = ntoh.(read!(io, Array{dtype}(undef, Tuple(volsz))))
  else						# Read a subset of slices/frames
    isnothing(frames) && (frames = 1:nframes)
    isnothing(slices) && (slices = 1:ndim3)

    nvslice = ndim1 * ndim2
    nvvol = nvslice * ndim3
    filepos0 = position(io)

    vol = zeros(dtype, ndim1, ndim2, length(slices), length(frames))

    for iframe in 1:length(frames)
      frame = frames(iframe)

      for islice in 1:length(slices)
        slice = slices(islice)

        filepos = ((frame-1)*nvvol + (slice-1)*nvslice)*nbytespervox + filepos0
        seek(io, filepos)

        vol[:, :, islice, iframe] =
          ntoh.(read!(io, Array{dtype}(undef, Tuple(volsz[1:2]))))
      end
    end

    # Seek to just beyond the last slice/frame
    filepos = (nframes*nvvol)*nbytespervox + filepos0;
    seek(io, filepos)
  end

  if !eof(io)
    mr_parms = ntoh.(read!(io, Vector{Float32}(undef, 4)))
  end

  close(io)

  if gzipped		# Clean up
    cmd = `rm -f $tmpfile`
    run(cmd)
  end

  return vol, M, mr_parms, volsz
end

"""
    load_nifti_hdr(fname::String) -> NIfTIheader

Load the header of a .nii volume from disk.

Return a `NIfTIheader` structure, where units have been converted to mm and
msec, the sform and qform matrices have been computed and stored in the .sform
and .qform fields, and the .vox2ras field has been set to the sform (if valid),
then the qform (if valid).

Assume that the input file is uncompressed (compression is handled in the
wrapper load_nifti()).

Handle data with more than 32k cols by looking at the .dim field of the header.
If dim[2] = -1, then the .glmin field contains the number of columns. This is
FreeSurfer specific, for handling surfaces. When the total number of spatial
voxels equals 163842, then reshape the volume to 163842x1x1xnframes. This is
for handling the 7th order icosahedron used by FS group analysis.
"""
function load_nifti_hdr(fname::String)

  nifti_header_size = 348

  # Read header as vector of bytes
  io = open(fname, "r")
  buffer = read!(io, Vector{UInt8}(undef, nifti_header_size))
  close(io)

  typelist = NIfTIheader.types[1:end-4]
  plist = pointer(buffer) .+ cumsum(sizeof.(typelist))

  # Check endian-ness of first entry to decide if byte order should be reversed
  headsize = unsafe_load(convert(Ptr{typelist[1]}, pointer(buffer)))

  if headsize == nifti_header_size
    do_bswap = false
    f = identity
  elseif headsize == bswap(Int32(nifti_header_size))
    do_bswap = true
    f = bswap
  else
    error("Invalid header size " * string(headsize) * " found in NIfTI header")
  end

  # Convert vector of bytes to array of header fields
  header = f.([(unsafe_load(convert(Ptr{typelist[i]}, plist[i-1])))
               for i = 2:length(typelist)])

  pushfirst!(header, f.(headsize))

  # "dim" field:
  # This is to accomodate structures with more than 32k cols
  # FreeSurfer specific (see also mriio.c)
  i_dim = findfirst(fieldnames(NIfTIheader) .== :dim)

  if header[i_dim][2] < 0
    i_glmin = findfirst(fieldnames(NIfTIheader) .== :glmin)
    header[i_dim] = (header[i_dim][1], header[i_glmin], header[i_dim][3:end]...)
    header[i_glmin] = 0
  end

  nspatial = prod(Int32.(header[i_dim][2:4]))
  if nspatial == 163842		# Ico7 surface
    header[i_dim] = (header[i_dim][1], 163842, 1, 1, header[i_dim][5:end]...)
  end

  # "xyz_units" field:
  # Find physical and time units of header fields
  i_xyzt_units = findfirst(fieldnames(NIfTIheader) .== :xyzt_units)

  xyzunits = header[i_xyzt_units] & Int8(7)	# Bitwise AND with 00000111
  if xyzunits == 1
    xyzscale = Float32(1000)	# meters
  elseif xyzunits == 2
    xyzscale = Float32(1)	# mm
  elseif xyzunits == 3
    xyzscale = Float32(.001)	# microns
  else
    println("WARNING: xyz units code " * string(xyzunits) *
            " is unrecognized, assuming mm")
    xyzscale = Float32(1)	# just assume mm
  end

  tunits = header[i_xyzt_units] & Int8(56)	# Bitwise AND with 00111000
  if tunits == 8
    tscale = Float32(1000)	# seconds
  elseif tunits == 16
    tscale = Float32(1)		# msec
  elseif tunits == 32
    tscale = Float32(.001)	# microsec
  else
    tscale = Float32(0)		# no time scale
  end

  # "pixdim", "srow_x", "srow_y", "srow_z" fields:
  # Convert to physical units to mm and time units to msec
  i_pixdim = findfirst(fieldnames(NIfTIheader) .== :pixdim)
  header[i_pixdim] = (header[i_pixdim][1],
                      header[i_pixdim][2:4] .* xyzscale...,
                      header[i_pixdim][5] * tscale,
                      header[i_pixdim][6:8]...)

  i_srow_x = findfirst(fieldnames(NIfTIheader) .== :srow_x)
  header[i_srow_x] = header[i_srow_x] .* xyzscale

  i_srow_y = findfirst(fieldnames(NIfTIheader) .== :srow_y)
  header[i_srow_y] = header[i_srow_y] .* xyzscale

  i_srow_z = findfirst(fieldnames(NIfTIheader) .== :srow_z)
  header[i_srow_z] = header[i_srow_z] .* xyzscale

  # Change value in xyzt_units to reflect scale change
  header[i_xyzt_units] = Int8(2) | Int8(16)	# Bitwise OR of 2=mm, 16=msec

  # Sform matrix
  sform = [collect(header[i_srow_x])';
	   collect(header[i_srow_y])';
	   collect(header[i_srow_z])';
	   Float32.([0 0 0 1])]

  # Qform matrix
  # (From DG's load_nifti_hdr.m: not quite sure how all this works,
  # mainly just copied CH's code from mriio.c)
  i_quatern_b = findfirst(fieldnames(NIfTIheader) .== :quatern_b)
  i_quatern_c = findfirst(fieldnames(NIfTIheader) .== :quatern_c)
  i_quatern_d = findfirst(fieldnames(NIfTIheader) .== :quatern_d)
  i_quatern_x = findfirst(fieldnames(NIfTIheader) .== :quatern_x)
  i_quatern_y = findfirst(fieldnames(NIfTIheader) .== :quatern_y)
  i_quatern_z = findfirst(fieldnames(NIfTIheader) .== :quatern_z)
  b = header[i_quatern_b]
  c = header[i_quatern_c]
  d = header[i_quatern_d]
  x = header[i_quatern_x]
  y = header[i_quatern_y]
  z = header[i_quatern_z]
  a = Float32(1) - (b*b + c*c + d*d)
  if abs(a) < 1.0e-7
    a = Float32(1) / sqrt(b*b + c*c + d*d)
    b = b*a
    c = c*a
    d = d*a
    a = Float32(0)
  else
    a = sqrt(a)
  end
  r11 = a*a + b*b - c*c - d*d
  r12 = 2*b*c - 2*a*d
  r13 = 2*b*d + 2*a*c
  r21 = 2*b*c + 2*a*d
  r22 = a*a + c*c - b*b - d*d
  r23 = 2*c*d - 2*a*b
  r31 = 2*b*d - 2*a*c
  r32 = 2*c*d + 2*a*b
  r33 = a*a + d*d - c*c - b*b
  if header[i_pixdim][1] < 0.0
    r13 = -r13
    r23 = -r23
    r33 = -r33
  end
  qMdc = [r11 r12 r13; r21 r22 r23; r31 r32 r33]
  D = Diagonal(collect(header[i_pixdim][2:4]))
  P0 = [x y z]'
  qform = [qMdc*D P0; 0 0 0 1]

  # "sform_code", "qform_code" fields:
  # Determine which matrix to use as vox2ras
  i_sform_code = findfirst(fieldnames(NIfTIheader) .== :sform_code)
  i_qform_code = findfirst(fieldnames(NIfTIheader) .== :qform_code)

  if header[i_sform_code] != 0
    # Use sform first
    vox2ras = sform
  elseif header[i_qform_code] != 0
    # Then use qform first
    vox2ras = qform
  else
    println("WARNING: neither sform or qform are valid in " * fname)
    D = Diagonal(collect(header[i_pixdim][2:4]))
    P0 = [0 0 0]'
    vox2ras = [D P0; 0 0 0 1]
  end

  return NIfTIheader(header..., do_bswap, sform, qform, vox2ras)
end

"""
    load_nifti(fname::String, hdronly::Bool=false) -> NIfTIheader, Array

Load a NIfTI (.nii or .nii.gz) volume from disk and return an array containing
the image data and a `NIfTIheader` structure and the image data as an array.

Handle compressed NIfTI (nii.gz) by issuing an external Unix call to
uncompress the file to a temporary file, which is then deleted.

The output `NIfTIheader` structure contains:
- the units for each dimension of the volume [mm or msec], in the .pixdim field
- the sform and qform matrices, in the .sform and .qform fields
- the vox2ras matrix, which is the sform (if valid), otherwise the qform, in
  the .vox2ras field
"""
function load_nifti(fname::String, hdronly::Bool=false)

  # Unzip if it is compressed
  ext = lowercase(fname[(end-1):end])

  if cmp(ext, "gz") == 0
    # Create unique temporary file name
    tmpfile = tempname(get_tmp_path()) * ".load_nifti.nii"
    gzipped = true

    if Sys.isapple()
      cmd = `gunzip -c $fname`
    else
      cmd = `zcat $fname`
    end
    run(pipeline(cmd, stdout=tmpfile))
  else
    tmpfile = fname
    gzipped = false
  end

  # Read NIfTI header
  hdr = load_nifti_hdr(tmpfile)

  if isempty(hdr.vox2ras)
    if gzipped		# Clean up
      cmd = `rm -f $tmpfile`
      run(cmd)
    end
    return hdr
  end

  # If only header is desired, return now
  if hdronly
    if gzipped		# Clean up
      cmd = `rm -f $tmpfile`
      run(cmd)
    end
    return hdr
  end

  # Get volume dimensions
  dim = hdr.dim[2:findlast(hdr.dim .!= 0)]

  # Open to read the pixel data
  io = open(tmpfile, "r")

  # Get past the header
  seek(io, Int64(round(hdr.vox_offset)))

  if hdr.datatype == 2
    dtype = UInt8
  elseif hdr.datatype == 4
    dtype = Int16
  elseif hdr.datatype == 8
    dtype = Int32
  elseif hdr.datatype == 16
    dtype = Float32
  elseif hdr.datatype == 64
    dtype = Float64
  elseif hdr.datatype == 256
    dtype = Int8
  elseif hdr.datatype == 512
    dtype = UInt16
  elseif hdr.datatype == 768
    dtype = UInt32
  else
    close(io)
    if gzipped		# Clean up
      cmd = `rm -f $tmpfile`
      run(cmd)
    end
    error("Data type " * string(hdr.datatype) * " not supported")
  end

  vol = read!(io, Array{dtype}(undef, dim))

  close(io)
  if gzipped		# Clean up
    cmd = `rm -f $tmpfile`
    run(cmd)
  end

  # Check if end-of-file was reached
  if !eof(io)
    error(tmpfile * ", read a " * string(size(vol)) *
          " volume but did not reach end of file")
  end

  # If needed, reverse order of bytes to correct endian-ness
  if hdr.do_bswap
    vol .= bswap.(vol)
  end

  if hdr.scl_slope != 0 && !(hdr.scl_inter == 0 && hdr.scl_slope == 1)
    # Rescaling is not needed if the slope==1 and intersect==0,
    # skipping this preserves the numeric class of the data
    vol .= dtype.(vol .* hdr.scl_slope .+ hdr.scl_inter)
  end

  return hdr, vol
end

"""
    mri_write(mri::MRI, outfile::String, datatype::DataType=Float32)

Write an MRI volume to disk. Return true is an error occurred (i.e., the
number of bytes written were not as expected based on the size of the volume).

# Arguments
- `mri::MRI`: A structure like that returned by `mri_read()`. The geometry
  (i.e., direction cosines, voxel resolution, and P0) are all recomputed from
  mri.vox2ras0. So, if a method has changed one of the other fields, e.g.,
  mri.x_r, this change will not be reflected in the output volume.

- `outfile::String`: Path to the output file, which can be:
  1. An MGH file, e.g., f.mgh or f.mgz (uncompressed or compressed)
  2. A NIfTI file, e.g., f.nii or f.nii.gz (uncompressed or compressed).

- `datatype::DataType=eltype(mri.vol)`: Only applies to NIfTI and can be UInt8,
  Int16, Int32, Float32, Float64, Int8, UInt16, UInt32. By default, the native
  data type of the volume array is preserved when writing to disk.
"""
function mri_write(mri::MRI, outfile::String, datatype::DataType=eltype(mri.vol))

  err = true

  if isempty(mri.vol)
    error("Input structure has empty vol field")
  end

  vsz = collect(size(mri.vol))
  nvsz = length(vsz)
  if nvsz < 4
    vsz = [vsz; ones(eltype(vsz), 4-nvsz)]
  end

  if isempty(mri.volsize)
    mri.volsize = vsz[1:3]
  end

  if mri.nframes == 0
    mri.nframes = vsz[4]
  end

  if isempty(mri.vox2ras0)
    mri.vox2ras0 = Diagonal(ones(4))
  end

  if isempty(mri.volres)
    mri.volres = sqrt.((ones(3)' * (mri.vox2ras0[1:3,1:3].^2))')
  end

  (fname, fstem, fext) = mri_filename(outfile, false)	# false = no checkdisk
  if isempty(fname)
    error("Cannot determine format of " * outfile)
  end

  if any(cmp.(fext, ["mgh", "mgz"]) .== 0)		#-------- MGH --------#
    M = mri.vox2ras0
    mr_parms = [mri.tr, mri.flip_angle, mri.te, mri.ti]

    if mri.ispermuted
      err = save_mgh(permutedims(mri.vol, [2; 1; 3:ndims(mri.vol)]),
                     fname, M, mr_parms)
    else
      err = save_mgh(mri.vol, fname, M, mr_parms)
    end
  elseif any(cmp.(fext, ["nii", "nii.gz"]) .== 0)	#------- NIfTI -------#
    hdr_sizeof_hdr    = Int32(348)
    hdr_data_type     = Tuple(zeros(UInt8, 10))
    hdr_db_name       = Tuple(zeros(UInt8, 18))
    hdr_extents       = Int32(0)
    hdr_session_error = Int16(0)
    hdr_regular       = UInt8(0)
    hdr_dim_info      = UInt8(0)

    dim = ones(Int16, 8)
    dim[1] = mri.nframes > 1 ? 4 : 3
    dim[2:4] = mri.ispermuted ? mri.volsize[[2,1,3]] : mri.volsize[1:3]
    dim[5] = mri.nframes

    # This is to accomodate structures with more than 32k cols
    # FreeSurfer specific. See also mriio.c.
    if dim[2] > 2^15
      hdr_glmin = dim[2]
      dim[2] = -1
    end

    hdr_dim = Tuple(dim)

    hdr_intent_p1   = Float32(0)
    hdr_intent_p2   = Float32(0)
    hdr_intent_p3   = Float32(0)
    hdr_intent_code = Int16(0)

    if datatype == UInt8
      hdr_datatype = Int16(2)
      hdr_bitpix   = Int16(8)
    elseif datatype == Int16
      hdr_datatype = Int16(4)
      hdr_bitpix   = Int16(16)
    elseif datatype == Int32
      hdr_datatype = Int16(8)
      hdr_bitpix   = Int16(32)
    elseif datatype == Float32
      hdr_datatype = Int16(16)
      hdr_bitpix   = Int16(32)
    elseif datatype == Float64
      hdr_datatype = Int16(64)
      hdr_bitpix   = Int16(64)
    elseif datatype == Int8
      hdr_datatype = Int16(256)
      hdr_bitpix   = Int16(8)
    elseif datatype == UInt16
      hdr_datatype = Int16(512)
      hdr_bitpix   = Int16(16)
    elseif datatype == UInt32
      hdr_datatype = Int16(768)
      hdr_bitpix   = Int16(32)
    else
      error("Data type " * string(datatype) * " not supported")
    end

    hdr_slice_start = Int16(0)

    if mri.ispermuted
      hdr_pixdim = Float32.((0, mri.volres[[2,1,3]]..., mri.tr, 0, 0, 0))
    else
      hdr_pixdim = Float32.((0, mri.volres[1:3]..., mri.tr, 0, 0, 0))
    end

    hdr_vox_offset = Float32(352)
    hdr_scl_slope  = mri.niftihdr.scl_slope
    hdr_scl_inter  = mri.niftihdr.scl_inter
    hdr_slice_end  = Int16(0)
    hdr_slice_code = Int8(0)

    hdr_xyzt_units = Int8(2) | Int8(16)   # Bitwise OR of 2=mm, 16=msec

    hdr_cal_max        = maximum(mri.vol)
    hdr_cal_min        = minimum(mri.vol)
    hdr_slice_duration = Float32(0)
    hdr_toffset        = Float32(0)
    hdr_glmax          = Int32(0)
    if dim[2] != -1
      hdr_glmin        = Int32(0)
    end
    hdr_descrip        = UInt8.(Tuple(@sprintf("%-80s","FreeSurfer julia")))
    hdr_aux_file       = Tuple(zeros(UInt8, 24))
    hdr_qform_code     = Int16(1)		# 1=NIFTI_XFORM_SCANNER_ANAT
    hdr_sform_code     = Int16(1)		# 1=NIFTI_XFORM_SCANNER_ANAT

    # Qform (must have 6 DOF)
    (b, c, d, x, y, z, qfac) = vox2ras_to_qform(mri.vox2ras0)
    hdr_pixdim    = (Float32(qfac), hdr_pixdim[2:end]...)
    hdr_quatern_b = Float32(b)
    hdr_quatern_c = Float32(c)
    hdr_quatern_d = Float32(d)
    hdr_quatern_x = Float32(x)
    hdr_quatern_y = Float32(y)
    hdr_quatern_z = Float32(z)

    # Sform (can be any affine)
    hdr_srow_x = Tuple(Float32.(mri.vox2ras0[1,:]))
    hdr_srow_y = Tuple(Float32.(mri.vox2ras0[2,:]))
    hdr_srow_z = Tuple(Float32.(mri.vox2ras0[3,:]))

    hdr_intent_name = (UInt8.(collect("huh?"))..., zeros(UInt8, 12)...)
    hdr_magic       = Tuple(collect("n+1\0"))

    hdr = NIfTIheader( hdr_sizeof_hdr,
                       hdr_data_type,
                       hdr_db_name,
                       hdr_extents,
                       hdr_session_error,
                       hdr_regular,
                       hdr_dim_info,
                       hdr_dim,
                       hdr_intent_p1,
                       hdr_intent_p2,
                       hdr_intent_p3,
                       hdr_intent_code,
                       hdr_datatype,
                       hdr_bitpix,
                       hdr_slice_start,
                       hdr_pixdim,
                       hdr_vox_offset,
                       hdr_scl_slope,
                       hdr_scl_inter,
                       hdr_slice_end,
                       hdr_slice_code,
                       hdr_xyzt_units,
                       hdr_cal_max,
                       hdr_cal_min,
                       hdr_slice_duration,
                       hdr_toffset,
                       hdr_glmax,
                       hdr_glmin,
                       hdr_descrip,
                       hdr_aux_file,
                       hdr_qform_code,
                       hdr_sform_code,
                       hdr_quatern_b,
                       hdr_quatern_c,
                       hdr_quatern_d,
                       hdr_quatern_x,
                       hdr_quatern_y,
                       hdr_quatern_z,
                       hdr_srow_x,
                       hdr_srow_y,
                       hdr_srow_z,
                       hdr_intent_name,
                       hdr_magic,
                       UInt8(0),
                       Matrix{Float32}(undef, 0, 0),
                       Matrix{Float32}(undef, 0, 0),
                       Matrix{Float32}(undef, 0, 0) )

    if mri.ispermuted
      vol = permutedims(mri.vol, [2; 1; 3:ndims(mri.vol)])
    else
      vol = mri.vol
    end

    err = save_nifti(hdr, vol, fname)
  else
    error("File extension " * fext * " not supported")
  end

  if err
    println("WARNING: Problem saving " * outfile)
  end

  # Optional DWI tables -----------------------------------------------#

  if !isempty(mri.bval)
    bfile = fstem * ".bvals"
    writedlm(bfile, mri.bval, ' ')
  end

  if !isempty(mri.bvec)
    gfile = fstem * ".bvecs"
    writedlm(gfile, mri.bvec, ' ')
  end

  return err
end

"""
    save_mgh(vol::Array, fname::String, M::Matrix=Diagonal(ones(4)), mr_parms::Vector=zeros(4))

Write an MRI volume to a .mgh or .mgz file. Return true is an error occurred
(i.e., the number of bytes written were not as expected based on the size of
the volume).

# Arguments
- `vol::Array`: the image data

- `fname::String`: path to the output file

- `M::Matrix`: the 4x4 vox2ras transform such that [x y z]' = M * [i j k 1]',
  where the voxel indices (i, j, k) are 0-based.

- `mr_parms::Vector`: a vector of MRI parameters, [tr, flip_angle, te, ti]
"""
function save_mgh(vol::Array, fname::String, M::Matrix=Diagonal(ones(4)), mr_parms::Vector=zeros(4))

  if size(M) != (4, 4)
    error("M size=" * string(size(M)) * ", must be (4, 4)")
  end

  if length(mr_parms) != 4
    error("mr_parms length=" * string(length(mr_parms)) * ", must be 4")
  end

  (ndim1, ndim2, ndim3, frames) = size(vol)

  dtype = eltype(vol)

  # The MATLAB version ingores these and always writes the data as float
  # Here we use them properly
  MRI_UCHAR =  0
  MRI_INT =    1
  MRI_LONG =   2
  MRI_FLOAT =  3
  MRI_SHORT =  4
  MRI_BITMAP = 5
  MRI_TENSOR = 6
  MRI_USHRT  = 10

  # Determine number of bytes per voxel
  if dtype == Float32
    type = MRI_FLOAT
  elseif dtype == UInt8
    type = MRI_UCHAR
  elseif dtype == Int32
    type = MRI_INT
  elseif dtype == Int64
    type = MRI_LONG
  elseif dtype == Int16
    type = MRI_SHORT
  elseif dtype == UInt16
    type = MRI_USHRT
  end

  io = open(fname, "w")

  nb = 0

  # Write everything as big-endian
  nb += write(io, hton(Int32(1)))		# magic number
  nb += write(io, hton(Int32(ndim1)))
  nb += write(io, hton(Int32(ndim2)))
  nb += write(io, hton(Int32(ndim3)))
  nb += write(io, hton(Int32(frames)))
  nb += write(io, hton(Int32(type)))

  nb += write(io, hton(Int32(1)))		# dof (not used)

  UNUSED_SPACE_SIZE = 256
  USED_SPACE_SIZE = (3*4+4*3*4)		# Space for RAS transform

  MdcD = M[1:3,1:3]

  delta = sqrt.(sum(MdcD.^2; dims=1))
  Mdc = MdcD ./ repeat(delta, 3)

  Pcrs_c = [ndim1/2, ndim2/2, ndim3/2, 1]
  Pxyz_c = M*Pcrs_c
  Pxyz_c = Pxyz_c[1:3]

  nb += write(io, hton(Int16(1)))		# ras_good_flag = 1
  nb += write(io, hton.(Float32.(delta)))
  nb += write(io, hton.(Float32.(Mdc)))
  nb += write(io, hton.(Float32.(Pxyz_c)))

  unused_space_size = UNUSED_SPACE_SIZE-2
  unused_space_size = unused_space_size - USED_SPACE_SIZE
  nb += write(io, UInt8.(zeros(unused_space_size)))

  nb += write(io, hton.(vol))

  nb += write(io, hton.(Float32.(mr_parms)))

  close(io)

  err = (nb != (sizeof(Int32) * 7 +
                sizeof(Int16) +
                sizeof(UInt8) * unused_space_size +
                sizeof(Float32) * 19 +
                sizeof(dtype) * length(vol)))

  ext = lowercase(fname[(end-1):end])

  if cmp(ext, "gz") == 0
    cmd = `gzip -f $fname`
    run(cmd)
    cmd = `mv $fname.gz $fname`
    run(cmd)
  end

  return err
end

"""
    save_nifti(hdr::NIfTIheader, vol::Array{T}, fname::String) where T<:Number

Write an MRI volume to a .nii or .nii.gz file. Return true is an error occurred
(i.e., the number of bytes written were not as expected based on the size of
the volume).

# Arguments
- `hdr::NIfTIheader`: a NIfTI header structure

- `vol::Array{T}`: an array that contains the image data

- `fname::String`: path to the output file

Handle data structures with more than 32k cols by setting hdr.dim[2] = -1 and
hdr.glmin = ncols. This is FreeSurfer specific, for handling surfaces. The
exception to this is when the total number of spatial voxels equals 163842,
then the volume is reshaped to 27307x1x6xnframes. This is for handling the 7th
order icosahedron used by FS group analysis.
"""
function save_nifti(hdr::NIfTIheader, vol::Array{T}, fname::String) where T<:Number

  ext = lowercase(fname[(end-1):end])
  if cmp(ext, "gz") == 0
    gzip_needed = true
    fname = fname[1:(end-3)]
  else
    gzip_needed = false
  end

  # Check for ico7
  sz = size(vol)
  if sz[1] == 163842
    dim = (27307, 1, 6, size(vol,4))
    vol = reshape(vol, dim)
  end

  io = open(fname, "w")

  nb = 0

  nb += write(io, hdr.sizeof_hdr)
  nb += write(io, collect(hdr.data_type))
  nb += write(io, collect(hdr.db_name))
  nb += write(io, hdr.extents)
  nb += write(io, hdr.session_error)
  nb += write(io, hdr.regular)
  nb += write(io, hdr.dim_info)
  nb += write(io, collect(hdr.dim))
  nb += write(io, hdr.intent_p1)
  nb += write(io, hdr.intent_p2)
  nb += write(io, hdr.intent_p3)
  nb += write(io, hdr.intent_code)
  nb += write(io, hdr.datatype)
  nb += write(io, hdr.bitpix)
  nb += write(io, hdr.slice_start)
  nb += write(io, collect(hdr.pixdim))
  nb += write(io, hdr.vox_offset)
  nb += write(io, hdr.scl_slope)
  nb += write(io, hdr.scl_inter)
  nb += write(io, hdr.slice_end)
  nb += write(io, hdr.slice_code)
  nb += write(io, hdr.xyzt_units)
  nb += write(io, hdr.cal_max)
  nb += write(io, hdr.cal_min)
  nb += write(io, hdr.slice_duration)
  nb += write(io, hdr.toffset)
  nb += write(io, hdr.glmax)
  nb += write(io, hdr.glmin)
  nb += write(io, collect(hdr.descrip))
  nb += write(io, collect(hdr.aux_file))
  nb += write(io, hdr.qform_code)
  nb += write(io, hdr.sform_code)
  nb += write(io, hdr.quatern_b)
  nb += write(io, hdr.quatern_c)
  nb += write(io, hdr.quatern_d)
  nb += write(io, hdr.quatern_x)
  nb += write(io, hdr.quatern_y)
  nb += write(io, hdr.quatern_z)
  nb += write(io, collect(hdr.srow_x))
  nb += write(io, collect(hdr.srow_y))
  nb += write(io, collect(hdr.srow_z))
  nb += write(io, collect(hdr.intent_name))
  nb += write(io, collect(hdr.magic))

  # Pad to get to 352 bytes (header size is 348)
  nb += write(io, zeros(UInt8,4))

  if hdr.datatype == 2
    dtype = UInt8
  elseif hdr.datatype == 4
    dtype = Int16
  elseif hdr.datatype == 8
    dtype = Int32
  elseif hdr.datatype == 16
    dtype = Float32
  elseif hdr.datatype == 64
    dtype = Float64
  elseif hdr.datatype == 256
    dtype = Int8
  elseif hdr.datatype == 512
    dtype = UInt16
  elseif hdr.datatype == 768
    dtype = UInt32
  else
    println("WARNING: data type " * string(hdr.datatype) *
            " not supported, but writing as float")
    dtype = Float32
  end

  nb += write(io, dtype.(vol))

  close(io)

  err = (nb != (sizeof(Float32) * 36 +
                sizeof(Int32) * 4 +
                sizeof(Int16) * 16 +
                sizeof(Int8) * 2 +
                sizeof(UInt8) * 158 +
                sizeof(dtype) * length(vol)))

  if gzip_needed
    cmd = `gzip -f $fname`
    run(cmd)
  end

  return err
end

"""
    mri_read_bfiles(infile1::String, infile2::String)

Read a DWI b-value table and gradient table from text files `infile1` and
`infile2`. The two input files can be specified in any order. The gradient
table file must contain 3 times as many entries as the b-value table file.

Return the b-value table as a vector of size n and the gradient table as a
matrix of size (n, 3).
"""
function mri_read_bfiles(infile1::String, infile2::String)

  tab = []

  for infile in (infile1, infile2)
    if ~isfile(infile)
      error("Could not open " * infile)
    end

    push!(tab, readdlm(infile))

    if ~all(isa.(tab[end], Number))
      error("File " * infile * " contains non-numeric entries")
    end

    if size(tab[end], 2) > size(tab[end], 1)
      tab[end] = permutedims(tab[end], 2:1)
    end

    if size(tab[end], 2) == 1
      tab[end] = tab[end][:,1]
    end
  end

  if size(tab[1], 1) != size(tab[2], 1)
    error("Dimension mismatch between tables in " *
          infile1 * " " * string(size(tab[1])) * " and " *
          infile2 * " " * string(size(tab[2])))
  end

  if (size(tab[1], 2) == 1 && size(tab[2], 2) == 3) ||
     (size(tab[1], 2) == 3 && size(tab[2], 2) == 1)
    return tab[1], tab[2]
  else
    error("Wrong number of entries in tables " *
          infile1 * " " * string(size(tab[1])) * " and " *
          infile2 * " " * string(size(tab[2])))
  end
end

"""
    mri_read_bfiles!(dwi::MRI, infile1::String, infile2::String)

Set the .bval and .bvec fields of the MRI structure `dwi`, by reading a DWI
b-value table and gradient table from the text files `infile1` and `infile2`.
The two input files can be specified in any order. The gradient table file
must contain 3 times as many entries as the b-value table file.

Return the b-value table as a vector of size n and the gradient table as a
matrix of size (n, 3).
"""
function mri_read_bfiles!(dwi::MRI, infile1::String, infile2::String)

  (tab1, tab2) = mri_read_bfiles(infile1, infile2)

  if size(tab1, 1) != size(dwi.vol, 4)
    error("Number of frames in volume (" * string(size(dwi.vol, 4)) *
          ") does not match dimensions of table in " *
          infile1 * " " * string(size(tab1)))
  end

  if size(tab1, 2) == 1
    dwi.bval = tab1
    dwi.bvec = tab2
  else
    dwi.bval = tab2
    dwi.bvec = tab1
  end

  # Normalize gradient vectors
  dwi.bvec = dwi.bvec ./ sqrt.(sum(dwi.bvec.^2, dims=2))
  dwi.bvec[isnan.(dwi.bvec)] .= Float32(0)

  return tab1, tab2
end

"""
    mri_read(inbase::String, type::DataType, headeronly::Bool=false, permuteflag::Bool=false)

Read a set of image files with base name `inbase`, which were generated by
an analysis, into a struct of a custom type (e.g., `DTI`)
"""
function mri_read(inbase::String, type::DataType, headeronly::Bool=false, permuteflag::Bool=false)

  if !isstructtype(type)
    error("Type " * type * " is not a structure")
  end

  absbase = abspath(inbase)
  flist = readdir(dirname(absbase), join=true)

  inputs = []

  for var in fieldnames(type)
    ftype = fieldtype.(type, var)

    if ftype == MRI
      infile = absbase * "_" * string(var) * ".nii.gz"
      push!(inputs, mri_read(infile, headeronly, permuteflag))
    elseif ftype == Vector{MRI}
      inpat = Regex("^" * absbase * "_" * string(var) * "[0-9]*.nii.gz\$")
      infiles = flist[.!isnothing.(match.(inpat, flist))]
      push!(inputs, mri_read.(infiles, headeronly, permuteflag))
    else
      infile = absbase * "_" * string(var) * ".txt"
      mat = Float32.(readdlm(infile))
      if length(mat) == 1
        push!(inputs, mat[1])
      else
        push!(inputs, mat)
      end
    end
  end

  return type(inputs...)
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
function DofInterpolate(basis::LagrangeBasis, field)

    num_bfs = numfunctions(basis)

    res = Vector(undef, num_bfs)

    for  b in 1 : num_bfs
        bfs = basis.fns[b]

        shape = bfs[1]

        cellid = shape.cellid
        refid = shape.refid

        cell = cells(basis.geo)[cellid]

        tria = chart(basis.geo, cell)

        vid =  cell[refid]

        if refid == 1
            v = neighborhood(tria, [1, 0])
        elseif refid == 2
            v = neighborhood(tria, [0, 1])
        else
            v = neighborhood(tria, [0, 0])
        end

        res[b] = field(v)
    end

    return res
end

function DofInterpolate(basis::RTBasis, field)

    num_bfs = numfunctions(basis)

    res = Vector(undef, num_bfs)

    for b in 1 : num_bfs

        bfs = basis.fns[b]

        shape = bfs[2]

        cellid = shape.cellid
        refid = shape.refid
        coeff = shape.coeff

        cell = cells(basis.geo)[cellid]

        e = refid

        v1 = cell[mod1(e+1,3)]
        v2 = cell[mod1(e+2,3)]

        edge = simplex(basis.geo.vertices[[v1,v2]]...)
        t = tangents(center(edge),1)
        tria = chart(basis.geo, cell)

        n = normal(center(tria))

        bn = normalize(cross(n,t))

        v = center(edge)

        res[b] = volume(edge)*coeff*dot(field(v),bn)

    end

    return res

end

function DofInterpolate(basis::BDMBasis, field)

    num_bfs = numfunctions(basis)

    res = Vector(undef, num_bfs)

    for b in 1 : num_bfs

        bfs = basis.fns[b]

        shape = bfs[2]

        cellid = shape.cellid
        refid = shape.refid
        coeff = shape.coeff

        cell = cells(basis.geo)[cellid]

        e = (refid-1)÷2+1

        v1 = cell[mod1(e+1,3)]
        v2 = cell[mod1(e+2,3)]

        edge = simplex(basis.geo.vertices[[v1,v2]]...)
        t = tangents(center(edge),1)
        tria = chart(basis.geo, cell)

        n = normal(center(tria))

        bn = normalize(cross(n,t))

        if refid == 1
            v = neighborhood(tria, [0, 1])
        elseif refid == 2
            v = neighborhood(tria, [0, 0])
        elseif refid == 3
            v = neighborhood(tria, [0, 0])
        elseif refid == 4
            v = neighborhood(tria, [1, 0])
        elseif refid == 5
            v = neighborhood(tria, [1, 0])
        else
            v = neighborhood(tria, [0, 1])
        end

        res[b] = volume(edge)*coeff*dot(field(v),bn)

    end

    return res

end

function DofInterpolate(basis::NDLCDBasis, field)

    num_bfs = numfunctions(basis)

    res = Vector(undef, num_bfs)

    for b in 1 : num_bfs

        bfs = basis.fns[b]

        shape = bfs[1]

        cellid = shape.cellid
        refid = shape.refid
        coeff = shape.coeff

        cell = cells(basis.geo)[cellid]

        f = refid

        v1 = cell[mod1(f+1,4)]
        v2 = cell[mod1(f+2,4)]
        v3 = cell[mod1(f+3,4)]

        v4 = cell[mod1(f,4)]

        face = simplex(basis.geo.vertices[[v1,v2,v3]]...)

        n = normal(center(face))

        inside = dot(n, basis.geo.vertices[v4]-basis.geo.vertices[v1])

        n *= -sign(inside)

        v = center(face)

        res[b] = volume(face)*coeff*dot(field(v),n)

    end

    return res

end

function DofInterpolate(basis::BEAST.BDM3DBasis, field)

    num_bfs = numfunctions(basis)

    res = Vector(undef, num_bfs)

    for b in 1 : num_bfs

        bfs = basis.fns[b]

        shape = bfs[1]

        cellid = shape.cellid
        refid = shape.refid
        coeff = shape.coeff

        cell = cells(basis.geo)[cellid]

        f = (refid-1)÷3+1

        v1 = cell[mod1(f+1,4)]
        v2 = cell[mod1(f+2,4)]
        v3 = cell[mod1(f+3,4)]

        v4 = cell[mod1(f,4)]

        face = simplex(basis.geo.vertices[[v1,v2,v3]]...)

        n = normal(center(face))

        inside = dot(n, basis.geo.vertices[v4]-basis.geo.vertices[v1])

        n *= -sign(inside)

        tet = simplex(basis.geo.vertices[cell]...)

        if refid == 1
            v = neighborhood(tet, [0, 1, 0])
        elseif refid == 2
            v = neighborhood(tet, [0, 0, 1])
        elseif refid == 3
            v = neighborhood(tet, [0, 0, 0])
        elseif refid == 4
            v = neighborhood(tet, [0, 0, 1])
        elseif refid == 5
            v = neighborhood(tet, [0, 0, 0])
        elseif refid == 6
            v = neighborhood(tet, [1, 0, 0])
        elseif refid == 7
            v = neighborhood(tet, [0, 0, 0])
        elseif refid == 8
            v = neighborhood(tet, [1, 0, 0])
        elseif refid == 9
            v = neighborhood(tet, [0, 1, 0])
        elseif refid == 10
            v = neighborhood(tet, [1, 0, 0])
        elseif refid == 11
            v = neighborhood(tet, [0, 1, 0])
        else
            v = neighborhood(tet, [0, 0, 1])
        end

        res[b] = volume(face)*coeff*dot(field(v),n)

    end

    return res

end

function Centervalues(mesh::Mesh,f::Functional)
    num_cells = numcells(mesh)

    res = Vector{SVector{3,Complex{Float64}}}(undef, num_cells)

    for i in 1:num_cells
        cell = cells(mesh)[i]

        tet = chart(mesh, cell)

        v = center(tet)

        res[i] = f(v)

    end

    return res

end

function EvalCenter(basis::Space, coeff)

    mesh = basis.geo
    num_cells = numcells(mesh)
    num_bfs = numfunctions(basis)
    ref_space = refspace(basis)

    res = Vector{SVector{3,Complex{Float64}}}(undef, num_cells)

    for i in 1:num_cells
        res[i] = [0,0,0]
    end

    for b in 1 : num_bfs

        bfs = basis.fns[b]

        for shape in bfs

            cellid = shape.cellid
            refid = shape.refid
            a = shape.coeff

            cell = cells(mesh)[cellid]

            tet = chart(mesh, cell)

            v = center(tet)

            local_bfs = ref_space(v)

            res[cellid] += a*coeff[b]*local_bfs[refid].value
        end

    end

    return res

end


text:
using Robotlib, Test, Random
@testset "Calib Force" begin
    N = 100
    Rf = Robotlib.orthonormal(randn(3, 3))
    mf = 1
    POSES = cat([Robotlib.orthonormal(randn(3, 3)) for i = 1:N]..., dims = 3)
    forceoffs = zeros(3)
    forces =
        hcat([
            Rf' * (POSES[1:3, 1:3, i]' * [0, 0, mf * -9.82] - forceoffs)
            for i = 1:N
        ]...)' |> copy

    R, m = calib_force(
        POSES,
        forces;
        offset = false,
        verbose = false,
    )

    @test R ≈ Rf
    @test m ≈ mf

    R = Robotlib.calib_force2(POSES, forces, mf)
    @test R ≈ Rf

    R,g,m = Robotlib.calib_force_iterative(POSES, forces, randn(3))
    @test R ≈ Rf
    @test m ≈ mf
    @test g ≈ [0,0,-9.82*mf]

    R,g = Robotlib.calib_force_eigen(POSES, forces)
    @test R ≈ Rf
    @test g ≈ [0,0,-9.82*mf]

    forceoffs = randn(3)
    forces =
        hcat([
            Rf' * (POSES[1:3, 1:3, i]' * [0, 0, mf * -9.82] - forceoffs)
            for i = 1:N
        ]...)'

    R, m, offs = calib_force(
        POSES,
        forces;
        offset = true,
        verbose = false,
    )

    @test R ≈ Rf
    @test m ≈ mf
    @test offs ≈ forceoffs

end

# using Plots, DSP, Robotlib, Random, LaTeXStrings
# pgfplots()
# σr = exp10.(range(-4,stop=1,length=1000))
# res = map(σr) do σ
#
#     N         = 100
#     Rf        = expω(1randn(3))
#     mf        = 1
#     POSES     = cat([expω(randn(3)) for i = 1:N]..., dims=3)
#     forceoffs = zeros(3)
#     forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*[0, 0, mf*-9.82] - forceoffs) for i = 1:N]...)'
#     forces  .+= σ*randn(size(forces))
#
#     R1,m = calib_force(POSES,forces,1.1mf; offset=false)
#
#     R2 = calib_force2(POSES,forces,1.1mf)
#     sum(abs2,R1-Rf), sum(abs2,R2-Rf)
# end
#
# e1,e2 = getindex.(res,1),getindex.(res,2)
# plot(σr,filtfilt(ones(20),[20],[e1 e2]), xscale=:log10, yscale=:log10, lab=["Relaxation" "Cayley"], xlabel=L"$\sigma$ [N]", ylabel=L"$||R_e||_F^2$",legend=:bottomright)
# hline!([3, 3], lab="", l=:dash)
# savefig2("/home/fredrikb/phdthesis/calibpaper/figs/calib_force2.tex")

max_forks_repo_name: baggepinnen/Robotlib


text:
# --------------- #
# Display methods #
# --------------- #

function show_extra(me::MomentBasedEstimator{T}) where T<:GMMEstimator
    j, p = J_test(me)
    "\nJ-test: $(round(j; digits = 3)) (P-value: $(round(p, digits = 3)))\n"
end

show_extra(me::MomentBasedEstimator{T}) where T<:MDEstimator = "\n"

function rshow(me::MomentBasedEstimator{T}) where T<:GMMEstimator
    println("GMM")
end

function rshow(me::MomentBasedEstimator{T}) where T<:MDEstimator
    println("MD")
end

function Base.show(io::IO, ::MIME"text/plain", me::T) where T<:MomentBasedEstimator
  s = symbolfy(me)
  println(io, "\n$s [k = $(npar(me)), m = $(nmom(me))]\n")

  if me.status[1] == :Solved
    # Only gives information if the solver
    # converged
    if (status(me) == :Optimal)
      # get coef table and j-test
      ct = coeftable(me)
      # print coefficient table
      println(io, "Coefficients:\n")
      show(io, ct)
      # Then show extra information for this type
      println(io, show_extra(me))
    else
      println("The optimization did not converge to a local solution")
    end
  end
end

function symbolfy(me::T) where T<:MomentBasedEstimator
    s = split(string(T), "MomentBasedEstimator")[2]
    v = split(s, "{")[2]
    Symbol(v)

    # if (v == "GMMEstimator{")
    #     Symbol(v, split(s[4], "{")[1], ", ", split(s[10], ",")[1], ", ", split(s[11], ",")[1], "}")
    # elseif (s[3] == "MDEstimator{")
    #     if (s[4] == "AnaGradMomFun{" || s[4] == "AnaFullMomFun{")
    #       Symbol(s[3], split(s[4], "{")[1], ", ", split(s[7], "},Divergences.")[2], " ", split(s[8], ",")[1], "}")
    #     else
    #       Symbol(s[3], split(s[4], "{")[1], ", ", split(s[9], "},Divergences.")[2], " ", split(s[10], ",")[1], "}")
    #     end
    # end
end


text:
using RandomizedLinAlg
using Test, LinearAlgebra, Random

@testset "rsvd_fnkz" begin

    Random.seed!(1234321)

    @testset "Rank 2 matrix" begin
        A = reshape(1 : 64, 8, 8)
        S = rsvd_fnkz(A, 2, ϵ=1e-9)
        @test norm(A - S.U * Diagonal(S.S) * S.Vt) ≤ 1e-9
    end

    @testset "Linearly dependent matrix" begin
        A = [collect(1 : 8) collect(1 : 8) collect(1 : 8) collect(1 : 8)]
        S = rsvd_fnkz(A, 4)
        @test norm(A - S.U * Diagonal(S.S) * S.Vt) ≤ 1e-9
    end

    #SVD of a matrix of zeros should be empty
    @testset "Zero matrix" begin
        S = rsvd_fnkz(zeros(10, 10), 1)
        @test S.U == zeros(10, 0)
        @test S.S == []
        @test S.Vt == zeros(0, 10)
    end

end