text:
module ParXRec

using AxisArrays
using Dates
using LightXML
using Printf

export readrec, readpar, readxml, echotimes, voxelsize, recdims

###
### Types
###

const UNKNOWN_KEY = :Other

"""
    mutable struct SeriesInfo
        Patient_Name            :: String            = ""
        Examination_Name        :: String            = ""
        Protocol_Name           :: String            = ""
        Examination_Date        :: Date              = Date(1970, 1, 1)
        Examination_Time        :: Time              = Time(0, 0, 0)
        Series_Data_Type        :: String            = ""
        Aquisition_Number       :: Int               = 0
        Reconstruction_Number   :: Int               = 0
        Scan_Duration           :: Float64           = 0
        Max_No_Phases           :: Int               = 1
        Max_No_Echoes           :: Int               = 1
        Max_No_Slices           :: Int               = 1
        Max_No_Dynamics         :: Int               = 1
        Max_No_Mixes            :: Int               = 1
        Max_No_B_Values         :: Int               = 1
        Max_No_Gradient_Orients :: Int               = 1
        No_Label_Types          :: Int               = 0
        Patient_Position        :: String            = ""
        Preparation_Direction   :: String            = ""
        Technique               :: String            = ""
        Scan_Resolution_X       :: Int               = 0
        Scan_Resolution_Y       :: Int               = 0
        Scan_Mode               :: String            = ""
        Repetition_Times        :: NTuple{2,Float64} = (0, 0)
        FOV_AP                  :: Float64           = 0
        FOV_FH                  :: Float64           = 0
        FOV_RL                  :: Float64           = 0
        Water_Fat_Shift         :: Float64           = 0
        Angulation_AP           :: Float64           = 0
        Angulation_FH           :: Float64           = 0
        Angulation_RL           :: Float64           = 0
        Off_Center_AP           :: Float64           = 0
        Off_Center_FH           :: Float64           = 0
        Off_Center_RL           :: Float64           = 0
        Flow_Compensation       :: Bool              = false
        Presaturation           :: Bool              = false
        Phase_Encoding_Velocity :: NTuple{3,Float64} = (0, 0, 0)
        MTC                     :: Bool              = false
        SPIR                    :: Bool              = false
        EPI_factor              :: Int               = 0
        Dynamic_Scan            :: Bool              = false
        Diffusion               :: Bool              = false
        Diffusion_Echo_Time     :: Float64           = 0
        Other                   :: Dict{String,Any}  = Dict()
    end

Struct containing .PAR or .XML series info.
"""
Base.@kwdef mutable struct SeriesInfo
    Patient_Name            :: String            = ""
    Examination_Name        :: String            = ""
    Protocol_Name           :: String            = ""
    Examination_Date        :: Date              = Date(1970, 1, 1)
    Examination_Time        :: Time              = Time(0, 0, 0)
    Series_Data_Type        :: String            = ""
    Aquisition_Number       :: Int               = 0
    Reconstruction_Number   :: Int               = 0
    Scan_Duration           :: Float64           = 0
    Max_No_Phases           :: Int               = 1
    Max_No_Echoes           :: Int               = 1
    Max_No_Slices           :: Int               = 1
    Max_No_Dynamics         :: Int               = 1
    Max_No_Mixes            :: Int               = 1
    Max_No_B_Values         :: Int               = 1
    Max_No_Gradient_Orients :: Int               = 1
    No_Label_Types          :: Int               = 0
    Patient_Position        :: String            = ""
    Preparation_Direction   :: String            = ""
    Technique               :: String            = ""
    Scan_Resolution_X       :: Int               = 0
    Scan_Resolution_Y       :: Int               = 0
    Scan_Mode               :: String            = ""
    Repetition_Times        :: NTuple{2,Float64} = (0, 0)
    FOV_AP                  :: Float64           = 0
    FOV_FH                  :: Float64           = 0
    FOV_RL                  :: Float64           = 0
    Water_Fat_Shift         :: Float64           = 0
    Angulation_AP           :: Float64           = 0
    Angulation_FH           :: Float64           = 0
    Angulation_RL           :: Float64           = 0
    Off_Center_AP           :: Float64           = 0
    Off_Center_FH           :: Float64           = 0
    Off_Center_RL           :: Float64           = 0
    Flow_Compensation       :: Bool              = false
    Presaturation           :: Bool              = false
    Phase_Encoding_Velocity :: NTuple{3,Float64} = (0, 0, 0)
    MTC                     :: Bool              = false
    SPIR                    :: Bool              = false
    EPI_factor              :: Int               = 0
    Dynamic_Scan            :: Bool              = false
    Diffusion               :: Bool              = false
    Diffusion_Echo_Time     :: Float64           = 0
    Other                   :: Dict{String,Any}  = Dict()
end

"""
    mutable struct ImageKey
        Slice       :: Int    = 1
        Echo        :: Int    = 1
        Dynamic     :: Int    = 1
        Phase       :: Int    = 1
        BValue      :: Int    = 1
        Grad_Orient :: Int    = 1
        Label_Type  :: String = "" # Enumeration
        Type        :: String = "" # Enumeration
        Sequence    :: String = "" # Enumeration
        Index       :: Int    = 0
    end

Struct containing .PAR or .XML image keys.
"""
Base.@kwdef mutable struct ImageKey
    Slice       :: Int    = 1
    Echo        :: Int    = 1
    Dynamic     :: Int    = 1
    Phase       :: Int    = 1
    BValue      :: Int    = 1
    Grad_Orient :: Int    = 1
    Label_Type  :: String = "" # Enumeration
    Type        :: String = "" # Enumeration
    Sequence    :: String = "" # Enumeration
    Index       :: Int    = 0
end

"""
    mutable struct ImageInfo
        Key                        :: ImageKey
        Pixel_Size                 :: Int               = 0
        Scan_Percentage            :: Float64           = 0
        Resolution_X               :: Int               = 0
        Resolution_Y               :: Int               = 0
        Rescale_Intercept          :: Float64           = 0
        Rescale_Slope              :: Float64           = 0
        Scale_Slope                :: Float64           = 0
        Window_Center              :: Float64           = 0
        Window_Width               :: Float64           = 0
        Slice_Thickness            :: Float64           = 0
        Slice_Gap                  :: Float64           = 0
        Display_Orientation        :: String            = ""      # Enumeration
        fMRI_Status_Indication     :: Int               = 0
        Image_Type_Ed_Es           :: String            = ""      # Enumeration
        Pixel_Spacing              :: NTuple{2,Float64} = (0, 0)
        Echo_Time                  :: Float64           = 0
        Dyn_Scan_Begin_Time        :: Float64           = 0
        Trigger_Time               :: Float64           = 0
        Diffusion_B_Factor         :: Float64           = 0
        No_Averages                :: Float64           = 0
        Image_Flip_Angle           :: Float64           = 0
        Cardiac_Frequency          :: Int               = 0
        Min_RR_Interval            :: Int               = 0
        Max_RR_Interval            :: Int               = 0
        TURBO_Factor               :: Int               = 0
        Inversion_Delay            :: Float64           = 0
        Contrast_Type              :: String            = ""
        Diffusion_Anisotropy_Type  :: String            = ""
        Diffusion_AP               :: Float64           = 0
        Diffusion_FH               :: Float64           = 0
        Diffusion_RL               :: Float64           = 0
        Angulation_AP              :: Float64           = 0
        Angulation_FH              :: Float64           = 0
        Angulation_RL              :: Float64           = 0
        Offcenter_AP               :: Float64           = 0
        Offcenter_FH               :: Float64           = 0
        Offcenter_RL               :: Float64           = 0
        Slice_Orientation          :: String            = ""      # Enumeration
        Image_Planar_Configuration :: Int               = 0
        Samples_Per_Pixel          :: Int               = 0
        Other                      :: Dict{String,Any}  = Dict()
    end

Struct containing .PAR or .XML image info.
"""
Base.@kwdef mutable struct ImageInfo
    Key                        :: ImageKey
    Pixel_Size                 :: Int               = 0
    Scan_Percentage            :: Float64           = 0
    Resolution_X               :: Int               = 0
    Resolution_Y               :: Int               = 0
    Rescale_Intercept          :: Float64           = 0
    Rescale_Slope              :: Float64           = 0
    Scale_Slope                :: Float64           = 0
    Window_Center              :: Float64           = 0
    Window_Width               :: Float64           = 0
    Slice_Thickness            :: Float64           = 0
    Slice_Gap                  :: Float64           = 0
    Display_Orientation        :: String            = ""      # Enumeration
    fMRI_Status_Indication     :: Int               = 0
    Image_Type_Ed_Es           :: String            = ""      # Enumeration
    Pixel_Spacing              :: NTuple{2,Float64} = (0, 0)
    Echo_Time                  :: Float64           = 0
    Dyn_Scan_Begin_Time        :: Float64           = 0
    Trigger_Time               :: Float64           = 0
    Diffusion_B_Factor         :: Float64           = 0
    No_Averages                :: Float64           = 0
    Image_Flip_Angle           :: Float64           = 0
    Cardiac_Frequency          :: Int               = 0
    Min_RR_Interval            :: Int               = 0
    Max_RR_Interval            :: Int               = 0
    TURBO_Factor               :: Int               = 0
    Inversion_Delay            :: Float64           = 0
    Contrast_Type              :: String            = ""
    Diffusion_Anisotropy_Type  :: String            = ""
    Diffusion_AP               :: Float64           = 0
    Diffusion_FH               :: Float64           = 0
    Diffusion_RL               :: Float64           = 0
    Angulation_AP              :: Float64           = 0
    Angulation_FH              :: Float64           = 0
    Angulation_RL              :: Float64           = 0
    Offcenter_AP               :: Float64           = 0
    Offcenter_FH               :: Float64           = 0
    Offcenter_RL               :: Float64           = 0
    Slice_Orientation          :: String            = ""      # Enumeration
    Image_Planar_Configuration :: Int               = 0
    Samples_Per_Pixel          :: Int               = 0
    Other                      :: Dict{String,Any}  = Dict()
end

"""
    struct RecHeader
        series :: SeriesInfo
        images :: Vector{ImageInfo}
    end

Struct containing .PAR or .XML header.
"""
struct RecHeader
    series :: SeriesInfo
    images :: Vector{ImageInfo}
end

"""
    struct Rec{A<:Union{Nothing, AbstractArray}}
        hdr  :: RecHeader
        data :: A
    end

Struct containing .PAR or .XML header and .REC data array (or nothing).
"""
struct Rec{A<:Union{Nothing, AbstractArray}}
    hdr  :: RecHeader
    data :: A
end

###
### Interface
###

"""
    load(path; kwargs...) = load(Float32, path; kwargs...)

    load(
        ::Type{T},
        path::AbstractString;
        load_data::Bool = true,
        scale_data::Bool = true
    ) where {T<:AbstractFloat} -> Rec{Union{Nothing, AxisArray}}

Load Philips PAR/REC or XML/REC header and data array.

### Arguments
- `path::AbstractString`: path to .REC or .XML or .PAR file

### Keywords
- `load_data::Bool = true`: load data and header (true) or header only (false)
- `scale_data::Bool = true`: scale data array with
    `(data + rescale slope * rescale intercept) / scale slope` (true) or
    return raw array (false)

### Returns
- `Rec{Union{Nothing, AxisArray}}`: struct containing header and data array
    if `load_data = true` or nothing if `load_data = false`
"""
load(path; kwargs...) = load(Float32, path; kwargs...)

function load(
    ::Type{T},
    path::AbstractString;
    load_data::Bool = true,
    scale_data::Bool = true,
) where {T<:AbstractFloat}
    hdrpath, recpath = parxrec_paths(path, load_data)
    _, ext = splitext(hdrpath)

    hdr = lowercase(ext) == ".par" ?
        readpar(hdrpath) :
        readxml(hdrpath)

    data = load_data ?
        readrec(T, recpath, hdr, scale_data=scale_data) :
        nothing

    return Rec(hdr, data)
end

###
### .REC
###

const _RECPREC = Dict{Int, DataType}(
    8  => UInt8,
    16 => UInt16,
    32 => Float32,
    64 => Float64,
)

"""
    readrec(path, hdr; kwargs...) = readrec(Float32, path, hdr; kwargs...)

    function readrec(
        ::Type{T},
        path::AbstractString,
        hdr::RecHeader;
        scale_data::Bool = true
    ) where {T<:AbstractFloat} -> AxisArray

Load .REC file.

### Arguments
- `path::AbstractString`: path to .REC file
- `hdr::RecHeader`: struct containing .PAR or .XML header

### Keywords
- `scale_data::Bool = true`: scale data array with
    `(data + rescale slope * rescale intercept) / scale slope` (true) or
    return raw array (false)

### Returns
- `AxisArray`: .REC data array
"""
readrec(path, hdr; kwargs...) = readrec(Float32, path, hdr; kwargs...)

function readrec(
    ::Type{T},
    path::AbstractString,
    hdr::RecHeader;
    scale_data::Bool = true
) where {T<:AbstractFloat}
    validate(hdr)

    imgs = sort(hdr.images, by = (i -> i.Key.Index))

    ps = Int(imgs[1].Pixel_Size)
    TR = _RECPREC[ps]

    dims = recdims(hdr)

    ax0 = (
        nx       = maximum(dims[:nx]),
        ny       = maximum(dims[:ny]),
        nz       = maximum(dims[:nz]),
        echoes   = maximum(dims[:echoes]),
        dynamics = maximum(dims[:dynamics]),
        phases   = maximum(dims[:phases]),
        bvals    = maximum(dims[:bvals]),
        grads    = maximum(dims[:grads]),
        labels   = length(dims[:labels]),
        types    = length(dims[:types]),
        seqs     = length(dims[:seqs]),
    )

    # remove singleton dimensions but keep nx, ny, nz
    ax = (; [(a, n) for (a, n) in pairs(ax0) if n > 1 || a ∈ (:nx, :ny, :nz)]...)

    nx = ax[:nx]
    ny = ax[:ny]

    labels = Dict(k => i for (i, k) in enumerate(dims[:labels]))
    types = Dict(k => i for (i, k) in enumerate(dims[:types]))
    seqs = Dict(k => i for (i, k) in enumerate(dims[:seqs]))

    L = LinearIndices(values(ax0))
    data = zeros(T, values(ax))
    slice = zeros(TR, nx * ny)
    slicet = zeros(TR, nx, ny)

    nb = nx * ny * ps÷8

    rsz = filesize(path)
    w = (nb * length(imgs)) - rsz
    w < 0 && @warn "REC file larger than expected from hdr"
    w > 0 && @warn "REC file smaller than expected from hdr"
    rsz -= nb

    open(path, "r") do io
        @inbounds for i in imgs
            pos = i.Key.Index * nb
            pos > rsz && continue

            seek(io, pos)
            read!(io, slicet)
            copyto!(slice, transpose(slicet))

            i3  = i.Key.Slice
            i4  = i.Key.Echo
            i5  = i.Key.Dynamic
            i6  = i.Key.Phase
            i7  = i.Key.BValue
            i8  = i.Key.Grad_Orient
            i9  = labels[i.Key.Label_Type]
            i10 = types[(i.Key.Type, i.Rescale_Intercept)]
            i11 = seqs[i.Key.Sequence]

            I = L[1,1,i3,i4,i5,i6,i7,i8,i9,i10,i11] : L[nx,ny,i3,i4,i5,i6,i7,i8,i9,i10,i11]

            if scale_data
                ss = i.Scale_Slope ≈ 0.0 ? 1.0 : inv(i.Scale_Slope)
                rs = i.Rescale_Slope ≈ 0.0 ? 1.0 : inv(i.Rescale_Slope)
                ri = rs * i.Rescale_Intercept
                data[I] .= T(ss) .* (slice .+ T(ri))
            else
                data[I] .= slice
            end
        end
    end

    return AxisArray(data, keys(ax)...)
end

###
### .XML
###

const _XML_HEADER           = "PRIDE_V5"
const _XML_SERIES_HEADER    = "Series_Info"
const _XML_IMAGE_ARR_HEADER = "Image_Array"
const _XML_IMAGE_HEADER     = "Image_Info"
const _XML_IMAGE_KEY_HEADER = "Key"
const _XML_ATTRIB_HEADER    = "Attribute"

"""
    readxml(path::AbstractString) -> RecHeader

Load .XML header.

### Arguments
- `path::AbstractString`: path to .XML file

### Returns
- `RecHeader`: struct containing .XML header
"""
function readxml(path::AbstractString)
    series = SeriesInfo()
    images = Vector{ImageInfo}()

    Ts = NamedTuple{fieldnames(SeriesInfo)}(SeriesInfo.types)
    Tk = NamedTuple{fieldnames(ImageKey)}(ImageKey.types)
    Ti = NamedTuple{fieldnames(ImageInfo)}(ImageInfo.types)

    xdoc = parse_file(path)
    xroot = root(xdoc)
    name(xroot) != _XML_HEADER &&
        error("Unknown XML/REC header: $(name(xroot)) (expected <$(_XML_HEADER)>)")

    # Series Info
    cur = find_element(xroot, _XML_SERIES_HEADER)
    cur === nothing &&
        error("<$(_XML_SERIES_HEADER)> tag not found")

    for e in child_elements(cur)
        name(e) != _XML_ATTRIB_HEADER && continue
        sethdr!(series, Ts, e)
    end

    # Image Info
    cur = find_element(xroot, _XML_IMAGE_ARR_HEADER)
    cur === nothing &&
        error("<$(_XML_IMAGE_ARR_HEADER)> tag not found")

    for c in child_elements(cur)
        name(c) != _XML_IMAGE_HEADER && continue

        ks = find_element(c, _XML_IMAGE_KEY_HEADER)
        ks === nothing &&
            error("<$(_XML_IMAGE_KEY_HEADER)> tag not found")

        imagekey = ImageKey()
        for e in child_elements(ks)
            name(e) != _XML_ATTRIB_HEADER && continue
            sethdr!(imagekey, Tk, e)
        end

        image = ImageInfo(Key = imagekey)
        for e in child_elements(c)
            name(e) != _XML_ATTRIB_HEADER && continue
            sethdr!(image, Ti, e)
        end

        push!(images, image)
    end

    free(xdoc)
    return RecHeader(series, images)
end

function sethdr!(h, T::NamedTuple, e::XMLElement)
    name = attribute(e, "Name")
    name === nothing && return nothing

    k = Symbol(replace(name, " " => "_"))
    v = content(e)

    if !haskey(T, k)
        _u = getfield(h, UNKNOWN_KEY)
        _u[name] = v
    elseif v != ""
        setfield!(h, k, _parse(T[k], v))
    end

    return h
end

###
### .PAR
###

const _PARSERIES = Dict{String, Union{Symbol, Tuple}}(
    "Patient name"                        => :Patient_Name,
    "Examination name"                    => :Examination_Name,
    "Protocol name"                       => :Protocol_Name,
    "Examination date/time"               => (:Examination_Date, :Examination_Time),
    "Series Type"                         => :Series_Data_Type,
    "Acquisition nr"                      => :Aquisition_Number,
    "Reconstruction nr"                   => :Reconstruction_Number,
    "Scan Duration [sec]"                 => :Scan_Duration,
    "Max. number of cardiac phases"       => :Max_No_Phases,
    "Max. number of echoes"               => :Max_No_Echoes,
    "Max. number of slices/locations"     => :Max_No_Slices,
    "Max. number of dynamics"             => :Max_No_Dynamics,
    "Max. number of mixes"                => :Max_No_Mixes,
    "Patient position"                    => :Patient_Position,
    "Preparation direction"               => :Preparation_Direction,
    "Technique"                           => :Technique,
    "Scan resolution  (x, y)"             => (:Scan_Resolution_X, :Scan_Resolution_Y),
    "Scan mode"                           => :Scan_Mode,
    "Repetition time [ms]"                => :Repetition_Times,
    "FOV (ap,fh,rl) [mm]"                 => (:FOV_AP, :FOV_FH, :FOV_RL),
    "Water Fat shift [pixels]"            => :Water_Fat_Shift,
    "Angulation midslice(ap,fh,rl)[degr]" => (:Angulation_AP, :Angulation_FH, :Angulation_RL),
    "Off Centre midslice(ap,fh,rl) [mm]"  => (:Off_Center_AP, :Off_Center_FH, :Off_Center_RL),
    "Flow compensation <0=no 1=yes> ?"    => :Flow_Compensation,
    "Presaturation     <0=no 1=yes> ?"    => :Presaturation,
    "Phase encoding velocity [cm/sec]"    => :Phase_Encoding_Velocity,
    "MTC               <0=no 1=yes> ?"    => :MTC,
    "SPIR              <0=no 1=yes> ?"    => :SPIR,
    "EPI factor        <0,1=no EPI>"      => :EPI_factor,
    "Dynamic scan      <0=no 1=yes> ?"    => :Dynamic_Scan,
    "Diffusion         <0=no 1=yes> ?"    => :Diffusion,
    "Diffusion echo time [ms]"            => :Diffusion_Echo_Time,
    "Max. number of diffusion values"     => :Max_No_B_Values,
    "Max. number of gradient orients"     => :Max_No_Gradient_Orients,
    "Number of label types   <0=no ASL>"  => :No_Label_Types,
)

const _PARSERIES3 = Dict{String, Union{Symbol, Tuple}}(
    "Image pixel size [8 or 16 bits]"     => :Pixel_Size,
    "Scan percentage"                     => :Scan_Percentage,
    "Recon resolution (x, y)"             => (:Resolution_X, :Resolution_Y),
    "Slice thickness [mm]"                => :Slice_Thickness,
    "Slice gap [mm]"                      => :Slice_Gap,
    "Number of averages"                  => :No_Averages,
    "Cardiac frequency"                   => :Cardiac_Frequency,
    "Min. RR interval"                    => :Min_RR_Interval,
    "Max. RR interval"                    => :Max_RR_Interval,
    "TURBO factor      <0=no turbo>"      => :TURBO_Factor,
    "Inversion delay [ms]"                => :Inversion_Delay,
)

const _PARKEY = Dict{Int, Symbol}(
    1  => :Slice,
    2  => :Echo,
    3  => :Dynamic,
    4  => :Phase,
    5  => :Type,
    6  => :Sequence,
    7  => :Index,
    42 => :BValue,
    43 => :Grad_Orient,
    49 => :Label_Type,
)

const _PARIMAGE = Dict{Int, Symbol}(
    8  => :Pixel_Size,
    9  => :Scan_Percentage,
    10 => :Resolution_X,
    11 => :Resolution_Y,
    12 => :Rescale_Intercept,
    13 => :Rescale_Slope,
    14 => :Scale_Slope,
    15 => :Window_Center,
    16 => :Window_Width,
    17 => :Angulation_AP,
    18 => :Angulation_FH,
    19 => :Angulation_RL,
    20 => :Offcenter_AP,
    21 => :Offcenter_FH,
    22 => :Offcenter_RL,
    23 => :Slice_Thickness,
    24 => :Slice_Gap,
    25 => :Display_Orientation,
    26 => :Slice_Orientation,
    27 => :fMRI_Status_Indication,
    28 => :Image_Type_Ed_Es,
    29 => :Pixel_Spacing,  # [29, 30]
    31 => :Echo_Time,
    32 => :Dyn_Scan_Begin_Time,
    33 => :Trigger_Time,
    34 => :Diffusion_B_Factor,
    35 => :No_Averages,
    36 => :Image_Flip_Angle,
    37 => :Cardiac_Frequency,
    38 => :Min_RR_Interval,
    39 => :Max_RR_Interval,
    40 => :TURBO_Factor,
    41 => :Inversion_Delay,
    44 => :Contrast_Type,
    45 => :Diffusion_Anisotropy_Type,
    46 => :Diffusion_AP,
    47 => :Diffusion_FH,
    48 => :Diffusion_RL,
)

const _PARIMAGE3 = Dict{Int, Symbol}(
    8  => :Rescale_Intercept,
    9  => :Rescale_Slope,
    10 => :Scale_Slope,
    11 => :Window_Center,
    12 => :Window_Width,
    13 => :Angulation_AP,
    14 => :Angulation_FH,
    15 => :Angulation_RL,
    16 => :Offcenter_AP,
    17 => :Offcenter_FH,
    18 => :Offcenter_RL,
    19 => :Display_Orientation,
    20 => :Slice_Orientation,
    21 => :fMRI_Status_Indication,
    22 => :Image_Type_Ed_Es,
    23 => :Pixel_Spacing,  # [23, 24]
    25 => :Echo_Time,
    26 => :Dyn_Scan_Begin_Time,
    27 => :Trigger_Time,
    28 => :Diffusion_B_Factor,
    29 => :Image_Flip_Angle,
)

"""
    readpar(path::AbstractString) -> RecHeader

Load .PAR header.

### Arguments
- `path::AbstractString`: path to .PAR file

### Returns
- `RecHeader`: struct containing .PAR header
"""
function readpar(path::AbstractString)
    version = ""

    series = SeriesInfo()
    images = Vector{ImageInfo}()

    Ts = NamedTuple{fieldnames(SeriesInfo)}(SeriesInfo.types)
    Tk = NamedTuple{fieldnames(ImageKey)}(ImageKey.types)
    Ti = NamedTuple{fieldnames(ImageInfo)}(ImageInfo.types)

    doc = eachline(path)
    iter = iterate(doc)

    # PAR version number
    while iter !== nothing
        line, state = iter
        first(line) == '.' && break

        if occursin("CLINICAL TRYOUT", line)
            m = match(r"(?<=V)[0-9]\.?[0-9]*$", strip(line))
            if m !== nothing
                version = m.match
            end
            iter = iterate(doc, state)
            break
        end

        iter = iterate(doc, state)
    end

    v = tryparse(Float64, version)
    if v === nothing
        error("PAR version not found")

    elseif v < 3
        doc = nothing
        iter = nothing
        error("PAR V$(version) not supported")

    elseif v < 4
        v3 = true
        image3 = ImageInfo(Key = ImageKey())
        parimage = _PARIMAGE3

    else
        v3 = false
        parimage = _PARIMAGE
    end

    # Series Info / General Information
    while iter !== nothing
        line, state = iter

        if line == ""
            iter = iterate(doc, state)
            break
        end

        if line[1] != '.'
            (line[1] == ' ' || isdigit(line[1])) && break
            iter = iterate(doc, state)
            continue
        end

        if line[end] == ':'
            iter = iterate(doc, state)
            continue
        end

        ks, vs = strip.(split(line[2:end], ": "))

        if v3 && occursin("[msec]", ks)
            ks = replace(ks, "[msec]" => "[ms]")
        end

        if haskey(_PARSERIES, ks)
            k = _PARSERIES[ks]

            if k == :Repetition_Times
                if length(split(vs)) == 1
                    vs *= " 0.0"
                end
            end

            if k isa Symbol
                setfield!(series, k, _parse(Ts[k], vs))

            elseif k isa Tuple
                if occursin("Examination", ks)
                    x = split(vs, "/")
                    # in case month and/or day is 00
                    x[1] = replace(x[1], ".00" => ".01")
                else
                    x = split(vs)
                end

                x = strip.(x)

                for (i, _k) in enumerate(k)
                    setfield!(series, _k, _parse(Ts[_k], x[i]))
                end

            else
                error("Not implemented")
            end

        elseif v3 && haskey(_PARSERIES3, ks)
            k = _PARSERIES3[ks]

            if k isa Symbol
                setfield!(image3, k, _parse(Ti[k], vs))

            elseif k isa Tuple
                x = strip.(split(vs))
                for (i, _k) in enumerate(k)
                    setfield!(image3, _k, _parse(Ti[_k], x[i]))
                end

            else
                error("Not implemented")
            end

        else
            _u = getfield(series, UNKNOWN_KEY)
            _u[ks] = vs
        end

        iter = iterate(doc, state)
    end

    # Image Info
    while iter !== nothing
        line, state = iter

        if line == "" || (line[1] != ' ' && !isdigit(line[1]))
            iter = iterate(doc, state)
            continue
        end

        xs = split(line)
        N = length(xs)

        imagekey = ImageKey()
        for (i, k) in _PARKEY
            i > N && continue
            setfield!(imagekey, k, _parse(Tk[k], xs[i]))
        end

        image = ImageInfo(Key = imagekey)
        for (i, k) in parimage
            i > N && continue
            k == :Pixel_Spacing ?
                setfield!(image, k, _parse(Ti[k], xs[i]*" "*xs[i+1])) :
                setfield!(image, k, _parse(Ti[k], xs[i]))
        end

        if v3
            for v in values(_PARSERIES3)
                if v isa Tuple
                    for _v in v
                        setfield!(image, _v, getfield(image3, _v))
                    end
                else
                    setfield!(image, v, getfield(image3, v))
                end
            end
        end

        push!(images, image)
        iter = iterate(doc, state)
    end

    return RecHeader(series, images)
end

###
### Utility
###

"""
    echotimes(R::Rec)
    echotimes(H::RecHeader)
    echotimes(I::AbstractArray{ImageInfo}) -> NTuple{N, Float64}

Get echo times (in seconds).
"""
echotimes(R::Rec) = echotimes(R.hdr)
echotimes(H::RecHeader) = echotimes(H.images)
echotimes(I::AbstractArray{ImageInfo}) =
    1e-3 .* (sort!(unique!((i -> i.Echo_Time).(I)))...,)

"""
    voxelsize(R::Rec)
    voxelsize(H::RecHeader)
    voxelsize(I::ImageInfo) -> NTuple{3, Float64}
    voxelsize(I::AbstractArray{ImageInfo}) -> Union{NTuple{3, Float64}, Array{NTuple{3, Float64}}}

Get voxel size.
"""
voxelsize(R::Rec) = voxelsize(R.hdr)
voxelsize(H::RecHeader) = voxelsize(H.images)
voxelsize(I::ImageInfo) = (I.Pixel_Spacing..., I.Slice_Thickness + I.Slice_Gap)
voxelsize(I::AbstractArray{ImageInfo}) =
    unique!(voxelsize.(I)) |> v -> length(v) == 1 ? first(v) : v

"""
    recdims(H::RecHeader) -> Dict{Symbol, Int}

Array dimensions for REC data array from .PAR/.XML header. Dict is unordered.

### Keys
- `:nx`
- `:ny`
- `:nz`
- `:echoes`
- `:dynamics`
- `:phases`
- `:bvals`
- `:grads`
- `:labels`
- `:types`
- `:seqs`
"""
function recdims(H::RecHeader)
    K = (i -> i.Key).(H.images)

    @inline fun(f, k; by=identity) = sort!(unique!(f.(k)), by=by)

    nx       = fun(i -> i.Resolution_X, H.images)
    ny       = fun(i -> i.Resolution_Y, H.images)
    nz       = fun(k -> k.Slice, K)

    echoes   = fun(k -> k.Echo, K)
    dynamics = fun(k -> k.Dynamic, K)
    phases   = fun(k -> k.Phase, K)
    bvals    = fun(k -> k.BValue, K)
    grads    = fun(k -> k.Grad_Orient, K)

    labels   = fun(k -> k.Label_Type, K)
    seqs     = fun(k -> k.Sequence, K)

    # derived images have type = -1 in .PAR.
    # differentiate based on rescale intercept.
    types    = fun(
        i -> (i.Key.Type, i.Rescale_Intercept),
        H.images,
        by = t -> (t[1], abs(t[2])) # phase -3142, mag 0 -> flip order
    )

    return Dict(
        :nx       => nx,
        :ny       => ny,
        :nz       => nz,

        :echoes   => echoes,
        :dynamics => dynamics,
        :phases   => phases,
        :bvals    => bvals,
        :grads    => grads,

        :labels   => labels,
        :types    => types,
        :seqs     => seqs,
    )
end

###
### Misc
###

function parxrec_paths(path::AbstractString, load_data::Bool)
    !isfile(path) && error("File `$path` does not exist")

    path1, ext = splitext(path)
    ext = lowercase(ext)

    if ext == ".xml" || ext == ".par"
        hdrpath = path
        recpath = load_data ? (
            isfile(path1 * ".REC") ? path1 * ".REC" :
            isfile(path1 * ".rec") ? path1 * ".rec" :
            error(".REC file $path1 does not exist")
        ) : nothing

    elseif ext == ".rec"
        recpath = path
        hdrpath =
            isfile(path1 * ".XML") ? path1 * ".XML" :
            isfile(path1 * ".PAR") ? path1 * ".PAR" :
            isfile(path1 * ".xml") ? path1 * ".xml" :
            isfile(path1 * ".par") ? path1 * ".par" :
            error(".PAR/.XML file $path1 does not exist")
    else
        error("File $path does not have valid PAR/XML/REC extension")
    end

    return hdrpath, recpath
end

function validate(hdr::RecHeader)
    series = hdr.series
    images = hdr.images
    keys   = (i -> i.Key).(images)

    smax = (s, k) -> maximum(x -> getfield(x, k), s)
    slength = (s, k) -> length(unique(x -> getfield(x, k), s))

    vals = [
        (series.Max_No_Slices, :Slice, "slices"),
        (series.Max_No_Echoes, :Echo, "echoes"),
        (series.Max_No_Dynamics, :Dynamic, "dynamics"),
        (series.Max_No_Phases, :Phase, "phases"),
        (series.Max_No_B_Values, :BValue, "bvalues"),
        (series.Max_No_Gradient_Orients, :Grad_Orient, "gradient orientations"),
    ]

    # Check whether data size in series info is consistent with image info
    ws1 = (f, s, i) ->
        "Maximum number of $f not matching: $s (series info) vs. $i (image info)"

    # restrict to mag, real, imag, phas
    ks = filter(k -> k.Type ∈ ("0", "1", "2", "3", "M", "R", "I", "P"), keys)

    if !isempty(ks)
        for (s, k, f) in vals
            i = smax(ks, k)
            i != s && @warn ws1(f, s, i)
        end

        nl = series.No_Label_Types == 0 ? 1 : series.No_Label_Types
        i = slength(ks, :Label_Type)
        i != nl && @warn ws1("label types", series.No_Label_Types, i)
    end

    # Non-exhaustive image info check
    ws2 = (f, m, l) ->
        "Inconsistent number of $f in image info: $m (maximum) vs. $l (length(unique))"

    # image keys
    for (_, k, f) in vals
        m = smax(keys, k)
        l = slength(keys, k)
        l != m && @warn ws2(f, m, l)
    end

    m = smax(keys, :Index)
    l = slength(keys, :Index)
    l != m + 1 && @warn ws2("rec indices", m, l)

    # image info
    slength(images, :Slice_Thickness) != 1 &&
        @warn "Multiple values for `Slice Thickness` found"

    slength(images, :Slice_Gap) != 1 &&
        @warn "Multiple values for `Slice Gap` found"

    slength(images, :Pixel_Spacing) != 1 &&
        @warn "Multiple values for `Pixel Spacing` found"

    slength(images, :Pixel_Size) != 1 &&
        error("Multiple values for `Pixel Size` found")

    slength(images, :Resolution_X) != 1 &&
        error("Multiple values for `Resolution X` found")

    slength(images, :Resolution_Y) != 1 &&
        error("Multiple values for `Resolution Y` found")

    return nothing
end

_parse(type, str) = parse(type, str)

_parse(::Type{String}, str::AbstractString) = String(str)
_parse(::Type{Time}, str::AbstractString) = Time(str, dateformat"H:M:S")
_parse(::Type{Date}, str::AbstractString) = Date(str, dateformat"y.m.d")

_parse(::Type{Bool}, str::AbstractString) =
    str[1] == '0' || uppercase(str[1]) == 'N' ? false :
    str[1] == '1' || uppercase(str[1]) == 'Y' ? true :
    parse(Bool, str)

_parse(::Type{T}, str::AbstractString) where {T<:Tuple} =
    Tuple(_parse.(T.types, split(str)))

function Base.show(io::IO, s::T) where {T<:Union{SeriesInfo, ImageInfo, ImageKey}}
    ks = fieldnames(T)

    println(io, string(T))

    for k in ks
        v = getfield(s, k)
        @printf(io, "    %-27s : ", string(k))
        println(io, v)
    end
end

function Base.show(io::IO, s::RecHeader)
    N = length(s.images)

    println(io, typeof(s))
    println(io)
    print(io, s.series)

    if N > 0
        print(io, "\nImage 1/$(N)\n")
        print(io, s.images[1])
    end

    if N > 1
        println(io)
        print(io, "    .\n")
        print(io, "    .\n")
        print(io, "    .\n")
        println(io)
        print(io, "Image $(N)/$(N)\n")
        print(io, s.images[N])
    end
end

end # module

ext: jl


text:
using FINUFFT
using Test

@testset "Error handling and dumb inputs" begin

    xj = zeros(10)
    cj = complex(zeros(10))
    iflag = 1
    tol = 1e-15
    ms = 10

    nufft1d1(xj, cj, iflag, tol, ms) # This should pass

    @info("Testing error handling")

    # Tolerance too small
    try
        nufft1d1(xj, cj, iflag, 1e-100, ms)
    catch e
        @test e.errno == FINUFFT.ERR_EPS_TOO_SMALL
    end

    # Allocate too much
    opts = finufft_default_opts()
    spread_kerevalmeth = 0
    upsampfac = maxintfloat(typeof(opts.upsampfac))   # hack to alloc a lot
    try
        nufft1d1(xj, cj, iflag, tol, ms, spread_kerevalmeth=spread_kerevalmeth, upsampfac=upsampfac)
    catch e
        @test e.errno == FINUFFT.ERR_MAXNALLOC
    end

    # Too small upsampfac
    spread_kerevalmeth = 0
    upsampfac = 0.9                     # note 0 is auto-choice
    try
        nufft1d1(xj, cj, iflag, tol, ms, spread_kerevalmeth=spread_kerevalmeth, upsampfac=upsampfac)
    catch e
        @test e.errno == FINUFFT.ERR_UPSAMPFAC_TOO_SMALL
    end

    # Wrong beta (ie upsampfac not known for Horner poly eval rule)
    upsampfac = 1.5
    try
        nufft1d1(xj, cj, iflag, tol, ms, upsampfac=upsampfac)
    catch e
        @test e.errno == FINUFFT.HORNER_WRONG_BETA
    end

    # Test immediate destroy and double-destroy, and their status codes...
    p = finufft_makeplan(2,10,+1,1,1e-6);
    @test finufft_destroy(p)==0   # 0 signifies success.
    @test finufft_destroy(p)==1   # 1 since already destroyed; watch for crash
end


text:
# perfctr_gpu.jl
using LIKWID
using LinearAlgebra
using CUDA

@assert CUDA.functional()

GPUMarker.init()

# Note: CUDA defaults to Float32
Agpu = CUDA.rand(128, 64)
Bgpu = CUDA.rand(64, 128)
Cgpu = CUDA.zeros(128, 128)

GPUMarker.startregion("matmul")
for _ in 1:100
    mul!(Cgpu, Agpu, Bgpu)
end
GPUMarker.stopregion("matmul")

GPUMarker.close()


text:
"""
Based on:
Wilhelm Burger: Zhang’s Camera Calibration Algorithm: In-Depth Tutorial and Implementation,
Technical Report HGB16-05, University of Applied Sciences Upper Austria, School of
Informatics, Communications and Media, Dept. of Digital Media, Hagenberg, Austria, May
2016. DOI: 10.13140/RG.2.1.1166.1688, http://staff.fh-hagenberg.at/burger/
"""
module CameraCalibration
using Statistics, LinearAlgebra
export estimateHomography, calibrate, getHomographies, getCameraIntrinsics, getExtrinsics, getCameraIntrinsicsB
"""
estimateHomography - returns the estimated homography matrix H, such that b_j = H*a_j.
'''
A and B are in the form 3xN (each column is a point)
'''
"""
function estimateHomography(a::Array{<:Real,2}, b::Array{<:Real,2})
    if size(a) != size(b)
        error("src and dst must have same dimensions")
    end
    n = size(a)[2] #number of points
    #normalisation matrices
    nA = getNormalisationMatrix(a)
    nB = getNormalisationMatrix(b)
    m = zeros((2*n, 9))

    #construct matrix for solving
    for j  = 1:n
        k = 2*j
        ap = nA*a[:,j]
        bp = nB*b[:,j]
        # M_{k,*} = (x′_a, y′_a, 1, 0, 0, 0, −x′_ax′_b , −y′_ax′_b, −x′_b )
        m[k-1,:] = [ap[1] ap[2] 1 0 0 0 -ap[1]*bp[1] -ap[2]*bp[1] -bp[1]]
        # M_{k+1,*} = (0, 0, 0, x′_a , y′_a , 1, −x′_ay′_b, −y′_ay′_b , −y′_b )
        m[k,:] = [0 0 0 ap[1] ap[2] 1 -ap[1]*bp[2] -ap[2]*bp[2] -bp[2]]
    end

    #solve
    u, d, v = svd(m, full = true)
    h = v[:,9] #last column of v
    H = [h[1] h[2] h[3];
         h[4] h[5] h[6]
         h[7] h[8] h[9]]

    H = inv(nB)*H*nA
    return H./H[3,3]
end
"""
getNormalisationMatrix - returns normalisation matrix of points x
'''

'''
"""
function getNormalisationMatrix(x::Array{<:Real, 2})
    m = mean(x[1:2, :], dims = 2)
    v = var(x[1:2, :], dims = 2, corrected = false, mean = m)
    #scale values (182)
    sx = sqrt(2.0/v[1])
    sy = sqrt(2.0/v[2])
    #construct normalisation matrix (181)
    return [sx 0 -sx*m[1]; 0 sy -sy*m[2]; 0 0 1]
end

function refineHomography(h::Array{<:Real, 2}, a::Array{<:Real, 2}, b::Array{<:Real, 2})
    n = size(a)[2]
    X = zeros((2,2*n))
    Y = zoeros((2,2*n))
    for i = 1:n

    end
end

function homographyVal(x::Array{<:Real,2}, h::Array{<:Real, 1})
    n = size(x)[2]
    y = zeros(2*n)
    for j = 1:n
        w = h[7]*x[1,j] + h[8]*x[2,j]+h[9]
        y[j*2-1] = (h[1]*x[1,j] + h[2]*x[2,j]+h[3])/w
        y[j*2] = (h[4]*x[1,j] + h[5]*x[2,j]+h[6])/w
    end
    return y
end

function homographyJac(x::Array{<:Real,2}, h::Array{<:Real, 1})
    n = size(x)[2]
    J = zeros((2*n, 9))
    for j = 1:n
        sx = h[1]*x[1,j]+h[2]*x[2,j] + h[3]
        sy = h[4]*x[1,j]+h[5]*x[2,j] + h[6]
        w = h[7]*x[1,j]+h[8]*x[2,j] + h[9]
        J[2j-1,:]=[x[1,j]/w x[2,j]/w 1/w 0 0 0 -sx*x[1,j]/(w^2) -sx*x[2,j]/(w^2) -sx/(w^2)]
        J[2j,:]=[0 0 0 x[1,j]/w x[2,j]/w 1/w -sy*x[1,j]/(w^2) -sy*x[2,j]/(w^2) -sy/(w^2)]
    end
    return J
end

function calibrate(x::Array{<:Real,2}, u::Array{Array{T1,2},1}) where T1 <: Real
    hListInit = getHomographies(x,u)
    aInit = getCameraIntrinsicsB(hListInit)
    wInit = getExtrinsics(aInit,hListInit)
end

function getHomographies(x::Array{<:Real,2}, u::Array{Array{T1,2},1}) where T1 <: Real
    m = size(u)[1]
    hList = Array{Float64,2}[]
    for i = 1:m
        hInit = estimateHomography(x, u[i])
        # TODO refineHomography
        push!(hList, hInit)
    end
    return hList
end

function getCameraIntrinsics(hs::Array{Array{T1,2},1}) where T1 <: Real
    m = size(hs)[1]
    V = zeros((2*m, 6))

    for i = 1:m
        V[2*i-1,:] = getIntrinsicRowVector(1,2,hs[i])
        V[2*i,:] = getIntrinsicRowVector(1,1,hs[i]) - getIntrinsicRowVector(2,2,hs[i])
    end

    #solve
    u, s, v = svd(V, full = true)
    print(size(v))
    b = v[:,6]

    w = b[1]*b[3]*b[6] - (b[2]^2)*b[6] - b[1]*(b[5]^2) + 2*b[2]*b[4]*b[5] - b[3]*(b[4]^2) # (104)
    d = b[1]*b[3] - b[2]^2 # (105)
    α = sqrt(w/(d*b[1])) # (99)
    β = sqrt(w/d^2 *b[1]) # (100)
    γ = sqrt(w/(d^2 * b[1])) * b[2] # (101)
    uc = (b[2]*b[5] - b[3]*b[4])/d
    vc = (b[2]*b[4] - b[1]*b[5])/d
    return [α γ uc
            0 β vc
            0 0 1]
end

function getCameraIntrinsicsB(hs::Array{Array{T1,2},1}) where T1 <: Real
    m = size(hs)[1]
    V = zeros((2*m, 6))

    for i = 1:m
        V[2*i-1,:] = getIntrinsicRowVector(1,2,hs[i])
        V[2*i,:] = getIntrinsicRowVector(1,1,hs[i]) .- getIntrinsicRowVector(2,2,hs[i])
    end

    #solve
    u, s, v = svd(V, full = true)
    b = v[:,6]
    if b[1] < 0 || b[3] < 0 || b[6] < 0
        b = -b
    end
    bInit = [b[1] b[2] b[4]
             b[2] b[3] b[5]
             b[4] b[5] b[6]]
    C = cholesky(bInit)
    return inv(C.L)' * C.L[3,3]
end

function getIntrinsicRowVector(p::Int64, q::Int64,  h::Array{<:Real,2})
    return [h[1,p]*h[1,q]
            h[1,p]*h[2,q]+h[2,p]*h[1,q]
            h[2,p]*h[2,q]
            h[3,p]*h[1,q]+h[1,p]*h[3,q]
            h[3,p]*h[2,q]+h[2,p]*h[3,q]
            h[3,p]*h[3,q]]'
end

"""
a is intrinsic matrix hs is list of homographies

"""
function getExtrinsics(a::Array{<:Real,2},hs::Array{Array{T1,2},1}) where T1 <: Real
    wList = Array{Float64,2}[]
    m = size(hs)[1]

    for i = 1:m
        push!(wList, estimateViewTransorm(a, hs[i]))
    end
    return wList
end

function estimateViewTransorm(a::Array{<:Real,2}, h::Array{<:Real,2})
    κ = 1/norm(inv(a)*h[:,1])
    r0 = κ*inv(a)*h[:,1]
    r1 = κ*inv(a)*h[:,2]
    r2= cross(r0,r1)
    R = zeros((3, 3))
    R[:,1] = r0
    R[:,2] = r1
    R[:,3] = r2
    # make true rotation matrix
    u, s, v = svd(R, full = true)
    R = v*u'
    t = κ*inv(a)*h[:,3]
    return hcat(R,t)
end
end  # module CameraCalib

"""
tests
asample = [182 535 171 537; 350 358 553 563; 1 1 1 1]
bsample = [0 888 0 888; 0 0 500 500; 1 1 1 1]

testobj = [0. 1. 2. 3. 4. 5. 6. 0. 1. 2. 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6;
           4 4 4 4 4 4 4 3 3 3 3 3 3 3 2 2 2 2 2 2 2 1 1 1 1 1 1 1 0 0 0 0 0 0 0;
           1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

using CSV, DataFrames
imgpoints = Array{Float64, 2}[]
    for i = 0:10
        push!(imgpoints, vcat(convert(Array{Float64, 2}, CSV.File(string("./test/imgpoints", i, ".csv"), header = false) |> DataFrame)', ones(1, 35)))
    end

using Plots
cord = [-ex[i][:,1:3]'*ex[i][:,4] for i in 1:11]
p = scatter([cord[i][1] for i = 1:11], [cord[i][2] for i = 1:11], [cord[i][3] for i = 1:11],marker=:circle,linewidth=0, group = 1:11)
plot(p, xlabel="X",ylabel="Y",zlabel="Z", size = (800, 800))
"""


text:
using DifferentialEquations, Plots

𝞪 = 1
u0 = 1/2
f(t,u) = 𝞪*u
tspan = (0.0, 1.0)

prob = ODEProblem(f, u0, tspan)

sol1 = solve(prob, Euler(), dt=1/2^1)
sol2 = solve(prob, Euler(), dt=1/2^2)
sol3 = solve(prob, Euler(), dt=1/2^4)

#[t+2u for (t,u) in zip(sol1.t,sol.u)]

x = 0.0:0.05:1.0
p = plot(x, sol1(x))
plot!(p, x, sol2(x))
plot!(p, x, sol3(x))

gui()


text:
using Turing

#Generate some data
const data = [randn(10)+6, randn(10)-6]

@sf mixture_model() = begin
    #The model parameters
    K = 2
    ms = [normal(0, 10) for i = 1:K]
    vs = [gamma(2, 2) for i = 1:K]
    ps = dirichlet(ones(K))

    #Which mixture component does each data item belong to?
    ks = [categorical(ps) for i = 1:length(data)]

    for i = 1:length(data)
      #Condition on the data.
      normal(ms[ks[i]], vs[ks[i]]; condition=data[i])
    end
    (ms, vs, ps)
end
print(sample(mixture_model, 10^4)[end])


text:
using Cassette, InteractiveUtils

Cassette.@context TraceurCtx

struct Trace
  seen::Set
  warn
end

Trace(w) = Trace(Set(), w)

isprimitive(f) = f isa Core.Builtin || f isa Core.IntrinsicFunction

const ignored_methods = Set([@which((1,2)[1])])
const ignored_functions = Set([getproperty, setproperty!])

function Cassette.posthook(ctx::TraceurCtx, out, f, args...)
  C, T = DynamicCall(f, args...), typeof.((f, args))
  tra = ctx.metadata
  (f ∈ ignored_functions || T ∈ tra.seen || isprimitive(f) ||
    method(C) ∈ ignored_methods || method(C).module ∈ (Core, Core.Compiler)) && return nothing

  push!(tra.seen, T)
  analyse((a...) -> tra.warn(Warning(a...)), C)
  return nothing
end

trace(w, f) = Cassette.overdub(TraceurCtx(metadata=Trace(w)), f)

warntrace(f) = trace(warning_printer(), f)

function warnings(f)
  warnings = Warning[]
  trace(w -> push!(warnings, w), f)
  return warnings
end

macro trace(ex)
  :(warntrace(() -> $(esc(ex))))
end


text:
"""
    StateCheck

Module with a minimal set of functions for getting statistics
and basic I/O from ClimateMachine DG state arrays (`MPIStateArray` type).
Created for regression testing and code change tracking and debugging.
StateCheck functions iterate over named variables in an `MPIStateArray`,
calculate and report their statistics and/or write values for all or
some subset of points at a fixed frequency.

# Functions

 - [`sccreate`](@ref) Create a StateCheck callback variable.
 - [`scdocheck`](@ref) Check StateCheck variable values against reference values.
 - [`scprintref`](@ref) Print StateCheck variable in format for creating reference values.
"""
module StateCheck

# Imports from standard Julia packages
using Formatting
using MPI
using ..MPIStateArrays: vars
using PrettyTables
using Printf
using StaticArrays
using Statistics
using ..VariableTemplates:
    flattenednames, flattened_tup_chain, RetainArr, varsindex, varsize

# Imports from ClimateMachine core
import ..GenericCallbacks: EveryXSimulationSteps
import ..MPIStateArrays: MPIStateArray

"""
    VStat

Type for returning variable statistics.
"""
struct VStat
    max::Any
    min::Any
    mean::Any
    std::Any
end

# Global functions to expose
export sccreate # Create a state checker callback
export scdocheck
export scprintref

const nt_freq_def = 10 # default frequency (in time steps) for output.
const prec_def = 15    # default precision used for formatted output table

# TODO: this should use the new callback interface

"""
    sccreate(
        io::IO,
        fields::Array{<:Tuple{<:MPIStateArray, String}},
        nt_freq::Int = $nt_freq_def;
        prec = $prec_def
    )

Create a "state check" call-back for one or more `MPIStateArray` variables
that will report basic statistics for the fields in the array.

  -  `io` an IO stream to use for printed output
  -  `fields` a required first argument that is an array of one or more
                        `MPIStateArray` variable and label string pair tuples.
                        State array statistics will be reported for the named symbols
                        in each `MPIStateArray` labeled with the label string.
  -  `nt_freq` an optional second argument with default value of
                        $nt_freq_def that sets how frequently (in time-step counts) the
                        statistics are reported.
  -  `prec` a named argument that sets number of decimal places to print for
                        statistics, defaults to $prec_def.

# Examples
```julia
using ClimateMachine.VariableTemplates
using StaticArrays
using ClimateMachine.MPIStateArrays
using MPI
MPI.Init()
FT=Float64
F1=@vars begin; ν∇u::SMatrix{3, 2, FT, 6}; κ∇θ::SVector{3, FT}; end
F2=@vars begin; u::SVector{2, FT}; θ::SVector{1, FT}; end
Q1=MPIStateArray{Float32,F1}(MPI.COMM_WORLD,ClimateMachine.array_type(),4,9,8);
Q2=MPIStateArray{Float64,F2}(MPI.COMM_WORLD,ClimateMachine.array_type(),4,6,8);
cb=ClimateMachine.StateCheck.sccreate([(Q1,"My gradients"),(Q2,"My fields")],1; prec=$prec_def);
```
"""
function sccreate(
    io::IO,
    fields::Array{<:Tuple{<:MPIStateArray, String}},
    nt_freq::Int = nt_freq_def;
    prec = prec_def,
    print_head = true,
)
    mpicomm = first(first(fields)).mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    mpirank == 0 && println(io, "# SC Start: creating state check callback")

    ####
    # Print fields that the callback created by this call will query
    ####
    for (Q, lab) in fields
        V = vars(Q)
        ftc = flattened_tup_chain(V, RetainArr())
        nss = length(flattenednames(V))
        if length(flattenednames(V)) == 0 && mpirank == 0
            println(
                io,
                "# SC  MPIStateArray labeled \"$lab\" has no named symbols.",
            )
        elseif mpirank == 0
            println(
                io,
                "# SC Creating state check callback labeled \"$lab\" for symbols",
            )
            for s in join.(ftc, "")
                println(io, "# SC ", s)
            end
        end
    end

    ###
    # Initialize total calls counter for the callback
    ###
    n_cb_calls = 0

    ###
    # Create holder for most recent stats
    ###
    nvars = sum([length(flattenednames(vars(Q))) for (Q, lab) in fields])
    cur_stats_flat = Array{Any}(undef, nvars)

    # Save io pointer
    iosave = io

    ######
    # Create the callback
    ######
    cb = EveryXSimulationSteps(nt_freq) do
        # Track which timestep this is
        n_cb_calls = n_cb_calls + 1
        n_step = (n_cb_calls - 1) * nt_freq + 1
        ns_str = @sprintf("%7.7d", n_step - 1)
        io = iosave
        # Obscure trick to do with running in cells in notebook that close Base.stdout
        # between different stages. This affects parsing Literate/Documenter examples.
        if !isopen(io)
            io = Base.stdout
        end

        ## Print header
        nprec = min(max(1, prec), 20)
        if mpirank == 0
            println(
                io,
                "# SC +++++++++++ClimateMachine StateCheck call-back start+++++++++++++++++",
            )
            println(io, "# SC Step = $ns_str")
        end

        labs = vcat([
            map(x -> lab, flattenednames(vars(Q))) for (Q, lab) in fields
        ]...)

        varnames = vcat([
            map(x -> x, flattenednames(vars(Q))) for (Q, lab) in fields
        ]...)

        header = ["Label" "Field" "min()" "max()" "mean()" "std()"]

        stats_all = vcat([
            [scstats(Q, i, nprec)[5] for i in 1:varsize(vars(Q))] for (Q, lab) in fields
        ]...)

        cur_stats_flat .= vcat([
            [
                begin
                    scstats_i = scstats(Q, i, nprec)[5]
                    [
                        lab,
                        fn,
                        scstats_i.max,
                        scstats_i.min,
                        scstats_i.mean,
                        scstats_i.std,
                    ]
                end for (i, fn) in enumerate(flattenednames(vars(Q)))
            ] for (Q, lab) in fields
        ]...)

        # TODO: Verify, not sure why we are swapping these columns:
        # data_min = map(x->x.min, stats_all)
        # data_max = map(x->x.max, stats_all)

        data_max = map(x -> x.min, stats_all)
        data_min = map(x -> x.max, stats_all)

        data_mean = map(x -> x.mean, stats_all)
        data_std = map(x -> x.std, stats_all)
        data = hcat(labs, varnames, data_min, data_max, data_mean, data_std)
        pretty_table(
            io,
            data,
            header;
            formatters = ft_printf("%.16e", 3:6),
            header_crayon = crayon"yellow bold",
            subheader_crayon = crayon"green bold",
            crop = :none,
        )

        if mpirank == 0
            println(
                io,
                "# SC +++++++++++ClimateMachine StateCheck call-back end+++++++++++++++++++",
            )
        end
    end

    if mpirank == 0
        println(io, "# SC Finish: creating state check callback")
    end

    return cb
end

function scstats(Q, ivar, nprec)

    # Get number of MPI procs
    mpicomm = Q.mpicomm
    nproc = MPI.Comm_size(mpicomm)

    npr = nprec
    fmt = @sprintf("%%%d.%de", npr + 8, npr)

    getByField(Q, ivar) = (Q.realdata[:, ivar, :])

    # Get local and global field sizes (degrees of freedom).
    n_size_loc = length(getByField(Q, ivar))
    n_size_tot = MPI.Allreduce(n_size_loc, +, mpicomm)

    # Min
    phi_loc = minimum(getByField(Q, ivar))
    phi_min = MPI.Allreduce(phi_loc, min, mpicomm)
    phi = phi_min
    # minVstr=@sprintf("%23.15e",phi)
    min_v_str = sprintf1(fmt, phi)

    # Max
    phi_loc = maximum(getByField(Q, ivar))
    phi_max = MPI.Allreduce(phi_loc, max, mpicomm)
    phi = phi_max
    # maxVstr=@sprintf("%23.15e",phi)
    max_v_str = sprintf1(fmt, phi)

    # Ave
    phi_loc = mean(getByField(Q, ivar))
    phi_loc = phi_loc * n_size_loc / n_size_tot
    phi_sum = MPI.Allreduce(phi_loc, +, mpicomm)
    phi_mean = phi_sum
    phi = phi_mean
    # aveVstr=@sprintf("%23.15e",phi)
    ave_v_str = sprintf1(fmt, phi)

    # Std
    phi_loc = (getByField(Q, ivar) .- phi_mean) .^ 2
    phi_loc = sum(phi_loc)  # Sum local data explicitly since GPU Allreduce
    # does not take arrays yet.
    phi_sum = MPI.Allreduce(phi_loc, +, mpicomm)
    n_val_sum = n_size_tot
    phi_std = sqrt(sum(phi_sum) / (n_val_sum - 1))
    phi = phi_std
    # stdVstr=@sprintf("%23.15e",phi)
    std_v_str = sprintf1(fmt, phi)

    vals = VStat(phi_min, phi_max, phi_mean, phi_std)

    return min_v_str, max_v_str, ave_v_str, std_v_str, vals
end

function sccreate(
    fields::Array{<:Tuple{<:MPIStateArray, String}},
    nt_freq::Int = nt_freq_def;
    prec = prec_def,
)
    return sccreate(Base.stdout, fields, nt_freq; prec = prec)
end

"""
    scprintref(cb)

Print out a "state check" call-back table of values in a format
suitable for use as a set of reference numbers for CI comparison.

 - `cb` callback variable of type ClimateMachine.GenericCallbacks.Every*
"""
function scprintref(cb)
    sc = cb.callback
    io = sc.iosave
    # Obscure trick to do with running in cells in notebook
    if !isopen(io)
        io = Base.stdout
    end
    mpirank = MPI.Comm_rank(MPI.COMM_WORLD)
    mpirank == 0 || return
    # Get print format lengths for cols 1 and 2 so they are aligned
    # for readability.
    phi = sc.cur_stats_flat
    a1l = maximum(length.(map(x -> x[1], phi)))
    a2l = maximum(length.(String.(map(x -> x[2], phi))))
    fmt1 = @sprintf("%%%d.%ds", a1l, a1l) # Column 1
    fmt2 = @sprintf("%%%d.%ds", a2l, a2l) # Column 2
    fmt3 = @sprintf("%%28.20e")         # All numbers at full precision
    # Create an string of spaces to be used for formatting
    sp = repeat(" ", 75)

    # Write header
    println(io, "# BEGIN SCPRINT")
    println(io, "# varr - reference values (from reference run)    ")
    println(io, "# parr - digits match precision (hand edit as needed) ")
    println(io, "#")
    println(io, "# [")
    println(
        io,
        "#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],",
    )
    println(
        io,
        "#  [         :                :          :        :      :          :           ],",
    )
    println(io, "# ]")
    #
    # Write tables
    #  Reference value and precision match tables are separate since it is more
    #  common to update reference values occasionally while precision values are
    #  typically changed rarely and the precision values are hand edited from experience.
    #
    # Write table of reference values
    println(io, "varr = [")
    for (s1, s2, s3′, s4′, s5′, s6′) in sc.cur_stats_flat
        s1 = sp[1:(a1l - length(s1))] * "\"" * s1 * "\""
        s2 = sp[1:(a2l - length(s2))] * "\"" * s2 * "\""
        s3 = sprintf1(fmt3, s3′)
        s4 = sprintf1(fmt3, s4′)
        s5 = sprintf1(fmt3, s5′)
        s6 = sprintf1(fmt3, s6′)
        srow = (s3, s4, s5, s6)
        println(io, " [ ", s1, ", ", s2, ", ", join(srow, ","), " ],")
    end
    println(io, "]")

    # Write table of reference match precisions using default precision that
    # can be hand updated.
    println(io, "parr = [")
    for (s1, s2, s3′, s4′, s5′, s6′) in sc.cur_stats_flat
        s1 = sp[1:(a1l - length(s1))] * "\"" * s1 * "\""
        s2 = sp[1:(a2l - length(s2))] * "\"" * s2 * "\""
        println(io, " [ ", s1, ", ", s2, ",", "    16,    16,    16,    16 ],")
    end
    println(io, "]")
    println(io, "# END SCPRINT")
end

# TODO: write unit test for process_stat
function process_stat(row)
    ## Debugging
    # println(row)
    # println(ref_dat_val)
    # println(ref_dat_prec)
    row_col_pass = 0
    row_col_na = 0
    row_col_fail = 0

    ## Make array copy for reporting

    ## Check MPIStateArrayName
    cval = row.lab.cur
    rval = row.lab.ref_val
    if cval != rval
        row_col_fail += 1
        lab = "N" * "(" * rval * ")"
    else
        lab = cval
        row_col_pass += 1
        lab = rval
    end

    ## Check term name
    cval = row.varname.cur
    rval = row.varname.ref_val
    if cval != rval
        varname = "N" * "(" * rval * ")"
        row_col_fail += 1
    else
        varname = cval
        row_col_pass += 1
    end

    res_dat = Dict()
    # Check numeric values
    for nv in (:min, :max, :mean, :std)
        fmt = @sprintf("%%28.20e")
        lfld = 28
        cval = getproperty(row, nv).cur
        cvalc = sprintf1(fmt, cval)
        rval = getproperty(row, nv).ref_val
        rvalc = sprintf1(fmt, rval)
        pcmp = getproperty(row, nv).ref_prec

        # Skip if compare digits set to 0
        if pcmp > 0

            # Check exponent part
            ep1 = cvalc[(lfld - 3):lfld]
            ep2 = rvalc[(lfld - 3):lfld]
            if ep1 != ep2
                nmatch = 0
            else
                # Now check individual digits left to right
                dp1 = cvalc[2:(3 + pcmp + 1)]
                dp2 = rvalc[2:(3 + pcmp + 1)]
                nmatch = 0
                imatch = 1
                for c in dp1
                    if c == dp2[imatch]
                        nmatch = imatch
                    else
                        break
                    end
                    imatch += 1
                end
            end
            # Check for trailing exact 0s number change numerically
            if nmatch < pcmp
                if rval != 0
                    e_diffl10 = round(log10(abs((rval - cval) / rval)))
                else
                    e_diffl10 = round(log10(abs(rval - cval)))
                end
                if e_diffl10 < -pcmp
                    nmatch = Int(-e_diffl10)
                end
            end
            if nmatch < pcmp
                res_dat[nv] = "N(" * string(nmatch) * ")"
                row_col_fail += 1
            else
                res_dat[nv] = string(nmatch)
                row_col_pass += 1
            end
        else
            res_dat[nv] = "0"
            row_col_na += 1
        end
    end
    na_count = row_col_na
    fail_count = row_col_fail
    pass_count = row_col_pass
    return (; lab, varname, res_dat..., pass_count, fail_count, na_count)
end

"""
    scdocheck(cb, ref_dat)

Compare a current State check call-back set of values with a
reference set and match precision table pair.

 - `cb` `StateCheck` call-back variables
 - `ref_dat` an array of reference values and precision to match tables.
"""
function scdocheck(cb, ref_dat)
    sc = cb.callback
    io = sc.iosave
    # Obscure trick to do with running in cells in notebook
    if !isopen(io)
        io = Base.stdout
    end
    mpirank = MPI.Comm_rank(MPI.COMM_WORLD)
    mpirank == 0 || return true
    println(
        io,
        "# SC +++++++++++ClimateMachine StateCheck ref val check start+++++++++++++++++",
    )
    println(
        io,
        "# SC \"N( )\" bracketing indicates field failed to match      ",
    )

    get_row(cur, ref_prec, ref_val, i) =
        (; cur = cur[i], ref_prec = ref_prec[i], ref_val = ref_val[i])

    Z = (sc.cur_stats_flat, ref_dat...)
    stats_all = [
        begin
            (;
                lab = get_row(row, row_ref_prec, row_ref_val, 1),
                varname = get_row(row, row_ref_prec, row_ref_val, 2),
                min = get_row(row, row_ref_prec, row_ref_val, 3),
                max = get_row(row, row_ref_prec, row_ref_val, 4),
                mean = get_row(row, row_ref_prec, row_ref_val, 5),
                std = get_row(row, row_ref_prec, row_ref_val, 6),
            )
        end for (row, row_ref_val, row_ref_prec) in zip(Z...)
    ]

    data_all = map(row -> process_stat(row), stats_all)

    labs = map(stats -> stats.lab, data_all)
    varnames = map(stats -> stats.varname, data_all)
    data_min = map(stats -> stats.min, data_all)
    data_max = map(stats -> stats.max, data_all)
    data_mean = map(stats -> stats.mean, data_all)
    data_std = map(stats -> stats.std, data_all)
    pass_count = map(stats -> stats.pass_count, data_all)
    fail_count = map(stats -> stats.fail_count, data_all)
    na_count = map(stats -> stats.na_count, data_all)
    all_pass = sum(fail_count) == 0

    data = hcat(
        labs,
        varnames,
        data_min,
        data_max,
        data_mean,
        data_std,
        pass_count,
        fail_count,
        na_count,
    )
    header = [
        "Label" "Field" "min()" "max()" "mean()" "std()" "Pass" "Fail" "Not checked"
        "" "" "" "" "" "" "count" "count" "count"
    ]

    pretty_table(
        io,
        data,
        header;
        header_crayon = crayon"yellow bold",
        subheader_crayon = crayon"green bold",
        crop = :none,
    )

    println(
        io,
        "# SC +++++++++++ClimateMachine StateCheck ref val check end+++++++++++++++++",
    )
    return all_pass

end

end # module


text:
import Molecules: translate, Cn, σ, Sn, i, transform, isequivalent

@testset "Transformations" begin
    @testset "Translate" begin
        atoms = [
            Atom(1, 1.008, [1.0, 2.0, 3.0])
        ]

        natoms = translate(atoms, [1.0, 2.0, 3.0])

        @test natoms[1].xyz == [0.0, 0.0, 0.0]
        end

    @testset "Operators" begin
        test = [[-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0],
                [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 1.0],
                [1.0 -0.0 -0.0; -0.0 1.0 -0.0; -0.0 -0.0 -1.0],
                [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 -1.0],
                Int8[-1 0 0; 0 -1 0; 0 0 -1]]
        C2 = Cn([0.0, 0.0, 1.0], 2)
        C3 = Cn([0.0, 0.0, 1.0], 3)
        σxy = σ([0.0, 0.0, 1.0])
        S3 = Sn([0.0, 0.0, 1.0], 3)
        inv = i()
        out = [C2, C3, σxy, S3, inv]
        for i in 1:5
            @test isapprox(out[i], test[i], rtol=1E-5)
        end

    end

    @testset "Transform" begin
        H2 = [Atom(1, 1.008, [-0.5, 0.0, 0.0]), Atom(1, 1.008, [0.5, 0.0, 0.0])]
        H2O = [Atom(1, 1.008, [-0.5, 0.0, 0.0]), Atom(1, 1.008, [0.5, 0.0, 0.0]), Atom(1, 1.008, [0.0, 1.0, 0.0])]
        ops = [Cn([0.0, 1.0, 0.0], 2), Cn([0.0, 1.0, 1.0], 2), Cn([1.0, 1.0, 0.0], 2), σ([1.0, 0.0, 0.0]), σ([1.0, 1.0, 0.0]), Sn([1.0, 0.0, 0.0], 5), Sn([1.0, 1.0, 0.0], 3), i()]
        H2ret = [true, true, false, true, false, true, false, true]
        H2Oret = [true, false, false, true, false, false, false, false]
        for i = 1:size(ops)[1]
            H2b = transform(H2, ops[i])
            H2Ob = transform(H2O, ops[i])
            @test isequivalent(H2, H2b) == H2ret[i]
            @test isequivalent(H2O, H2Ob) == H2Oret[i]
        end

    end
end


text:
@enum ResponseState NotStarted OnMessageBegin StatusComplete OnHeaderField OnHeaderValue HeadersDone OnBody BodyDone EarlyEOF UpgradeConnection

type ResponseStream{T<:IO} <: IO
    response::Response
    socket::T
    state::ResponseState
    buffer::IOBuffer
    parser::ResponseParser
    timeout::Float64
    current_header::Nullable{Compat.UTF8String}
    state_change::Condition
    cookie_buffer::IOBuffer
    ResponseStream() = new()
end

function ResponseStream{T}(response, socket::T)
    r = ResponseStream{T}()
    r.response = response
    r.socket = socket
    r.state = NotStarted
    r.buffer = IOBuffer()
    r.parser = ResponseParser(r)
    r.timeout = Inf
    r.current_header = Nullable()
    r.state_change = Condition()
    r.cookie_buffer = IOBuffer()
    r
end

function Base.show(io::IO, stream::ResponseStream)
    print(io, "ResponseStream($(requestfor(stream)))")
end

function send_headers(response_stream::ResponseStream)
    request = requestfor(response_stream)
    print(response_stream, request.method, " ", isempty(request.resource) ? "/" : request.resource,
          " HTTP/1.1", CRLF,
          map(h->string(h,": ",request.headers[h],CRLF), collect(keys(request.headers)))...,
          "", CRLF)
    response_stream
end

immutable TimeoutException <: Exception
    timeout::Float64
end

function Base.show(io::IO, err::TimeoutException)
    print(io, "TimeoutException: server did not respond for more than $(err.timeout) seconds. ")
end

function process_response(stream)
    rp = stream.parser
    last_received = now()
    status_channel = Channel{Symbol}(1)
    timeout = stream.timeout
    if timeout < Inf
        Timer(0, timeout) do timer
            delta = now() - last_received
            if timeout_in_sec(delta) > timeout
                close(timer)
                put!(status_channel, :timeout)
            end
        end
    end
    @schedule begin
        while stream.state < BodyDone && !eof(stream.socket)
            last_received = now()
            data = readavailable(stream.socket)
            if length(data) > 0
                add_data(rp, data)
            end
        end
        put!(status_channel, :success)
    end
    status = take!(status_channel)
    status == :timeout && throw(TimeoutException(timeout))
    if stream.state < BodyDone
        stream.state = EarlyEOF
        notify(stream.state_change)
    end
    stream
end

function tls_dbg(level, filename, number, msg)
    warn("MbedTLS emitted debug info: $msg in $filename:$number")
end

function get_default_tls_config(verify=true)
    conf = MbedTLS.SSLConfig()
    MbedTLS.config_defaults!(conf)

    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)
    MbedTLS.rng!(conf, rng)

    MbedTLS.authmode!(conf,
      verify ? MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED : MbedTLS.MBEDTLS_SSL_VERIFY_NONE)
    MbedTLS.dbg!(conf, tls_dbg)
    MbedTLS.ca_chain!(conf)

    conf
end

function Base.wait(stream::ResponseStream)
    stream.state >= BodyDone && return
    wait(stream.state_change)
end

function Base.eof(stream::ResponseStream)
    eof(stream.buffer) && (stream.state==BodyDone || eof(stream.socket))
end

for T in [BitArray, Vector{UInt8}, UInt8]
    @eval Base.write(stream::ResponseStream, data::$T) = write(stream.socket, data)
end

function Base.readbytes!(stream::ResponseStream, data::Vector{UInt8}, sz)
    while stream.state < BodyDone && nb_available(stream) < sz
        wait(stream)
    end
    readbytes!(stream.buffer, data, sz)
end

function Base.read(stream::ResponseStream)
    while stream.state < BodyDone
        wait(stream)
    end
    take!(stream.buffer)
end

function Base.read(stream::ResponseStream, ::Type{UInt8})
    while !eof(stream) && nb_available(stream) < 1
        wait(stream)
    end
    read(stream.buffer, UInt8)
end

function Base.readavailable(stream::ResponseStream)
    while nb_available(stream) == 0
        wait(stream)
    end
    read(stream, nb_available(stream))
end

Base.close(stream::ResponseStream) = close(stream.socket)
Base.nb_available(stream::ResponseStream) = nb_available(stream.buffer)

for getter in [:headers, :cookies, :statuscode, :requestfor, :history]
    @eval $getter(stream::ResponseStream) = $getter(stream.response)
end

Base.convert(::Type{Response}, stream::ResponseStream) = stream.response

immutable ProxyException <: Exception
    resp::Response
end

function Base.show(io::IO, err::ProxyException)
    print(io, "Failed to open CONNECT tunnel on proxy server. Proxy response was $(err.resp)")
end

function open_stream(req::Request, tls_conf=TLS_VERIFY, timeout=Inf,
                     http_proxy=Nullable{URI}(), https_proxy=Nullable{URI}())
    uri = req.uri
    connect_method = :direct
    if scheme(uri) == "http"
        if !isnull(http_proxy)
            uri = get(http_proxy)
            connect_method = :tunnel
        end
    elseif scheme(uri) == "https"
        if !isnull(https_proxy)
            uri = get(https_proxy)
            connect_method = :tunnel
        end
    else
        error("Unsupported scheme \"$(scheme(uri))\"")
    end
    if connect_method == :tunnel
        req.resource = "$(req.uri.scheme)://$(req.uri.host)$(resourcefor(req.uri))"
    end
    ip = Base.getaddrinfo(uri.host)
    if scheme(req.uri) == "http"
        stream = Base.connect(ip, http_port(uri))
    else
        # Initialize HTTPS
        if connect_method == :tunnel
            sock = Base.connect(ip, http_port(uri))
            tunnel_req = Request()
            tunnel_req.method = "CONNECT"
            tunnel_resp = Response()
            tunnel_resp.request = Nullable(tunnel_req)
            tunnel_resp_stream = ResponseStream(tunnel_resp, sock)
            tunnel_resp_stream.timeout = timeout
            write(sock, "CONNECT $(req.uri.host):$(https_port(req.uri)) HTTP/1.1\r\n\r\n")
            process_response(tunnel_resp_stream)
            if statuscode(tunnel_resp) ≠ 200
                throw(ProxyException(tunnel_resp))
            end
        else
            sock = Base.connect(ip, uri.port == 0 ? 443 : uri.port)
        end
        stream = MbedTLS.SSLContext()
        MbedTLS.setup!(stream, tls_conf)
        MbedTLS.set_bio!(stream, sock)
        MbedTLS.hostname!(stream, req.uri.host)
        MbedTLS.handshake(stream)
    end
    resp = Response()
    empty!(resp.headers)
    resp.request = Nullable(req)
    stream = ResponseStream(resp, stream)
    stream.timeout = timeout
    send_headers(stream)
    stream
end

function __init_streaming__()
    global const TLS_VERIFY = get_default_tls_config(true)
    global const TLS_NOVERIFY = get_default_tls_config(false)
end


text:
using NPZ
include("Lattice2D.jl")

data_path = "../../data"
if !isdir(data_path)
	mkpath(data_path)
end

matrix_to_tensor(matrix) = reshape(convert(Array{Int8}, matrix), size(matrix)..., 1)

function generate_data(
	L,
	n_samples,
	p_occupations;
	extra_sample_range=(0.55,0.65),
	extra_sample_rate=10,
	as_dict=false
	)

	N_samples = n_samples * length(p_occupations)
	extra = 0
	if extra_sample_rate > 1
		extra = Int(
			sum(extra_sample_range[1] .≤ p_occupations .≤ extra_sample_range[2]) * (extra_sample_rate - 1)
			)
	end
	N_samples += extra * n_samples

	X = zeros(Int8, L, L, 1, N_samples)
	Y = zeros(Int8, N_samples)
	p = zeros(N_samples)

	i = 1
	for p_occupation in p_occupations
		if 0.55 ≤ p_occupation ≤ 0.65
			n_samples_effective = n_samples * extra_sample_rate
		else
			n_samples_effective = n_samples
		end

		for sample in 1:n_samples_effective
			lattice = Lattice2D(L, p_occupation, i)
			X[:, :, :, i] = matrix_to_tensor(lattice.occupied)
			Y[i] = lattice.is_percolating
			p[i] = p_occupation

			i += 1
		end
	end

	if as_dict
		return Dict("X" => X, "Y" => Y, "p" => p)
	else
		return (X = X, Y = Y, p = p)
	end
end

# %% train and test data generation

L = 4

n_samples = 25
p_occupations = [
	range(0, stop=0.54, length=55);
	range(0.55, stop=0.65, length=51);
	range(0.66, stop=1, length=35)
	];

@time data_train_test  = generate_data(L, n_samples, p_occupations, as_dict=true);
size(data_train_test["X"])
npzwrite(joinpath(data_path, "data_train_test_$(L).npz"), data_train_test)

n_samples = 5000
p_occupations = collect(range(0.5, stop=0.65, step=1e-2));

@time data_transition = generate_data(L, n_samples, p_occupations, extra_sample_rate=1, as_dict=true);
size(data_transition["X"])
npzwrite(joinpath(data_path, "data_transition_$(L).npz"), data_transition)

# %%

L = 8

n_samples = 50
p_occupations = [
	range(0, stop=0.54, length=55);
	range(0.55, stop=0.65, length=51);
	range(0.66, stop=1, length=35)
	];

@time data_train_test  = generate_data(L, n_samples, p_occupations, as_dict=true);
size(data_train_test["X"])
npzwrite(joinpath(data_path, "data_train_test_$(L).npz"), data_train_test)

n_samples = 5000
p_occupations = collect(range(0.5, stop=0.65, step=1e-2));

@time data_transition = generate_data(L, n_samples, p_occupations, extra_sample_rate=1, as_dict=true);
size(data_transition["X"])
npzwrite(joinpath(data_path, "data_transition_$(L).npz"), data_transition)

# %%

L = 16

n_samples = 100
p_occupations = [
	range(0, stop=0.54, length=55);
	range(0.55, stop=0.65, length=51);
	range(0.66, stop=1, length=35)
	];

@time data_train_test  = generate_data(L, n_samples, p_occupations, as_dict=true);
size(data_train_test["X"])
npzwrite(joinpath(data_path, "data_train_test_$(L).npz"), data_train_test)

n_samples = 5000
p_occupations = collect(range(0.5, stop=0.65, step=1e-2));

@time data_transition = generate_data(L, n_samples, p_occupations, extra_sample_rate=1, as_dict=true);
size(data_transition["X"])
npzwrite(joinpath(data_path, "data_transition_$(L).npz"), data_transition)

# %%

L = 32

n_samples = 150
p_occupations = [
	range(0, stop=0.54, length=55);
	range(0.55, stop=0.65, length=51);
	range(0.66, stop=1, length=35)
	];

@time data_train_test  = generate_data(L, n_samples, p_occupations, as_dict=true);
size(data_train_test["X"])
npzwrite(joinpath(data_path, "data_train_test_$(L).npz"), data_train_test)

n_samples = 5000
p_occupations = collect(range(0.5, stop=0.65, step=1e-2));

@time data_transition = generate_data(L, n_samples, p_occupations, extra_sample_rate=1, as_dict=true);
size(data_transition["X"])
npzwrite(joinpath(data_path, "data_transition_$(L).npz"), data_transition)

max_stars_repo_name: cameronperot/percolation-classification

max_issues_repo_name: cameronperot/percolation-classification

max_forks_repo_name: cameronperot/percolation-classification


text:
"""
ДАНО: Робот - в произвольной клетке поля (с внутренними перегородками и маркерами)
РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы
"""
function new_perimeter(r::Robot)
    steps=0
    arr_of_steps=[]
    while ((isborder(r,Sud))&&(isborder(r,West)))==false
        push!(arr_of_steps,moves!(r,West))
        push!(arr_of_steps,moves!(r,Sud))
        steps+=2
    end
    side=Nord
    for a ∈ 1:4
        while isborder(r,side)==false
            putmarker!(r)
            move!(r,side)
        end
        side=inverse(next(side))
    end
    while (steps>0)==true
        side=isodd(steps) ? Ost : Nord
        for a ∈ 1:arr_of_steps[steps]
            move!(r,side)
        end
        steps-=1
    end
end
function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function next(side::HorizonSide)
    (HorizonSide(mod(Int(side)+1,4)))
end
function inverse(side::HorizonSide)
    (HorizonSide(mod(Int(side)+2,4)))
end
#Функции аналогичны тем, которые описаны в Задаче 5

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
#Demo 1: A full calculation/plotting example using EasyPlotInspect
#-------------------------------------------------------------------------------
module CMDimData_SampleUsage

using CMDimData
using CMDimData.MDDatasets
using CMDimData.EasyPlot
using CMDimData.EasyData
CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
WIDTH_LEGEND = 200
LBLAX_POS = "Position [m]"
LBLAX_TIME = "Time [s]"
pvst = cons(:a, labels = set(xaxis="Time (s)", yaxis="Position (m)"))
lstylesweep = cons(:a, line = set(width=3))
lstyle1 = cons(:a, line = set(color=:red, width=3, style=:solid))
lstyle2 = cons(:a, line = set(color=:blue, width=3))

#==Plot config (make legend more readable)
===============================================================================#
function adjust_legend(gtkplot)
	mplot = gtkplot.src
	for sp in mplot.subplots
		sp.layout[:halloc_legend] = WIDTH_LEGEND
	end
end
plotdisplay = EasyPlot.GUIDisplay(:InspectDR, postproc=adjust_legend)

#==Input data
===============================================================================#
x = DataF1(0:.1:20)
#NOTE: Both x & y coordinates of "x" object initialized as y = x = [supplied range]

#==Computations
===============================================================================#
#“Extract” maximum x-value from data:
xmax = maximum(x)

#Construct a “normalized” ramp dataset, unity_ramp:
unity_ramp = x/xmax

#Compute sin(x)
sinx = sin(x)

#Compute ramps with different slopes using unity_ramp (previously computed):
#NOTE: for Inner-most sweep, we need to specify leaf element type (DataF1 here):
ramp = fill(DataRS{DataF1}, PSweep("slope", [0, 0.5, 1, 1.5, 2])) do slope
	return unity_ramp * slope
end

#Merge two datasets with different # of sweeps (sinx & ramp):
r_sin = sinx+ramp

#Shift all ramped sin(x) waveforms to make them centered at their mid-points:
midval = (minimum(ramp) + maximum(ramp)) / 2

#Shift by midval (different for each swept slope of "ramp"):
c_sin = r_sin - midval

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=4,
	ystrip1 = set(axislabel=LBLAX_POS, striplabel="Plot: Building blocks"),
	ystrip2 = set(axislabel=LBLAX_POS, striplabel="Plot: Ramp with swept \"slope\""),
	ystrip3 = set(axislabel=LBLAX_POS, striplabel="Plot: sin(x) + ramp"),
	ystrip4 = set(axislabel=LBLAX_POS, striplabel="Plot: sin(x) + ramp - midval"),
	xaxis = set(label=LBLAX_TIME),
)

push!(plot,
	cons(:wfrm, unity_ramp, lstyle1, label="Unity Ramp", strip=1),
	cons(:wfrm, sinx, lstyle2, label="sin(x)", strip=1),
	cons(:wfrm, ramp, lstylesweep, label="swept ramp", strip=2),
	cons(:wfrm, r_sin, lstylesweep, label="r_sin", strip=3),
	cons(:wfrm, c_sin, lstylesweep, label="r_sin", strip=4),
)

pcoll = cons(:plot_collection, title="Simple Parametric Analysis", ncolumns=1)
	push!(pcoll, plot)

	#Save pcoll for later use:
	filename = basename(@__FILE__)
	savefile = joinpath("./", splitext(filename)[1] * ".hdf5")
	EasyData.writeplot(savefile, pcoll)

	#Display pcoll:
	#EasyPlot.displaygui(:InspectDR, pcoll)
	plotgui = display(plotdisplay, pcoll) #Tweak plot appearance with `plotdisplay`
	savefile = joinpath("./", splitext(filename)[1] * ".png")
	EasyPlot._write(:png, savefile, plotgui, dim=set(w=640, h=480))
end #module
:SampleCode_Executed


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    SemidiscretizationHyperbolic

A struct containing everything needed to describe a spatial semidiscretization
of a hyperbolic conservation law.
"""
struct SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition, BoundaryConditions,
                                    SourceTerms, Solver, Cache} <: AbstractSemidiscretization

  mesh::Mesh
  equations::Equations

  # This guy is a bit messy since we abuse it as some kind of "exact solution"
  # although this doesn't really exist...
  initial_condition::InitialCondition

  boundary_conditions::BoundaryConditions
  source_terms::SourceTerms
  solver::Solver
  cache::Cache
  performance_counter::PerformanceCounter

  function SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache}(
      mesh::Mesh, equations::Equations,
      initial_condition::InitialCondition, boundary_conditions::BoundaryConditions,
      source_terms::SourceTerms,
      solver::Solver, cache::Cache) where {Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache}
    @assert ndims(mesh) == ndims(equations)

    performance_counter = PerformanceCounter()

    new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, performance_counter)
  end
end

"""
    SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                 source_terms=nothing,
                                 boundary_conditions=boundary_condition_periodic,
                                 RealT=real(solver),
                                 uEltype=RealT,
                                 initial_cache=NamedTuple())

Construct a semidiscretization of a hyperbolic PDE.
"""
function SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                      source_terms=nothing,
                                      boundary_conditions=boundary_condition_periodic,
                                      # `RealT` is used as real type for node locations etc.
                                      # while `uEltype` is used as element type of solutions etc.
                                      RealT=real(solver), uEltype=RealT,
                                      initial_cache=NamedTuple())

  cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...)
  _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache)

  SemidiscretizationHyperbolic{typeof(mesh), typeof(equations), typeof(initial_condition), typeof(_boundary_conditions), typeof(source_terms), typeof(solver), typeof(cache)}(
    mesh, equations, initial_condition, _boundary_conditions, source_terms, solver, cache)
end

# Create a new semidiscretization but change some parameters compared to the input.
# `Base.similar` follows a related concept but would require us to `copy` the `mesh`,
# which would impact the performance. Instead, `SciMLBase.remake` has exactly the
# semantics we want to use here. In particular, it allows us to re-use mutable parts,
# e.g. `remake(semi).mesh === semi.mesh`.
function remake(semi::SemidiscretizationHyperbolic; uEltype=real(semi.solver),
                                                    mesh=semi.mesh,
                                                    equations=semi.equations,
                                                    initial_condition=semi.initial_condition,
                                                    solver=semi.solver,
                                                    source_terms=semi.source_terms,
                                                    boundary_conditions=semi.boundary_conditions
                                                    )
  # TODO: Which parts do we want to `remake`? At least the solver needs some
  #       special care if shock-capturing volume integrals are used (because of
  #       the indicators and their own caches...).
  SemidiscretizationHyperbolic(
    mesh,  equations, initial_condition, solver; source_terms, boundary_conditions, uEltype)
end

# general fallback
digest_boundary_conditions(boundary_conditions, mesh, solver, cache) = boundary_conditions

# general fallback
digest_boundary_conditions(boundary_conditions::BoundaryConditionPeriodic,
                           mesh, solver, cache) = boundary_conditions

# resolve ambiguities with definitions below
digest_boundary_conditions(boundary_conditions::BoundaryConditionPeriodic,
                           mesh::Union{TreeMesh{1}, StructuredMesh{1}}, solver, cache) = boundary_conditions

digest_boundary_conditions(boundary_conditions::BoundaryConditionPeriodic,
                           mesh::Union{TreeMesh{2}, StructuredMesh{2}}, solver, cache) = boundary_conditions

digest_boundary_conditions(boundary_conditions::BoundaryConditionPeriodic,
                           mesh::Union{TreeMesh{3}, StructuredMesh{3}}, solver, cache) = boundary_conditions

# allow passing a single BC that get converted into a tuple of BCs
# on (mapped) hypercube domains
function digest_boundary_conditions(boundary_conditions,
                                    mesh::Union{TreeMesh{1}, StructuredMesh{1}}, solver, cache)
  (; x_neg=boundary_conditions, x_pos=boundary_conditions)
end

function digest_boundary_conditions(boundary_conditions,
                                    mesh::Union{TreeMesh{2}, StructuredMesh{2}}, solver, cache)
  (; x_neg=boundary_conditions, x_pos=boundary_conditions,
     y_neg=boundary_conditions, y_pos=boundary_conditions)
end

function digest_boundary_conditions(boundary_conditions,
                                    mesh::Union{TreeMesh{3}, StructuredMesh{3}}, solver, cache)
  (; x_neg=boundary_conditions, x_pos=boundary_conditions,
     y_neg=boundary_conditions, y_pos=boundary_conditions,
     z_neg=boundary_conditions, z_pos=boundary_conditions)
end

# allow passing a tuple of BCs that get converted into a named tuple to make it
# self-documenting on (mapped) hypercube domains
function digest_boundary_conditions(boundary_conditions::NTuple{2, Any},
                                    mesh::Union{TreeMesh{1}, StructuredMesh{1}}, solver, cache)
  (; x_neg=boundary_conditions[1], x_pos=boundary_conditions[2])
end

function digest_boundary_conditions(boundary_conditions::NTuple{4, Any},
                                    mesh::Union{TreeMesh{2}, StructuredMesh{2}}, solver, cache)
  (; x_neg=boundary_conditions[1], x_pos=boundary_conditions[2],
     y_neg=boundary_conditions[3], y_pos=boundary_conditions[4])
end

function digest_boundary_conditions(boundary_conditions::NTuple{6, Any},
                                    mesh::Union{TreeMesh{3}, StructuredMesh{3}}, solver, cache)
  (; x_neg=boundary_conditions[1], x_pos=boundary_conditions[2],
     y_neg=boundary_conditions[3], y_pos=boundary_conditions[4],
     z_neg=boundary_conditions[5], z_pos=boundary_conditions[6])
end

# allow passing named tuples of BCs constructed in an arbitrary order
# on (mapped) hypercube domains
function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes},
                                    mesh::Union{TreeMesh{1}, StructuredMesh{1}}, solver, cache) where {Keys, ValueTypes<:NTuple{2,Any}}
  @unpack x_neg, x_pos = boundary_conditions
  (; x_neg, x_pos)
end

function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes},
                                    mesh::Union{TreeMesh{2}, StructuredMesh{2}}, solver, cache) where {Keys, ValueTypes<:NTuple{4,Any}}
  @unpack x_neg, x_pos, y_neg, y_pos = boundary_conditions
  (; x_neg, x_pos, y_neg, y_pos)
end

function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes},
                                    mesh::Union{TreeMesh{3}, StructuredMesh{3}}, solver, cache) where {Keys, ValueTypes<:NTuple{6,Any}}
  @unpack x_neg, x_pos, y_neg, y_pos, z_neg, z_pos = boundary_conditions
  (; x_neg, x_pos, y_neg, y_pos, z_neg, z_pos)
end

# sort the boundary conditions from a dictionary and into tuples
function digest_boundary_conditions(boundary_conditions::Dict, mesh, solver, cache)
  UnstructuredSortedBoundaryTypes(boundary_conditions, cache)
end

function digest_boundary_conditions(boundary_conditions::AbstractArray, mesh, solver, cache)
  throw(ArgumentError("Please use a (named) tuple instead of an (abstract) array to supply multiple boundary conditions (to improve performance)."))
end

function Base.show(io::IO, semi::SemidiscretizationHyperbolic)
  @nospecialize semi # reduce precompilation time

  print(io, "SemidiscretizationHyperbolic(")
  print(io,       semi.mesh)
  print(io, ", ", semi.equations)
  print(io, ", ", semi.initial_condition)
  print(io, ", ", semi.boundary_conditions)
  print(io, ", ", semi.source_terms)
  print(io, ", ", semi.solver)
  print(io, ", cache(")
  for (idx,key) in enumerate(keys(semi.cache))
    idx > 1 && print(io, " ")
    print(io, key)
  end
  print(io, "))")
end

function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperbolic)
  @nospecialize semi # reduce precompilation time

  if get(io, :compact, false)
    show(io, semi)
  else
    summary_header(io, "SemidiscretizationHyperbolic")
    summary_line(io, "#spatial dimensions", ndims(semi.equations))
    summary_line(io, "mesh", semi.mesh)
    summary_line(io, "equations", semi.equations |> typeof |> nameof)
    summary_line(io, "initial condition", semi.initial_condition)

    print_boundary_conditions(io, semi)

    summary_line(io, "source terms", semi.source_terms)
    summary_line(io, "solver", semi.solver |> typeof |> nameof)
    summary_line(io, "total #DOFs", ndofs(semi))
    summary_footer(io)
  end
end

# type alias for dispatch in printing of boundary conditions
const SemiHypMeshBCSolver{Mesh, BoundaryConditions, Solver} =
      SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition, BoundaryConditions,
                                   SourceTerms, Solver} where {Equations, InitialCondition, SourceTerms}

# generic fallback: print the type of semi.boundary_condition.
print_boundary_conditions(io, semi::SemiHypMeshBCSolver) = summary_line(io, "boundary conditions", typeof(semi.boundary_conditions))

function print_boundary_conditions(io, semi::SemiHypMeshBCSolver{<:Any, <:UnstructuredSortedBoundaryTypes})
  @unpack boundary_conditions = semi
  @unpack boundary_dictionary = boundary_conditions
  summary_line(io, "boundary conditions", length(boundary_dictionary))
  for (boundary_name, boundary_condition) in boundary_dictionary
    summary_line(increment_indent(io), boundary_name, typeof(boundary_condition))
  end
end

function print_boundary_conditions(io, semi::SemiHypMeshBCSolver{<:Any, <:NamedTuple})
  @unpack boundary_conditions = semi
  summary_line(io, "boundary conditions", length(boundary_conditions))
  bc_names = keys(boundary_conditions)
  for (i, bc_name) in enumerate(bc_names)
    summary_line(increment_indent(io), String(bc_name), typeof(boundary_conditions[i]))
  end
end

function print_boundary_conditions(io, semi::SemiHypMeshBCSolver{<:Union{TreeMesh, StructuredMesh}, <:Union{Tuple,NamedTuple,AbstractArray}})
  summary_line(io, "boundary conditions", 2*ndims(semi))
  bcs = semi.boundary_conditions

  summary_line(increment_indent(io), "negative x", bcs[1])
  summary_line(increment_indent(io), "positive x", bcs[2])
  if ndims(semi) > 1
    summary_line(increment_indent(io), "negative y", bcs[3])
    summary_line(increment_indent(io), "positive y", bcs[4])
  end
  if ndims(semi) > 2
    summary_line(increment_indent(io), "negative z", bcs[5])
    summary_line(increment_indent(io), "positive z", bcs[6])
  end
end

@inline Base.ndims(semi::SemidiscretizationHyperbolic) = ndims(semi.mesh)

@inline nvariables(semi::SemidiscretizationHyperbolic) = nvariables(semi.equations)

@inline Base.real(semi::SemidiscretizationHyperbolic) = real(semi.solver)

@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolic)
  @unpack mesh, equations, solver, cache = semi
  return mesh, equations, solver, cache
end

function calc_error_norms(func, u_ode, t, analyzer, semi::SemidiscretizationHyperbolic, cache_analysis)
  @unpack mesh, equations, initial_condition, solver, cache = semi
  u = wrap_array(u_ode, mesh, equations, solver, cache)

  calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis)
end

function compute_coefficients(t, semi::SemidiscretizationHyperbolic)
  # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl`
  compute_coefficients(semi.initial_condition, t, semi)
end

function compute_coefficients!(u_ode, t, semi::SemidiscretizationHyperbolic)
  compute_coefficients!(u_ode, semi.initial_condition, t, semi)
end

function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t)
  @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi

  u  = wrap_array(u_ode,  mesh, equations, solver, cache)
  du = wrap_array(du_ode, mesh, equations, solver, cache)

  # TODO: Taal decide, do we need to pass the mesh?
  time_start = time_ns()
  @trixi_timeit timer() "rhs!" rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache)
  runtime = time_ns() - time_start
  put!(semi.performance_counter, runtime)

  return nothing
end

end # @muladd


text:
struct TestSourceWithoutLength
end

function Base.eltype(iter::TestSourceWithoutLength)
    return NamedTuple{(:a, :b), Tuple{Int, Float64}}
end

Base.IteratorSize(::Type{T}) where {T <: TestSourceWithoutLength} = Base.SizeUnknown()

function Base.iterate(iter::TestSourceWithoutLength, state=1)
    if state==1
        return (a=1, b=1.), 2
    elseif state==2
        return (a=2, b=2.), 3
    else
        return nothing
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 88419b2e-5c5a-11eb-2885-fb367995002b
begin
	using DrWatson

end

# ╔═╡ 7dbb3cb4-5c5a-11eb-3191-db2a1e5dddd5
begin
	@quickactivate

	using PlutoUI
	using Serialization

	const pd = projectdir()
end

# ╔═╡ 75678b5a-5c5a-11eb-195e-2f7e93faacb1
dice_faces = [ load("$pd/Base/ui/dice/$fn") for fn in
	readdir("$pd/Base/ui/dice") if occursin("gif", fn) ];

# ╔═╡ 94e33342-5c5a-11eb-348a-01e6747337f8
face_preview = [ LocalResource("$pd/Base/ui/dice/$fn") for fn in
	readdir("$pd/Base/ui/dice") if occursin("gif", fn) ]

# ╔═╡ e807667e-5c5a-11eb-2175-7b84dffeb740
md"""
Save data to jls file? $(@bind save_data CheckBox())
"""

# ╔═╡ 80fa2680-5c5a-11eb-3e92-05a4be927f4e
if save_data
	serialize("$pd/Base/ui/dice/dice_faces.jls", dice_faces)
end

# ╔═╡ Cell order:
# ╠═88419b2e-5c5a-11eb-2885-fb367995002b
# ╠═7dbb3cb4-5c5a-11eb-3191-db2a1e5dddd5
# ╠═75678b5a-5c5a-11eb-195e-2f7e93faacb1
# ╠═94e33342-5c5a-11eb-348a-01e6747337f8
# ╟─e807667e-5c5a-11eb-2175-7b84dffeb740
# ╠═80fa2680-5c5a-11eb-3e92-05a4be927f4e