text:
"""Module Sario
"""

__precompile__(true)

module Sario

include("./demrsc.jl")

export load, save, DemRsc, load_stack, load_hdf5_stack, load_mask, load_geolist_from_h5,
       load_intlist_from_h5, load_dem_from_h5, get_file_ext, find_files, save_hdf5_stack

# TODO: python transfers:
# const LOAD_IN_PYTHON = vcat(IMAGE_EXTS, [".geojson", ".npy"])
# GeoJSON.geo2dict(GeoJSON.parsefile("../../bbox2.geojson"))["coordinates"]

import Glob
import JSON
using HDF5
using Dates
using Parameters

const DATE_FMT = "yyyymmdd"
const SENTINEL_EXTS = [".geo", ".cc", ".int", ".amp", ".unw", ".unwflat"]
const COMPLEX_EXTS = [".int", ".slc", ".geo", ".cc", ".unw", ".unwflat", ".mlc", ".grd"]
const REAL_EXTS = [".amp", ".cor", ".mlc", ".grd"]
const ELEVATION_EXTS = [".dem", ".hgt"]
# These file types are not simple complex matrices: see load_stacked_img for detail
# .unwflat are same as .unw, but with a linear ramp removed
const STACKED_EXTS = [".cc", ".unw", ".unwflat"]
const BOOL_EXTS = [".mask"]

const IMAGE_EXTS = [".png", ".tif", ".tiff", ".jpg"]
const NOT_IMPLEMENTED = vcat(IMAGE_EXTS, [".geojson", ".npy"])

# For HDF5 saving
const DEM_RSC_DSET = "dem_rsc"
const STACK_DSET = "stack"
const GEOLIST_DSET = "geo_dates"
const INTLIST_DSET = "int_dates"

RSC_KEY_TYPES = [
    ("width", Int),
    ("file_length", Int),
    ("x_first", Float64),
    ("y_first", Float64),
    ("x_step", Float64),
    ("y_step", Float64),
    ("x_unit", String),
    ("y_unit", String),
    ("z_offset", Int),
    ("z_scale", Int),
    ("projection", String),
]
RSC_KEYS = [tup[1] for tup in RSC_KEY_TYPES]

find_files(ext, directory = ".") = sort(Glob.glob("*" * ext, directory))

"""Extracts the file extension, including the "." (e.g.: .slc)"""
get_file_ext(filename::AbstractString) = splitext(filename)[end]

_is_h5(filename) = get_file_ext(filename) in (".h5", ".hdf5")

import Base.size
function Base.size(h5file::String, dset::String)
    h5open(h5file) do f
        return size(f[dset])
    end
end

"""
    load(filename::AbstractString;
         rsc_file::Union{AbstractString, Nothing}=nothing,
         looks::Tuple{Int, Int}=(1, 1), do_permute=true,
         dset_name::AbstractString="", return_amp::Bool=false)

Main entry point to loading files.

Examines file type of `filename` and runs appropriate load function.

"""
function load(filename::AbstractString;
              rsc_file::Union{AbstractString,Nothing} = nothing,
              looks::Tuple{Int,Int} = (1, 1), do_permute = true,
              dset_name::AbstractString = "", return_amp::Bool = false)
    ext = get_file_ext(filename)

    # For now, just pass through unimplemented extensions to Python
    if ext in NOT_IMPLEMENTED
        error("$ext is not yet implemented")
    elseif ext == ".rsc"
        return take_looks(_get_rsc_data(filename, filename), looks...)
    elseif _is_h5(filename)
        return take_looks(_load_hdf5(filename, dset_name; do_permute = do_permute), looks...)
    end

    # Sentinel files should have .rsc file: check for dem.rsc, or elevation.rsc
    data_type, demrsc, num_rows, num_cols = _file_info(filename, rsc_file)
    # If the filesize doesn't match the .rsc, error
    _check_filesize(filename, data_type, (num_rows, num_cols))

    if ext in ELEVATION_EXTS
        return take_looks(load_elevation(filename, demrsc, do_permute = do_permute), looks...)
    elseif ext in STACKED_EXTS
        return take_looks(load_stacked_img(filename, demrsc, do_permute = do_permute, return_amp = return_amp),
                          looks...)
    elseif ext in BOOL_EXTS
        return take_looks(load_bool(filename, demrsc, do_permute = do_permute), looks...)
    else
        return take_looks(load_complex(filename, demrsc, do_permute = do_permute), looks...)
    # having demrsc implies that this is not a UAVSAR file, so is complex
    # TODO: haven"t transferred over UAVSAR functions, so no load_real yet
    end

end

function _check_filesize(fname, data_type, dsize)
    fsize = Int(filesize(fname) / sizeof(data_type))
    # The stacked files have 2x the number of pixels of Floats- needs to be scaled to match dem.rsc
    fsize = get_file_ext(fname) in STACKED_EXTS ? fsize / 2 : fsize
    @assert prod(dsize) == fsize ".rsc size $dsize does not match $fname size $fsize"
end

"""Load one element of a file on disk (avoid reading in all of huge file"""
# TODO: Load a chunk of a file now?
function load(filename::AbstractString, row_col::Tuple{Int,Int}; rsc_file::Union{AbstractString,Nothing} = nothing)
    data_type, demrsc, num_rows, num_cols = _file_info(filename, rsc_file)

    row, col = row_col
    _check_bounds(row, col, num_rows, num_cols)

    seek_pos = _get_seek_position(row, col, num_cols, data_type)

    open(filename) do f
        seek(f, seek_pos)
        # This read syntax loads single `data_type`
        return read(f, data_type)
    end
end

"""For single element reading in binary files, seek to the right row, col"""
_get_seek_position(row, col, num_cols, data_type) = sizeof(data_type) * ((col - 1) + (num_cols * (row - 1)) )

function _check_bounds(row, col, num_rows, num_cols)
    if row < 1 || col < 1 || row > num_rows || col > num_cols
        throw(DomainError((row, col), " out of bounds for $filename of size ($num_rows, $num_cols)"))
    end
end

function _file_info(filename, rsc_file)
    @assert isfile(filename) "$filename does not exist."

    data_type = _get_data_type(filename)
    demrsc = _get_rsc_data(filename, rsc_file)
    @unpack rows, cols = demrsc
    return data_type, demrsc, rows, cols
end

# NOTE: the subset will only work for sequential data (complex), not the stacked filetypes
"""Load subset of a file on disk using range"""
RangeTuple = Tuple{T,T} where {T <: OrdinalRange}
function load(filename::AbstractString, idxs::RangeTuple;
              rsc_file::Union{AbstractString,Nothing} = nothing,  do_permute = true)
    data_type, demrsc, num_rows, num_cols = _file_info(filename, rsc_file)

    # Note: reversing these since julia uses column-major
    cols, rows = idxs
    _check_bounds(rows.start, cols.start, num_rows, num_cols)

    outrows = rows.stop  - rows.start + 1
    outcols = cols.stop  - cols.start + 1
    out = zeros(data_type, (outrows, outcols))
    buf = zeros(data_type, outcols)

    open(filename) do f
        for (idx, r) in enumerate(rows)
            seek_pos = _get_seek_position(r, cols.start, num_cols, data_type)
            seek(f, seek_pos)
            # This read syntax loads single `data_type`
            # read!(out[r, :], f, data_type)
            read!(f, buf)
            out[:, idx] .= buf
        end
    end
    return do_permute ? permutedims(out) : out
end

function _load_hdf5(h5file::AbstractString, dset_name::AbstractString = ""; do_permute = true)
    h5open(h5file, "r") do f
        if isempty(dset_name)
            if length(names(f)) == 1 && typeof(f[names(f)[1]]) == HDF5Dataset
                dset_name = names(f)[1]
                println("Loading \"$dset_name\" from $h5file")
            else
                error("More than 1 dset exists in $h5file, provide name to `load`")
            end
        end
        order = ndims(f[dset_name]) == 3 ? (2, 1, 3) : (2, 1)  # Handle 3D permuting
        return do_permute ? permutedims(read(f, dset_name), order) : read(f, dset_name)
        k
    end
end

"""Handle getting the .rsc data from either an image filename, or .rsc filename"""
function _get_rsc_data(filename, rsc_file::Union{AbstractString,Nothing})
    ext = get_file_ext(filename)

    demrsc = nothing
    if !isnothing(rsc_file)
        demrsc = load_dem_rsc(rsc_file)
    elseif ext in vcat(SENTINEL_EXTS, ELEVATION_EXTS, BOOL_EXTS)
        rsc_file = find_rsc_file(filename)
        demrsc = load_dem_rsc(rsc_file)
    end

    return demrsc
end

function find_rsc_file(filename = nothing; directory = nothing, verbose = false)
    if !isnothing(filename)
        dir_parts = splitpath(filename)[1:end - 1]
        directory = isempty(dir_parts) ? "." : joinpath(dir_parts...)
    end
    # Should be just elevation.dem.rsc (for .geo folder) or dem.rsc (for igrams)
    possible_rscs = find_files("*.rsc", directory)
    if verbose
        println("Searching $directory for rsc files")
        println("Possible rsc files:")
        println(possible_rscs)
    end
    if length(possible_rscs) < 1
        println("No .rsc file found in $directory")
        return nothing
    elseif length(possible_rscs) > 1
        pathlist = splitpath(filename)
        if length(pathlist) == 1
            filepath, fileonly = ".", pathlist[1]
        else
            filepath, fileonly = joinpath(pathlist[1:end - 1]...), pathlist[end]
        end
        rsc_nopaths = [splitpath(r)[end] for r in possible_rscs]
        if any(startswith(r, fileonly) for r in rsc_nopaths)
            possible_rscs = [joinpath(filepath, r) for r in rsc_nopaths if startswith(r, fileonly)]
        else
            error("$filename has multiple .rsc files in its directory: $possible_rscs")
        end
    end
    return abspath(possible_rscs[1])
end

"""Convert the text file into the DemRsc struct
Starts with a dict to gather all fiields, then unpacks using keywords args"""
function load_dem_rsc(filename, kwargs...)::DemRsc
    output_data = Dict{Symbol,Any}()

    for line in readlines(filename)
        for (field, valtype) in RSC_KEY_TYPES
            if startswith(line, uppercase(field))
                val = valtype <: AbstractString ? split(line)[2] : parse(valtype, split(line)[2])
                output_data[Symbol(lowercase(field))] = val
            end
        end
    end
    return DemRsc(;output_data...)
end

function _get_data_type(filename)
    ext = get_file_ext(filename)
    if ext in ELEVATION_EXTS
        return Int16
    elseif ext in STACKED_EXTS
        return Float32
    elseif ext in BOOL_EXTS
        return Bool
    else
        return ComplexF32
    end
end

"""Loads a digital elevation map from either .hgt file or .dem

.hgt is the NASA SRTM files given. Documentation on format here:
https://dds.cr.usgs.gov/srtm/version2_1/Documentation/SRTM_Topo.pdf
Key point: Big-endian 2 byte (16-bit) integers

.dem is format used by Zebker geo-coded and ROI-PAC SAR software
Only difference is data is stored little-endian (like other SAR data)

Note on both formats: gaps in coverage are given by INT_MIN -32768,
so either manually set data(data == np.min(data)) = 0,
or something like
    data = clamp(data, -10000, Inf)
"""
function load_elevation(filename, demrsc::DemRsc; do_permute = true)
    rows, cols = (demrsc.file_length, demrsc.width)
    data = Array{Int16,2}(undef, (cols, rows))
    read!(filename, data)

    # # TODO: Verify that the min real value will be above -1000
    # min_valid = -10000
    # # Set NaN values to 0
    # @. data[data < min_valid] = 0
    return do_permute ? permutedims(data) : data
end

function load_complex(filename::AbstractString, demrsc::DemRsc; do_permute = true)
    return _load_bin_matrix(filename, demrsc, ComplexF32, do_permute)
end

function load_bool(filename::AbstractString, demrsc::DemRsc; do_permute = true)
    return _load_bin_matrix(filename, demrsc, Bool, do_permute)
end

function _load_bin_matrix(filename, demrsc::DemRsc, dtype, do_permute::Bool)
    @unpack rows, cols = demrsc
    # Note: must be saved in c/row-major order, so loading needs a transpose
    out = Array{dtype,2}(undef, (cols, rows))

    read!(filename, out)
    return do_permute ? permutedims(out) : out
end

"""load_stacked_img is for the weird "ALT_LINE_DATA" formatted images:
https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/snaphu_man1.html#FILE%20FORMATS

Format is two stacked matrices:
[[first], [second]] where the first "cols" number of floats
are the first matrix, next "cols" are second, etc.

For .unw height files, the first is amplitude, second is phase (unwrapped)
For .cc correlation files, first is amp, second is correlation (0 to 1)
"""
function load_stacked_img(filename::AbstractString, demrsc::DemRsc; do_permute = true, return_amp = false)
    @unpack rows, cols = demrsc
    # Note: must be saved in c/row-major order, so loading needs a transpose
    # out_left usually is amplitude data
    # Usually we are interested in out_right
    #
    # First make container for all of data
    out = Array{Float32,2}(undef, (2cols, rows))
    read!(filename, out)

    out_amp = @view out[1:cols, :]
    out_data = @view out[cols + 1:end, :]

    if return_amp
        return do_permute ? permutedims(cat(out_amp, out_data, dims = 3), (2, 1, 3)) : cat(out_amp, out_data, dims = 3)
    else
        return do_permute ? permutedims(out_data) : out_data
    end
end

function load_geolist_from_h5(h5file::AbstractString; geolist_dset = GEOLIST_DSET)
    geo_strings = h5read(h5file, geolist_dset)
    return parse_geolist_strings(geo_strings)
end

function load_geolist_from_h5(h5file::AbstractString, dset::AbstractString; geolist_dset = GEOLIST_DSET)
    geo_strings = h5readattr(h5file, dset)[geolist_dset]
    return parse_geolist_strings(geo_strings)
end

function load_intlist_from_h5(h5file)
    # Note transpose, since it's stored as a N x 2 array
    # (which is flipped to 2 x N for julia
    int_strings = h5read(h5file, INTLIST_DSET)
    return parse_intlist_strings([Tuple(int_strings[:, ii]) for ii in 1:size(int_strings, 2)])
end

function load_intlist_from_h5(h5file, dset::AbstractString)
    int_strings = h5readattr(h5file, dset)[INTLIST_DSET]
    return parse_intlist_strings([Tuple(int_strings[:, ii]) for ii in 1:size(int_strings, 2)])
end

_parse(datestr) = Date(datestr, DATE_FMT)

function parse_intlist_strings(date_pairs::AbstractArray{<:AbstractString,1})
    # For passing list of filenames ["20150101_20160101.int", ...]
    # `collect` used to make into an array of chars for strip
    date_pairs = [split(strip(d, collect(".int")), "_")[1:2] for d in date_pairs]
    date_tups = [Tuple(d) for d in date_pairs]
    return parse_intlist_strings(date_tups)
end

parse_intlist_strings(date_pairs) = [_parse.(pair) for pair in date_pairs]

parse_geolist_strings(geolist_str::AbstractArray{String}) = _parse.(geolist_str)

import Base.string
string(arr::AbstractArray{Date}, fmt = "yyyymmdd") = Dates.format.(arr, fmt)
string(arr::AbstractArray{Tuple{Date,Date}}, fmt = "yyyymmdd") = [Dates.format.(tup, fmt) for tup in arr]

"""
    save_geolist_to_h5(h5file::String, geolist::AbstractArray{Date}; overwrite=false)

Save the geolist as a list of strings to an dataset `h5file`"""
function save_geolist_to_h5(h5file::String, geolist::AbstractArray{Date}; overwrite = false)
    !check_dset(h5file, GEOLIST_DSET, overwrite) && return
    h5write(h5file, GEOLIST_DSET, string(geolist))
end

"""In this version, save the geolist to an attribute of `object` already in `h5file`"""
function save_geolist_to_h5(h5file::String, object::String, geolist::AbstractArray{Date}; overwrite = false)
    h5open(h5file, "cw") do f
        obj = f[object]
        # Delete if exists
        # TODO: wrap this attr one into check_dset
        overwrite && exists(attrs(obj), GEOLIST_DSET) && a_delete(obj, GEOLIST_DSET)

        attrs(obj)[GEOLIST_DSET] = string(geolist)
    end
end

"""Save the geolist as a list of strings to an dataset `h5file`"""
function save_intlist_to_h5(h5file::String, intlist::AbstractArray{Tuple{Date,Date}}; overwrite = false)
    !check_dset(h5file, INTLIST_DSET, overwrite) && return
    h5write(h5file, INTLIST_DSET, _intlist_to_strings(intlist))
end
#
_intlist_to_strings(intlist::AbstractArray{Tuple{Date,Date}})::AbstractArray{String,2} = hcat([collect(tup) for tup in string(intlist)]...)

# In HDF5, DemRscs are stored as JSON dicts
load_dem_from_h5(h5file, dset = DEM_RSC_DSET) = DemRsc(JSON.parse(h5read(h5file, dset)))

function save_dem_to_h5(h5file, demrsc::DemRsc, dset_name = DEM_RSC_DSET; overwrite = true)
    !check_dset(h5file, dset_name, overwrite) && return
    h5write(h5file, dset_name, JSON.json(demrsc))
end
# TODO:
# save_dem_to_h5(h5file, demrsc::DemRsc, dset_name=DEM_RSC_DSET; overwrite=true)

function save_hdf5_stack(h5file::AbstractString, dset_name::AbstractString, stack; overwrite::Bool = false, do_permute = true)
    # TODO: is there a way to combine this with normal saving of files?
    mode = overwrite ? "w" : "cw"  # Append mode is 'cw'  (why different than "a"?)
    h5open(h5file, mode) do f
        do_permute ? write(f, dset_name, permutedims(stack, (2, 1, 3))) :
                     write(f, dset_name, stack)
    end
end

    """Wrapper around h5read to account for the transpose
necessary when reading from python-written stacks
"""
function load_hdf5_stack(h5file::AbstractString, dset_name::AbstractString; do_permute = true)
    h5open(h5file) do f
        return do_permute ? permutedims(read(f[dset_name]), (2, 1, 3)) : read(f[dset_name])
    end
end

# If loading only certain layers, don't read all into memory
function load_hdf5_stack(h5file::AbstractString, dset_name::AbstractString, valid_layer_idxs; do_permute = true)
    nrows, ncols, _ = size(h5file, dset_name)
    h5open(h5file) do f
        dset = f[dset_name]
        out = Array{eltype(dset),ndims(dset)}(undef, (nrows, ncols, length(valid_layer_idxs)))
        for (tot_idx, v_idx) in enumerate(valid_layer_idxs)
            out[:, :, tot_idx] = dset[:, :, v_idx]
        end
        return do_permute ? permutedims(out, (2, 1, 3)) : out
    end
end

"""Sum the 3rd dim (layers) of a stack without loading all into memory"""
function sum_hdf5_stack(h5file::AbstractString, dset_name::AbstractString, valid_layer_idxs)
    rows, cols, _ = size(h5file, dset_name)
    # For summing, using int instead of bool
    dtype = eltype(h5file, dset_name) <: Bool ? Int : eltype(h5file, dset_name)
    out = zeros(dtype, (rows, cols))

    h5open(h5file) do f
        d = f[dset_name]
        for ii in valid_layer_idxs
            out .+= @view d[:, :, ii][:, :, 1]
        end
    end
    return out
    end

"""Get the composite mask from the stack, true only where ALL pixels are masked"""
function load_mask(geolist::AbstractArray{Date,1}; do_permute::Bool = true, fname = "masks.h5", dset_name = "geo")
    geolist_full = load_geolist_from_h5(fname)
    idxs = indexin(geolist, geolist_full)
    out = convert(Array{Bool}, sum_hdf5_stack(fname, dset_name, idxs) .> 0)
    return do_permute ? permutedims(out) : out
end

# TODO: this sucks and doesnt work well... if one day masks out half the map,
# the geo_sum will be 0 on half, 1 on half. if we ignore that date... this still
# says "it equals the maximum of the sum, so it's masked
"""Get the composite mask from the stack, true only where ALL pixels are masked"""
function load_mask(;do_permute::Bool = true, fname = "masks.h5", dset = "geo_sum")
    mask = h5read(fname, dset)
    mval = max(1, maximum(mask))  # Make sure we dont mask all 0s
    return do_permute ? permutedims(mask .== mval) : mask .== mval
end

"""Returns false if the dataset exists and overwrite is false

If overwrite is set to true, will delete the dataset to make
sure a new one can be created"""
function check_dset(h5file, dset_name, overwrite)
    h5open(h5file, "cw") do f
        if dset_name in names(f)
            println("$dset_name already exists in $h5file")
            if overwrite
                println("Overwrite true: Deleting.")
                o_delete(f, dset_name)
            else
                println("Skipping $dset_name: not deleting")
                return false
            end
        end

        return true
    end
end

# TODO: probably a better way to do this.. but can't figure out how to
# without allocating so many arrays that it's slow as Python
function load_stack(; file_list::Union{AbstractArray{AbstractString},Nothing} = nothing,
                    directory::Union{AbstractString,Nothing} = nothing,
                    file_ext::Union{AbstractString,Nothing} = nothing)
    if isnothing(file_list)
        file_list = find_files(file_ext, directory)
    end

    # demrsc = load(find_rsc_file(basepath=directory))
    test_arr = load(file_list[1])
    rows, cols = size(test_arr)
    T = eltype(test_arr)

    stack, buffer = _return_array(T, rows, cols, length(file_list))
    for (idx, f) in enumerate(file_list)
         read!(f, buffer)
            stack[:, :, idx] = buffer
    end

    return _permute(stack, cols)
end

# Using multiple dispatch to avoid if statements for 2 types of images
function _return_array(T::Union{Type{ComplexF32},Type{Bool}}, rows::Int, cols::Int, file_len::Int)
    stack = Array{T,3}(undef, (cols, rows, file_len))
    buffer = Array{T,2}(undef, (cols, rows))
    stack, buffer
    end

function _return_array(T::Type{Float32}, rows::Int, cols::Int, file_len::Int)
    stack = Array{T,3}(undef, (2cols, rows, file_len))
    buffer = Array{T,2}(undef, (2cols, rows))
    stack, buffer
    end

# For normal .int, .geo complex/ masks, just transpose stack
function _permute(stack::AbstractArray{T,3}, cols::Int) where {T <: Number}
    return permutedims(stack, (2, 1, 3))
end

# For normal weird stacked types, pick just the right half
function _permute(stack::AbstractArray{Float32,3}, cols::Int)
    return permutedims(stack[cols + 1:end, :, :], (2, 1, 3))
end

# function _load_stack_complex(file_list::AbstractArray{AbstractString}, rows::Int, cols::Int)
#     stack = Array{ComplexF32, 3}(undef, (cols, rows, length(file_list)))
#     buffer = Array{ComplexF32, 2}(undef, (cols, rows))
#
#     for (idx, f) in enumerate(file_list)
#          read!(f, buffer)
#             stack[:, :, idx] = buffer
#     end
#
#     return permutedims(stack, (2, 1, 3))
# end
#
#
# function _load_stack_stacked(file_list::AbstractArray{AbstractString}, rows::Int, cols::Int)
#     stack = Array{Float32, 3}(undef, (2cols, rows, length(file_list)))
#     buffer = Array{Float32, 2}(undef, (2cols, rows))
#
#     for (idx, f) in enumerate(file_list)
#          read!(f, buffer)
#             stack[:, :, idx] = buffer
#     end
#     return permutedims(stack[cols+1:end, :, :], (2, 1, 3))
# end

"""
    save(filename::AbstractString, array; do_permute=true, kwargs...)

Saves a 2D/3D array as `filename`.
If the file was loaded with `do_permute=false`, it should also be
saved as `do_permute=false`

For stacked files, the expected format of `array` is a 3D stack with
`array[:,:,1]` being the amplitude image, `array[:,:,2]` being the data.
If only a 2D image is passed, an dummy amplitude image of all 1s is made.
"""
function save(filename::AbstractString, array; do_permute = true, kwargs...)
    ext = get_file_ext(filename)

    if ext == ".rsc"
        write(filename, string(array))  # array here should be a DemRsc
    elseif ext in BOOL_EXTS
        tofile(filename, array, do_permute = do_permute)
    elseif (ext in vcat(COMPLEX_EXTS, REAL_EXTS, ELEVATION_EXTS)) && (!(ext in STACKED_EXTS))
        tofile(filename, _force_float32(array), do_permute = do_permute)
    elseif ext in STACKED_EXTS
        # ndims(array) != 3 && throw(DimensionMismatch("array must be 3D [amp; data] to save as $filename"))
        if ndims(array) == 3
            amp = view(array, :, :, 1)
            data = view(array, :, :, 2)
        else
            println("!!! Warning: saving $filename with 1s for amplitude")
            data = array
            amp = ones(size(array))
        end
        # Handle the permuting for stacked here outside of function
        out = do_permute ? transpose(hcat(amp, data)) : vcat(amp, data)
        tofile(filename, _force_float32(out), do_permute = false)
    else
        error("Filetype not implemented: $filename")
    end
end

function tofile(filename, array; do_permute = true)  # , overwrite=true)
    # mode = overwrite ? "w" : "a"
    mode = "w"
    open(filename, mode) do f
        do_permute ? write(f, transpose(array)) : write(f, array)
    end
end

_iscomplex(array) = eltype(array) <: Complex
# _isfloat(array) = eltype(array) <: AbstractFloat

function _force_float32(array)
    if eltype(array) <: Complex
        return ComplexF32.(array)
    elseif eltype(array) <: AbstractFloat
        return Float32.(array)
    else
        return array
    end
end

"""Downsample a matrix by summing blocks of (row_looks, col_looks)

Cuts off values if the size isn't divisible by num looks
size = floor(rows / row_looks, cols / col_looks)
"""
function take_looks(image::AbstractArray{T}, row_looks, col_looks) where {T <: Number}
    (row_looks == 1 && col_looks == 1) && return image

    if ndims(image) > 2
        return cat((take_looks(image[:, :, i], row_looks, col_looks)
                    for i in 1:size(image, 3))..., dims = 3)
    end

    nrows, ncols = size(image)
    nr = div(nrows, row_looks)
    nc = div(ncols, col_looks)

    if T == Bool
        # Can't sum into a Bool array
        outtype = Float32
    elseif T <: Integer
        # Incase we take many looks for Int16 overflow, and to avoid InexactError
        outtype = Float64
    else
        outtype = T
    end
    out = zeros(outtype, nr, nc)
    return take_looks!(out, image, row_looks, col_looks)
end

function take_looks!(out::AbstractArray{S}, image::AbstractArray{T}, row_looks, col_looks) where {S,T <: Number}
    # Check that output is compatible
    nrows, ncols = size(image)
    nr = div(nrows, row_looks)
    nc = div(ncols, col_looks)
    (nr, nc) != size(out) && error("""
        Size of `out` $(size(out)) doesnt match taking ($row_looks, $col_looks)
        looks on image size $(size(image))""")
    return take_looks!(out, image, row_looks, col_looks, nr, nc)
end

function take_looks!(out::AbstractArray{S}, image::AbstractArray{T}, row_looks, col_looks, nr, nc) where {S,T <: Number}
    @inbounds Threads.@threads for j = 1:nc
        @inbounds for i = 1:nr
            indx_i = 1 + (i - 1) * row_looks:i * row_looks
            indx_j = 1 + (j - 1) * col_looks:j * col_looks
            @views out[i, j] = sum(image[indx_i, indx_j])
        end
    end

    out ./= (row_looks * col_looks)
    if T == Bool
        return Bool.(out .> 0)
    elseif T <: Integer
        return T.(round.(out))
    else
        return T.(out)
    end
end

function take_looks(demrsc::DemRsc, row_looks, col_looks)
    nrows, ncols = size(demrsc)
    newr = div(nrows, row_looks)
    newc = div(ncols, col_looks)
    new_x_step = demrsc.x_step * col_looks
    new_y_step = demrsc.y_step * row_looks
    return DemRsc(demrsc, file_length = newr, rows = newr, y_step = new_y_step,
                  width = newc, cols = newc, x_step = new_x_step)
end

# If we pass anything else:
# take_looks(other, row_looks, col_looks) = other

# For future:
# cc_patch(a, b) = real(abs(sum(a .* conj.(b))) / sqrt(sum(a .* conj.(a)) * sum(b .* conj.(b))))

"""
    find_igrams(;directory::AbstractString=".", parse::Bool=true, filename::AbstractString="")

Reads the list of igrams to returns Array of Igrams

Args:
    directory (str): path to the igram directory
    parse (bool): output as parsed datetime tuples. False returns the filenames
    filename (string): name of a file with .geo filenames

Returns:
    tuple(date, date) of (early, late) dates for all igrams (if parse=True)
        if parse=False: returns list[str], filenames of the igrams
"""
function find_igrams(;directory::AbstractString = ".", parse::Bool = true, filename::AbstractString = "")
    igram_file_list = isempty(filename) ? find_files("*.int", directory) : readlines(filename)
    intlist_strings = [splitpath(fname)[end] for fname in igram_file_list]
    isempty(intlist_strings) && return intlist_strings
    return parse ? parse_intlist_strings(intlist_strings) : intlist_strings
end

"""Reads in the list of .geo files used, in time order

Can also pass a filename containing .geo files as lines.

Args:
    directory (str): path to the geolist file or directory
    parse (bool): output as parsed datetime tuples. False returns the filenames
    filename (string): name of a file with .geo filenames

Returns:
    list[date]: the parse dates of each .geo used, in date order

"""
# TODO: parsers.Sentinel...
function find_geos(;directory::AbstractString = ".", parse::Bool = true, filename::AbstractString = "")
    geo_file_list = isempty(filename) ? find_files("*.geo", directory) : readlines(filename)

    !parse && return geo_file_list

    # Stripped of path for parser
    geolist_strings = [splitpath(fname)[end] for fname in geo_file_list]

    return sort(parse_geolist_strings(_strip_geoname.(geolist_strings)))

    # return sorted(parse_geolist_strings([_strip_geoname(geo) for geo in geolist_strings]))
    # if re.match(r'S1[AB]_\d{8}\.geo', geolist_strings[0]):  # S1A_YYYYmmdd.geo
    #     return sorted([_parse(_strip_geoname(geo)) for geo in geolist_strings])
    # elseif re.match(r'\d{8}', geolist_strings[0]):  # YYYYmmdd , just a date string
    #     return sorted([_parse(geo) for geo in geolist_strings])
    # else  # Full sentinel product name
    #     return sorted([apertools.parsers.Sentinel(geo).start_time.date() for geo in geolist_strings])
    # end
end

"""Leaves just date from format S1A_YYYYmmdd.geo"""
_strip_geoname(name) = reduce(replace, ["S1A_" => "", "S1B_" => "", ".geo" => ""], init = name)

end # module


text:
module PlotlyJS

using Base64
using Reexport
@reexport using PlotlyBase
using JSON
using REPL, Pkg, Pkg.Artifacts, DelimitedFiles  # stdlib

# need to import some functions because methods are meta-generated
import PlotlyBase:
    restyle!, relayout!, update!, addtraces!, deletetraces!, movetraces!,
    redraw!, extendtraces!, prependtraces!, purge!, to_image, download_image,
    restyle, relayout, update, addtraces, deletetraces, movetraces, redraw,
    extendtraces, prependtraces, prep_kwargs, sizes, savefig, _tovec,
    react, react!, add_trace!

using WebIO
using JSExpr
using JSExpr: @var, @new
using Blink
using Pkg.Artifacts
using Requires

export plot, dataset, list_datasets, make_subplots

# globals for this package
const _pkg_root = dirname(dirname(@__FILE__))
const _js_path = joinpath(artifact"plotly-artifacts", "plotly.min.js")
const _js_cdn_path = "https://cdn.plot.ly/plotly-latest.min.js"
const _mathjax_cdn_path =
    "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_SVG"

struct PlotlyJSDisplay <: AbstractDisplay end

# include the rest of the core parts of the package
include("display.jl")
include("util.jl")
include("kaleido.jl")

make_subplots(;kwargs...) = plot(Layout(Subplots(;kwargs...)))

@doc (@doc Subplots) make_subplots

function docs()
    schema_path = joinpath(dirname(dirname(@__FILE__)), "deps", "schema.html")
    if !isfile(schema_path)
        msg = "schema docs not built. Run `Pkg.build(\"PlotlyJS\")` to generate"
        error(msg)
    end
    w = Blink.Window()
    wait(w.content)
    Blink.content!(w, "html", open(f -> read(f, String), schema_path), fade=false, async=false)
end

@enum RENDERERS BLINK IJULIA BROWSER DOCS

const DEFAULT_RENDERER = Ref(BLINK)

function set_default_renderer(s::RENDERERS)
    global DEFAULT_RENDERER
    DEFAULT_RENDERER[] = s
end

@inline get_renderer() = DEFAULT_RENDERER[]

list_datasets() = readdir(joinpath(artifact"plotly-artifacts", "datasets"))
function check_dataset_exists(name::String)
    ds = list_datasets()
    name_ext = Dict(name => strip(ext, '.') for (name, ext) in splitext.(ds))
    if !haskey(name_ext, name)
        error("Unknown dataset $name, known datasets are $(collect(keys(name_ext)))")
    end
    ds_path = joinpath(artifact"plotly-artifacts", "datasets", "$(name).$(name_ext[name])")
    return ds_path
end

function dataset(name::String)::Dict{String,Any}
    ds_path = check_dataset_exists(name)
    if endswith(ds_path, "csv")
        # if csv, use DelimitedFiles and convert to dict
        data = DelimitedFiles.readdlm(ds_path, ',')
        return Dict(zip(data[1, :], data[2:end, i] for i in 1:size(data, 2)))
    elseif endswith(ds_path, "json")
        # use json
        return JSON.parsefile(ds_path)
    end
    error("should not ever get here!!! Please file an issue")
end

function __init__()
    _build_log = joinpath(_pkg_root, "deps", "build.log")
    if isfile(_build_log) && occursin("Warning:", read(_build_log, String))
        @warn("Warnings were generated during the last build of PlotlyJS:  please check the build log at $_build_log")
    end

    @async _start_kaleido_process()

    if !isfile(_js_path)
        @info("plotly.js javascript libary not found -- downloading now")
        include(joinpath(_pkg_root, "deps", "build.jl"))
    end

    # set default renderer
    # First check env var
    env_val = get(ENV, "PLOTLY_RENDERER_JULIA", missing)
    if !ismissing(env_val)
        env_symbol = Symbol(uppercase(env_val))
        options = Dict(v => k for (k, v) in collect(Base.Enums.namemap(PlotlyJS.RENDERERS)))
        renderer_int = get(options, env_symbol, missing)
        if ismissing(renderer_int)
            @warn "Unknown value for env var `PLOTLY_RENDERER_JULIA` \"$(env_val)\", known options are $(string.(keys(options)))"
        else
            set_default_renderer(RENDERERS(renderer_int))
        end
    else
        # we have no env-var
        # check IJULIA
        isdefined(Main, :IJulia) && Main.IJulia.inited && set_default_renderer(IJULIA)
    end

    # set up display
    insert!(Base.Multimedia.displays, findlast(x -> x isa Base.TextDisplay || x isa REPL.REPLDisplay, Base.Multimedia.displays) + 1, PlotlyJSDisplay())

    atreplinit(i -> begin
        while PlotlyJSDisplay() in Base.Multimedia.displays
            popdisplay(PlotlyJSDisplay())
        end
        insert!(Base.Multimedia.displays, findlast(x -> x isa REPL.REPLDisplay, Base.Multimedia.displays) + 1, PlotlyJSDisplay())
    end)

    @require IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" begin

        function IJulia.display_dict(p::SyncPlot)
            Dict(
                "application/vnd.plotly.v1+json" => JSON.lower(p),
                "text/plain" => sprint(show, "text/plain", p),
                "text/html" => let
                    buf = IOBuffer()
                    show(buf, MIME("text/html"), p)
                    String(resize!(buf.data, buf.size))
                end
            )
        end
    end

    @require CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" begin
        function dataset(::Type{CSV.File}, name::String)
            ds_path = check_dataset_exists(name)
            if !endswith(ds_path, "csv")
                error("Can only construct CSV.File from a csv data source")
            end
            CSV.File(ds_path)
        end
        @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
            dataset(::Type{DataFrames.DataFrame}, name::String) = DataFrames.DataFrame(dataset(CSV.File, name))
        end
    end
end

end # module

ext: jl


text:
""" ReactionSystem constructor """
function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires unique parameters (i.e. SBML must have been imported with localParameter promotion in libSBML)
    rxs = mtk_reactions(model)
    u0map = get_u0map(model)
    parammap = get_paramap(model)
    defs = ModelingToolkit._merge(Dict(u0map), Dict(parammap))

    algrules, obsrules, raterules = get_rules(model)
    for o in obsrules
        defs[o.lhs] = substitute(o.rhs, defs)
    end
    constraints_sys = ODESystem(vcat(algrules, raterules, obsrules), Catalyst.DEFAULT_IV; name = gensym(:CONSTRAINTS))

    ReactionSystem(rxs, Catalyst.DEFAULT_IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML),
        constraints = constraints_sys, kwargs...)
end

""" ODESystem constructor """
function ModelingToolkit.ODESystem(model::SBML.Model; include_zero_odes = false, kwargs...)
    rs = ReactionSystem(model; kwargs...)
    convert(ODESystem, rs; include_zero_odes = include_zero_odes)
end

""" Check if conversion to ReactionSystem is possible """
function checksupport(filename::String)
    not_implemented = ["listOfConstraints", "listOfEvents"]
    sbml = open(filename) do file
        read(file, String)
    end
    for item in not_implemented
        occursin(item, sbml) && throw(ErrorException("SBML models with $item are not yet implemented."))
    end
    occursin("<sbml xmlns:fbc=", sbml) && throw(ErrorException("This model was designed for constrained-based optimisation. Please use COBREXA.jl instead of SBMLToolkit."))
end

""" Convert intensive to extensive mathematical expression """
function to_extensive_math!(model::SBML.Model)
    function conv(x::SBML.MathApply)
        SBML.MathApply(x.fn, SBML.Math[conv(x) for x in x.args])
    end
    function conv(x::SBML.MathIdent)
        x_new = x
        if x.id in keys(model.species)
            specie = model.species[x.id]
            if !specie.only_substance_units
                compartment = model.compartments[specie.compartment]
                if isnothing(compartment.size)
                    @warn "Specie $(x.id) hasOnlySubstanceUnits but its compartment $(compartment.name) has no size. Cannot auto-correct the rate laws $(x.id) is involved in. Please check manually."
                else
                    x_new = SBML.MathApply("/", SBML.Math[x,
                        SBML.MathVal(compartment.size)])
                    specie.only_substance_units = true
                end
            end
        end
        x_new
    end
    conv(x::SBML.MathVal) = x
    conv(x::SBML.MathLambda) =
        throw(DomainError(x, "can't translate lambdas to extensive units"))
    for reaction in values(model.reactions)
        reaction.kinetic_math = conv(reaction.kinetic_math)
    end
    model
end

""" Get dictonary to change types in kineticLaw """
function _get_substitutions(model)
    subsdict = Dict()
    for k in keys(model.species)
        push!(subsdict, Pair(create_var(k), create_var(k, Catalyst.DEFAULT_IV)))
    end
    for k in keys(model.parameters)
        push!(subsdict, Pair(create_var(k), create_param(k)))
    end
    for k in keys(model.compartments)
        push!(subsdict, Pair(create_var(k), create_param(k)))
    end
    subsdict
end

""" Convert SBML.Reaction to MTK.Reaction """
function mtk_reactions(model::SBML.Model)
    subsdict = _get_substitutions(model)
    rxs = []
    for reaction in values(model.reactions)
        extensive_math = SBML.extensive_kinetic_math(
            model, reaction.kinetic_math,
            handle_empty_compartment_size = _ -> 1.0)
        symbolic_math = Num(convert(Num, extensive_math,
            convert_time = (x::SBML.MathTime) -> Catalyst.DEFAULT_IV))

        rstoich = reaction.reactants
        pstoich = reaction.products
        if reaction.reversible
            symbolic_math = getunidirectionalcomponents(symbolic_math)
            kl_fw, kl_rv = [substitute(x, subsdict) for x in symbolic_math]
            reagents = getreagents(rstoich, pstoich, model)
            isnothing(reagents[1]) && isnothing(reagents[2]) && continue
            kl_fw, our = use_rate(kl_fw, reagents[1], reagents[3])
            kl_rv = from_noncombinatoric(kl_rv, reagents[3], our)
            push!(rxs, Catalyst.Reaction(kl_fw, reagents...; only_use_rate = our))

            reagents = getreagents(rstoich, pstoich, model; rev = true)
            kl_rv, our = use_rate(kl_rv, reagents[1], reagents[3])
            kl_rv = from_noncombinatoric(kl_rv, reagents[3], our)
            push!(rxs, Catalyst.Reaction(kl_rv, reagents...; only_use_rate = our))
        else
            kl = substitute(symbolic_math, subsdict)
            reagents = getreagents(rstoich, pstoich, model)
            isnothing(reagents[1]) && isnothing(reagents[2]) && continue
            kl, our = use_rate(kl, reagents[1], reagents[3])
            kl = from_noncombinatoric(kl, reagents[3], our)
            push!(rxs, Catalyst.Reaction(kl, reagents...; only_use_rate = our))
        end
    end
    rxs
end

function from_noncombinatoric(rl::Num, stoich::Union{Vector{<:Real},Nothing}, only_use_rate::Bool)
    if !isnothing(stoich) && !only_use_rate
        coef = 1
        for s in stoich
            isone(s) && continue
            coef *= factorial(s)
        end
        !isone(coef) && (rl *= coef)
    end
    rl
end

""" Get kineticLaw for use in MTK.Reaction """
function use_rate(kl::Num, react::Union{Vector{Num},Nothing}, stoich::Union{Vector{<:Real},Nothing})
    rate_const = getmassaction(kl, react, stoich)
    if !isnan(rate_const)
        kl = rate_const
        our = false
    else
        our = true
    end
    return (kl, our)
end

""" Get reagents """
function getreagents(rstoichdict::Dict{String,<:Real}, pstoichdict::Dict{String,<:Real}, model::SBML.Model; rev = false)
    reactants = Num[]
    products = Num[]
    rstoich = Float64[]
    pstoich = Float64[]

    if rev
        tmp = rstoichdict
        rstoichdict = pstoichdict
        pstoichdict = tmp
    end

    for (k, v) in rstoichdict
        iszero(v) && @error("Stoichiometry of $k must be non-zero")
        push!(reactants, create_var(k, Catalyst.DEFAULT_IV))
        push!(rstoich, v)
        if model.species[k].boundary_condition == true
            push!(products, create_var(k, Catalyst.DEFAULT_IV))
            push!(pstoich, v)
        end
    end

    for (k, v) in pstoichdict
        iszero(v) && @error("Stoichiometry of $k must be non-zero")
        if model.species[k].boundary_condition != true
            push!(products, create_var(k, Catalyst.DEFAULT_IV))
            push!(pstoich, v)
        end
    end

    if (length(reactants) == 0)
        reactants = nothing
        rstoich = nothing
    end
    if (length(products) == 0)
        products = nothing
        pstoich = nothing
    end
    (reactants, products, rstoich, pstoich)
end

""" Infer forward and reverse components of bidirectional kineticLaw """
function getunidirectionalcomponents(bidirectional_math)
    err = "Cannot separate bidirectional kineticLaw `$bidirectional_math` to forward and reverse part. Please make reaction irreversible or rearrange kineticLaw to the form `term1 - term2`."
    bidirectional_math = Symbolics.tosymbol(bidirectional_math)
    bidirectional_math = simplify(bidirectional_math; expand = true)
    if (bidirectional_math isa Union{Real,Symbol}) || (SymbolicUtils.operation(bidirectional_math) != +)
        throw(ErrorException(err))
    end
    terms = SymbolicUtils.arguments(bidirectional_math)
    fw_terms = []
    rv_terms = []
    for term in terms
        if (term isa SymbolicUtils.Mul) && (term.coeff < 0)
            push!(rv_terms, Num(-term))  # PL: @Anand: Perhaps we should to create_var(term) or so?
        else
            push!(fw_terms, Num(term))  # PL: @Anand: Perhaps we should to create_var(term) or so?
        end
    end
    if (length(fw_terms) != 1) || (length(rv_terms) != 1)
        throw(ErrorException(err))
    end
    return (fw_terms[1], rv_terms[1])
end

""" Extract paramap from Model """
function get_paramap(model)
    paramap = Pair{Num,Float64}[]
    for (k, v) in model.parameters
        push!(paramap, Pair(create_param(k), v[1])) # [1] index drops unit
    end
    for (k, v) in model.compartments
        if !isnothing(v.size)
            push!(paramap, Pair(create_param(k), v.size))
        end
    end
    paramap
end

get_u0map(model) = [create_var(k, Catalyst.DEFAULT_IV) => v for (k, v) in SBML.initial_amounts(model, convert_concentrations = true)]
ModelingToolkit.defaults(model::SBML.Model) = ModelingToolkit._merge(get_u0map(model), get_paramap(model))

""" Get rate constant of mass action kineticLaws """
function getmassaction(kl::Num, reactants::Union{Vector{Num},Nothing}, stoich::Union{Vector{<:Real},Nothing})
    function check_args(x::SymbolicUtils.Symbolic{Real})
        for arg in SymbolicUtils.arguments(x)
            if isnan(check_args(arg)) || isequal(arg, Catalyst.DEFAULT_IV)
                return NaN
            end
        end
        return 0
    end
    check_args(x::Term{Real,Nothing}) = NaN  # Variable leaf node
    check_args(x::Sym{Real,Base.ImmutableDict{DataType,Any}}) = 0  # Parameter leaf node
    check_args(x::Real) = 0  # Real leaf node
    check_args(x) = throw(ErrorException("Cannot handle $(typeof(x)) types."))  # Unknow leaf node
    if isnothing(reactants) && isnothing(stoich)
        rate_const = kl
    elseif isnothing(reactants) | isnothing(stoich)
        throw(ErrorException("`reactants` and `stoich` are incosistent: `reactants` are $(reactants) and `stoich` is $(stoich)."))
    else
        rate_const = SymbolicUtils.simplify_fractions(kl / *((.^(reactants, stoich))...))
    end
    isnan(check_args(rate_const.val)) ? NaN : rate_const
end

function create_var(x)
    sym = Symbol(x)
    Symbolics.unwrap(first(@variables $sym))
end
function create_var(x, iv)
    sym = Symbol(x)
    Symbolics.unwrap(first(@variables $sym(iv)))
end
function create_param(x)
    sym = Symbol(x)
    Symbolics.unwrap(first(@parameters $sym))
end

function get_rules(model)
    subsdict = _get_substitutions(model)
    # these three go into `constraints` field of ReactionSystem
    obseqs = Equation[]
    algeqs = Equation[]
    raterules = Equation[]

    rules = model.rules
    for r in rules
        if r isa SBML.AlgebraicRule
            push!(algeqs, 0 ~ convert(Num, r.math))
        elseif r isa SBML.AssignmentRule
            push!(obseqs, assignmentrule_to_obseq(model, r))
        elseif r isa SBML.RateRule
            push!(raterules, raterule_to_diffeq(model, r))
        else
            error()
        end
    end
    algeqs, obseqs, raterules = map(x -> substitute(x, subsdict), (algeqs, obseqs, raterules))
    algeqs, obseqs, raterules
end

function assignmentrule_to_obseq(model, rule)
    if haskey(model.species, rule.id)
        sym = Symbol(rule.id)
        var = Symbolics.unwrap(first(@variables $sym(Catalyst.DEFAULT_IV)))
        assignment = Num(convert(Num, rule.math, convert_time = (x::SBML.MathTime) -> Catalyst.DEFAULT_IV))
        return var ~ assignment
    elseif haskey(model.compartments, rule.id)
        error("not handling setting compartment volumes rn")
    elseif haskey(model.parameters, rule.id)
        error("not handling setting non-constant parameters rn")
    else
        error()
    end
end

function raterule_to_diffeq(model, rule)
    # the rule.id can be a species, speciesRef, compartment, or param. currently not doingthis
    # for now im just doing species
    D = Differential(Catalyst.DEFAULT_IV)
    if haskey(model.species, rule.id)
        sym = Symbol(rule.id)
        var = Symbolics.unwrap(first(@variables $sym(Catalyst.DEFAULT_IV)))
        assignment = Num(convert(Num, rule.math, convert_time = (x::SBML.MathTime) -> Catalyst.DEFAULT_IV))
        return D(var) ~ assignment
    elseif haskey(model.comnpartments, rule.id)
        error("not handling setting compartment volumes rn")
    elseif haskey(model.parameters, rule.id)
        error("not handling setting non-constant parameters rn")
    else
        error()
    end
end


text:
module Util

using Nullables
using Genie, Genie.Loggers

export expand_nullable, _!!, _!_, get_nested_field, get_deepest_module, DynamicField, psst, time_to_unixtimestamp, reload

mutable struct DynamicField{T}
  field::T
end

"""
    add_quotes(str::String) :: String

Adds quotes around `str` and escapes any previously existing quotes.
"""
function add_quotes(str::String) :: String
  if ! startswith(str, "\"")
    str = "\"$str"
  end
  if ! endswith(str, "\"")
    str = "$str\""
  end

  str
end

"""
    strip_quotes(str::String) :: String

Unquotes `str`.
"""
function strip_quotes(str::String) :: String
  if is_quoted(str)
    str[2:end-1]
  else
    str
  end
end

"""
    is_quoted(str::String) :: Bool

Checks weather or not `str` is quoted.
"""
function is_quoted(str::String) :: Bool
  startswith(str, "\"") && endswith(str, "\"")
end

"""
    expand_nullable{T}(value::Nullable{T}, default::T) :: T

Returns `value` if it is not `null` - otherwise `default`.
"""
function expand_nullable(value::T)::T where T
  value
end
function expand_nullable(value::Nullable{T}, default::T)::T where T
  if isnull(value)
    default
  else
    Base.get(value)
  end
end

"""
    _!!{T}(value::Nullable{T}) :: T

Shortcut for `Base.get(value)`.
"""
function _!!(value::Nullable{T})::T where T
  Base.get(value)
end

"""
    _!_{T}(value::Nullable{T}, default::T) :: T

Shortcut for `expand_nullable(value, default)`.
"""
function _!_(value::Nullable{T}, default::T)::T where T
  expand_nullable(value, default)
end

"""
    file_name_without_extension(file_name, extension = ".jl") :: String

Removes the file extension `extension` from `file_name`.
"""
function file_name_without_extension(file_name, extension = ".jl") :: String
  file_name[1:end-length(extension)]
end

"""
    function walk_dir(dir, paths = String[]; only_extensions = ["jl"], only_files = true, only_dirs = false) :: Vector{String}

Recursively walks dir and `produce`s non directories. If `only_files`, directories will be skipped. If `only_dirs`, files will be skipped.
"""
function walk_dir(dir, paths = String[]; only_extensions = ["jl"], only_files = true, only_dirs = false) :: Vector{String}
  f = readdir(abspath(dir))
  for i in f
    full_path = joinpath(dir, i)
    if isdir(full_path)
      ! only_files || only_dirs && push!(paths, full_path)
      walk_dir(full_path, paths)
    else
      only_dirs && continue

      (last(split(i, ['.'])) in only_extensions) || isempty(only_extensions) && push!(paths, full_path)
    end
  end

  paths
end

"""
    get_nested_field(path::String, depth::Int = 1, parent::Module) :: DynamicField

Returns the deepest nested field in a `path` of form `Module.Module.Module.function`, wrapped into a `DynamicField`.
"""
function get_nested_field(path::String, depth::Int, parent::Module) :: DynamicField
  parts = split(path, ".")

  new_parent = getfield(parent, Symbol(parts[depth]))

  if length(parts)-1 > depth
    get_nested_field(path, depth+1, new_parent)
  else
    return getfield(new_parent, Symbol(parts[depth+1])) |> DynamicField
  end
end

"""
    get_deepest_module(path::String, depth::Int = 1, parent::Module) :: Module

Returns the deepest nested `module` in a `path` of form `Module.Module.Module`.
"""
function get_deepest_module(path::String, depth::Int, parent::Module) :: Module
  parts = split(path, ".")

  new_parent = getfield(parent, Symbol(parts[depth]))

  if length(parts)-2 > depth
    get_deepest_module(path, depth+1, new_parent)
  else
    return new_parent
  end
end

"""
    psst(f::Function)

Invokes `f` while supressing all debugging output for the duration of the invocation.
"""
function psst(f::Function)
  Genie.config.suppress_output = true
  result = f()
  Genie.config.suppress_output = false

  result
end

"""
    kill(t::Task) :: Nothing

Kills `Task` `t` by forcing it to throw an `InterruptException`.
"""
function kill(t::Task) :: Nothing
  Base.throwto(t, InterruptException())
end

"""
    reload_modules(dir::String, md::Module = current_module()) :: Bool

Reloads all the modules in the specified `dir` in the scope of `md`.
Returns `true` if any modules were reloaded, `false` otherwise.
"""
function reload_modules(dir::String, md::Module = current_module()) :: Bool
  status = false

  for file in readdir(dir)
    if isfile(joinpath(dir, file)) && endswith(file, ".jl") && isdefined(Symbol(file[1:end-3]))
      Core.eval(md, :(reload("$($file[1:end-3])")) )
      status = true
    end
  end

  status
end

"""
    reload_modules(dirs::Vector{String}, md::Module = current_module()) :: Bool

Reloads all the modules in all the directories specified in `dirs` in the scope of `md`.
Returns `true` if any modules were reloaded, `false` otherwise.
"""
function reload_modules(dirs::Vector{String}, md::Module = current_module()) :: Bool
  status = false

  for dir in dirs
    reload_modules(dir, md) && (status = true)
  end

  status
end

"""

"""
function time_to_unixtimestamp(t::Float64)
  floor(t) |> Int
end
function time_to_unixtimestamp()
  time_to_unixtimestamp(time())
end

"""
"""
function reload(m::Module)
  reload(string(m))
end

end


text:
inputs = [1.0, 2.0, 3.0, 2.5]

weights1 = [0.2, 0.8, -0.5, 1.0]
weights2 = [0.5, -0.91, 0.26, -0.5]
weights3 = [-0.26, -0.27, 0.17, 0.87]

bias1 = 2.0
bias2 = 3.0
bias3 = 0.5

output = [inputs[1]*weights1[1] + inputs[2]*weights1[2] + inputs[3]*weights1[3] + inputs[4]*weights1[4] + bias1,
          inputs[1]*weights2[1] + inputs[2]*weights2[2] + inputs[3]*weights2[3] + inputs[4]*weights2[4] + bias2,
          inputs[1]*weights3[1] + inputs[2]*weights3[2] + inputs[3]*weights3[3] + inputs[4]*weights3[4] + bias3]
println(output)


text:
using MPSLabXD
using Documenter

makedocs(;
    modules=[MPSLabXD],
    authors="PhysicsCodesLab",
    sitename="MPSLabXD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        mathengine = MathJax()
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => ["manual/MPS.md","manual/fDMRG.md","manual/MPO.md"],
        "Index" => ["index/index.md"]
    ],
)

deploydocs(;
    repo="github.com/PhysicsCodesLab/MPSLabXD.jl.git",
    devbranch="master",
)


text:
using Test

@testset "JuliaExpr.emit" begin
    include("julia/node.jl")
    include("julia/arg.jl")
    include("julia/dash.jl")
    include("julia/optional.jl")
    include("julia/options.jl")
    include("julia/vararg.jl")
    include("julia/exception.jl")
    include("julia/plugin.jl")
end

@testset "ZSHCompletions.emit" begin
    include("zsh/zsh.jl")
end


text:
"""
    ParseNexus(filename::String)

This function parses a NEXUS file which stores the input for the MCMC compuation.
The file should follow the conventions used for MrBayes.

Returns ntax, nchar, gap, and missing_representation values; returns Dataframe storing language names and data.

* `filename` : NEXUS file to be parsed.
"""
function ParseNexus(filename::String)
    open(filename, "r") do file
        global content = readlines(file)
    end # do

    if lowercase(content[1]) != "#nexus"
        throw("$filename is not a Nexus file!")
    end # if

    while true
        line = popfirst!(content)
        if lowercase(line) == "begin data;"
            # here comes the important information, so break the loop
            break
        end #if
    end # while
    # get meta info
    ntax, nchar, gap, missing_representation, symbols = extract_meta_info(content)
    df = create_nexusdf(content)

    out_symbols = symbols == "NOSYMBOLS" ? get_alphabet(df, gap, missing_representation) : [string(s) for s in symbols]

    return ntax, nchar, gap, missing_representation, out_symbols, df
end # function ParseNexus

function get_alphabet(df::DataFrame, gap::String, missing_representation::String)
    d = df[:, :Data]
    alphabet = Set{String}()
    for row in d
        for entry in row
            entry = string(entry)
            if !(entry == gap) && !(entry == missing_representation)
                push!(alphabet, entry)
            end
        end
    end
    sort([i for i in alphabet])
end

"""
    extract_meta_info(content::Array{String})

This function extracts some meta information from the content of the nexus file.
It "eats-up" the stack.

Returns values derived from metadata (ntax, nchar, gap, missing_representation). Used in ParseNexus().

* `content` : Array of Strings; Strings are read from NEXUS file in ParseNexus().
"""
function extract_meta_info(content::Array{String})
    lct::Int64 = 0
    ntax::Int64 = 0
    nchar::Int64 = 0
    gap::String = "-"
    symbols::String= "NOSYMBOLS"
    missing_representation::String = "?"
    while true
        line = popfirst!(content)
        if lowercase(line)=="matrix"
            # meta info is done, data begins now
            break
        else
            splitted_line = split(line)
            for entry in splitted_line
                info = split(entry, "=")
                info = lowercase.(info)
                if length(info) != 1
                    choped = info[2]
                    if endswith(choped, ";")
                        choped = chop(choped)
                    end #if
                    k_word = info[1]
                    if k_word == "ntax"
                        ntax = parse(Int64, choped)
                    elseif k_word == "nchar"
                        nchar = parse(Int64, choped)
                    elseif k_word == "gap"
                        gap = choped
                    elseif k_word == "missing"
                        missing_representation = choped
                    elseif k_word == "symbols"
                        symbols = strip(String(choped), '\"')
                    else
                        @warn "Keyword $k_word not understood, will be ignored"
                    end # if
                end # if
            end # for
        end # if
    end # while
    return ntax, nchar, gap, missing_representation, symbols
end # function extract_meta_info

"""
    create_nexusdf(filecontent::Array{String})::DataFrame

This function creates a DataFrame of the actual data. Used in ParseNexus().

Returns DataFrame of language names and data derived from NEXUS file.

* `filecontent` : Array of Strings; Strings are read from NEXUS file in ParseNexus().
"""
function create_nexusdf(filecontent::Array{String})::DataFrame
    df = DataFrame(Language=String[], Data=String[])
    while true
        line = popfirst!(filecontent)
        if line == ""
            continue
        elseif line == ";"
            break
        else
            push!(df, split(line))
        end # if
    end # while
    return df
end # function create_nexusdf


text:
module ARMA

using Polynomials, NLsolve, HDF5
using LinearAlgebra
using Statistics
using ToeplitzMatrices

export
    estimate_covariance,
    fitARMA,
    fit_exponentials,
    ARMAModel,
    generate_noise,
    model_covariance,
    model_psd,
    ARMASolver,
    whiten,
    unwhiten,
    mult_covariance,
    solve_covariance,
    inverse_covariance,
    toeplitz_whiten

# ARMA.jl includes the basic ARMA models and their use
# model_selection.jl includes tools for choosing a model (p,q) order
#   and for fitting the best model of that order.

include("model_selection.jl")

"""
    ARMAModel(...)

An Autoregressive Moving-Average model of order (p,q).

There are three alternative ways to construct the model:

1. `ARMAModel(thetacoef::Vector, phicoef::Vector)`
2. `ARMAModel(roots_::Vector, poles::Vector, variance::Float64)`
3. `ARMAModel(bases::Vector, amplitudes::Vector, covarianceIV::Vector)`

The first method specifies the model by the coefficients of its rational
function representation. The ordering of the coefficients is constant first,
high-order last: `c[1]+c[2]*z+c[3]*z^2...` Here, `thetacoef` represents the
numerator of the forward transform and `phicoef` the denominator. By convention,
the value `phicoef[1]` is 1.0, but both polynomials will be rescaled to achieve
this before calling the inner constructor. To ensure a stable, invertible ARMA
model, the roots of both polynomials should have absolute values greater than 1.

The second method specifies the model by the roots and poles of its rational
function representation, along with a variance (that is, covariance at lag=0).
The latter is necessary to specify the overall scale of the model, because roots
and poles are insensitive to scale. All roots and poles must have absolute value
greater than 1, to ensure a stable, invertible ARMA model.

The third specifies the model by its covariance function, represented as a sum
of (potentially complex) exponentials, with 0 or more initial values of the
covariance (`covarianceIV`) that violate the general pattern.  For `ARMAModel(b,
a, cIV)`, the covariance at lag `j-1` is either `cIV[j]` when `j<=length(cIV)`,
or else `sum(b .* (a.^t))`. For any base and amplitude pair (bi,ai) where either
has a non-zero imaginary part, it is required that (conj(bi),conj(ai)) also be
among the base and amplitude pairs. This requirement ensures that the covariance
is everywhere real. Also, we require `|b|<1` for all bases to ensure a stable
model.

Methods 1 and 2 can create ARMA(p,q) models of arbitrary order. A pure AR(p)
model with q=0 can be generated by passing `thetacoef` of length 1 or
`roots_=[]` to the methods 1 and 2 constructors, respectively. Similarly, a pure
MA(q) model with p=0 can be generated using `phicoef` of length 1 or `poles=[]`.
Method 3 generates an ARMA(p,q) model where p is the length of both the `bases`
and `amplitudes` vectors, and where q = p-1+length(covarianceIV).  While models
of MA order `q<p-1` are possible in principle, they are not well-matched to the
sum-of-exponentials representation, because setting `q<p-1` imposes complicated
restrictions on the relative amplitudes of the exponentials.

When an `ARMAModel` is constructed by any of these methods, all three
representations are computed and stored.

See also: [`ARMASolver`](@ref) for fast mutliplies and solves of either the
noise model covariance matrix or its Cholesky factor.
"""
mutable struct ARMAModel
    p         ::Int
    q         ::Int
    roots_    ::Vector{ComplexF64}
    poles     ::Vector{ComplexF64}
    thetacoef ::Vector{Float64}
    phicoef   ::Vector{Float64}
    covarIV   ::Vector{Float64}
    expbases  ::Vector{ComplexF64}
    expampls  ::Vector{ComplexF64}

    function ARMAModel(p,q,roots_,poles,thetacoef,phicoef,covarIV,expbases,expampls)
        @assert p == length(poles)
        @assert q == length(roots_)
        @assert all(abs2.(poles) .> 1)
        @assert all(abs2.(roots_) .> 1)
        @assert p+1 == length(phicoef)
        @assert q+1 == length(thetacoef)
        @assert p == length(expbases)
        @assert p == length(expampls)
        @assert length(covarIV) >= 1+q
        @assert length(covarIV) >= p
        @assert thetacoef[1] > 0
        @assert phicoef[1] == 1.0
        # Note that these consistency checks don't cover everything. Specifically, we
        # do not test the consistency of the 3 representations with each other. That's
        # done only in the 3 outer constructors.
        new(p,q,roots_,poles,thetacoef,phicoef,covarIV,expbases,expampls)
    end
end

function Base.show(io::IO, m::ARMAModel)
    print(io, "ARMAModel(p=$(m.p), q=$(m.q))\n")
    print(io, "[Amplitude, Exponential base]\n")
    for i=1:m.p
        print(io, "(A,B)[$i]: $(m.expampls[i]),  $(m.expbases[i])\n")
    end
    print(io, "Initial covar: $(m.covarIV)\n")
end

"Go from theta,phi polynomials to the sum-of-exponentials representation.
Returns (covar_initial_values, exponential_bases, exponential_amplitudes)."
function _covar_repr(thetacoef::Vector, phicoef::Vector)
    roots_ = roots(Polynomial(thetacoef))
    poles = roots(Polynomial(phicoef))
    expbases = 1.0 ./ poles
    q = length(roots_)
    p = length(poles)
    n = max(p,q)

    # Find the initial, exceptional values.

    # psi comes from Brockwell 3.3.3 or delta in Pollock 17.87.
    phi = fill(0.0, n+1)
    phi[1:p+1] = phicoef
    P = TriangularToeplitz(phi, :L) # lower-triangular toeplitz
    theta = fill(0.0, n+1)
    theta[1:q+1] = thetacoef
    psi = P \ theta

    T = fill(0.0, n+1, n+1) # A Hankel matrix containing theta in first col.
    for r=1:n+1
        for c=1:n+2-r
            T[r,c] = theta[c+r-1]
        end
    end
    P2 = fill(0.0, n+1, n+1)
    for r=1:n+1
        for i=1:p+1
            c = r-i+1
            if c<1; c = 2-c; end
            P2[r,c] += phi[i]
        end
    end
    gamma = P2 \ (T*psi)

    # Work in the range where exceptional values no longer hold.
    # XI[r,c] = expbases[c]^(r-1-lowestpower)
    # lowestpower is chosen to make the top row be beyond the 1+q-p expectional values.
    XI = Array{ComplexF64}(undef, p, p)
    lowestpower = p >= q ? 1 : 1+q-p
    for c=1:p
        XI[:,c] = expbases[c] .^ collect(lowestpower:p+lowestpower-1)
    end
    expampls = XI \ model_covariance(gamma, phicoef, p+lowestpower)[1+lowestpower:p+lowestpower]
    gamma, expbases, expampls
end

# Construct from theta and phi polynomial representation
function ARMAModel(thetacoef::Vector, phicoef::Vector)
    # Normalize the input coefficient vectors
    theta = thetacoef / phicoef[1]
    if theta[1] < 0
        theta = -theta
    end
    phi = phicoef / phicoef[1]
    roots_ = roots(Polynomial(theta))
    poles = roots(Polynomial(phi))
    @assert all(abs2.(roots_) .> 1)
    @assert all(abs2.(poles) .> 1)
    q = length(roots_)
    p = length(poles)

    covarIV, expbases, expampls = _covar_repr(theta, phi)
    ARMAModel(p,q,roots_,poles,theta,phi,covarIV,expbases,expampls)
end

WhiteModel() = ARMAModel([1.0], [1.0])

"Form the coefficients of a polynomial from the given roots `r`.
It is assumed that the coefficients are real, so only the real part is kept.
The zero-order term (the first) has coefficient +1. Fails if 0 is a root."
function polynomial_from_roots(r::AbstractVector)
    pr = prod(r)
    @assert abs(imag(pr)/real(pr)) < 1e-10
    coef = real(fromroots(r).coeffs)
    coef / coef[1]
end

# Construct from roots-and-poles representation. We also need the gamma_0 value
# (the process variance) to set the scale of the model, as roots-and-poles omits this.
function ARMAModel(roots_::AbstractVector, poles::AbstractVector, variance::Real)
    @assert all(abs2.(roots_) .> 1)
    @assert all(abs2.(poles) .> 1)
    # The product of the roots and the product of the poles needs to be real.
    pr = prod(roots_)
    pp = prod(poles)
    @assert abs(imag(pr)/real(pr)) < 1e-10
    @assert abs(imag(pp)/real(pp)) < 1e-10
    q = length(roots_)
    p = length(poles)

    # Construct normalized MA and AR polynomials.
    # Note that we will NOT have the proper scale at first.
    thetacoef = polynomial_from_roots(roots_)
    phicoef = polynomial_from_roots(poles)

    covarIV, expbases, expampls = _covar_repr(thetacoef, phicoef)

    # Now that we know the normalized model's covariance, rescale
    # theta, covariance, and the exponential amplitudes to correct it.
    scale = sqrt(variance / covarIV[1])
    covarIV .*= variance / covarIV[1]
    ARMAModel(p,q,roots_,poles,thetacoef*scale,phicoef,covarIV,expbases,expampls*scale)
end

function solveGammaMA(gammaMA::AbstractVector; maxiter=1000)
    # Solve the system gammaMA[i] = sum mu[j]*mu[j+i] for mu.
    # Goal: find a root ; if we cannot, at least find the lowest-cost
    # value of x (here, meaning the lowest mean absolute residual).

    # Strategy is step 1: iterative estimation of mu.
    # Step 2: full nonlinear solve for mu, if successful.

    function f!(residual, x)
        q = length(residual)-1
        for t = 0:q
            residual[t+1] = gammaMA[t+1]
            for j = 0:q-t
                residual[t+1] -= x[j+1]*x[j+t+1]
            end
        end
        residual
    end

    function cost(x)
        residual = fill(zero(eltype(x)), size(x))
        f!(residual, x)
        mean(abs.(residual))
    end

    n = length(gammaMA)
    theta = fill(0.0, n)
    theta[1] = 1
    var = gammaMA[1]
    besttheta = copy(theta)
    bestcost = Inf
    consecincreases = 0

    # Iterate until result gets worse 10 times in a row, the mean absolute deviation
    # from gamma is less than gammaMA[1]/10^6, or we reach the iteration limit.
    for iter = 1:maxiter
        theta[n] = gammaMA[n]/var
        for i=n-1:-1:2
            theta[i] = gammaMA[i]/var - dot(theta[2:n+1-i], theta[1+i:n])
        end
        var = gammaMA[1]/sum(theta.^2)
        c = cost(theta * sqrt(var))
        if c < 1e-6*gammaMA[1]
            besttheta[:] = theta[:]*sqrt(var)
            # println("Basic iteration succeeds")
            break
        end
        if c <= bestcost
            bestcost = c
            besttheta[:] = theta[:]*sqrt(var)
        else
            consecincreases += 1
            if consecincreases >= 10
                # println("10 consecutive increases as of iter = $(iter)")
                break
            end
        end
    end

    # Step 2: try a full nonlinear solver. Ignore result if it's worse or errors.
    try
        resid = 0*theta
        f!(resid, besttheta)
        results = nlsolve(f!, besttheta, iterations=maxiter)
        costnl = cost(results.zero)
        if costnl < bestcost
            # println("Solver succeeded: costnl=$(costnl) vs $(bestcost) ")
            return results.zero
        else
            # println("Solver costnl=$(costnl) vs $(bestcost)")
        end
    catch e
    end
    return besttheta
end

# Construct ARMAModel from a sum-of-exponentials representation, along with
# zero or more exceptional values of the initial covariance, covarIV.
# It is allowed for covarIV==[], but it can't be optional,
# else this constructor will get mixed up with the thetacoef,phicoef one.
# The model will have order ARMA(p,q) where p=length(bases)=length(amplitudes),
# and q=p-1+length(covarIV).

function ARMAModel(bases::AbstractVector, amplitudes::AbstractVector, covarIV::AbstractVector)
    p = length(bases)
    @assert p == length(amplitudes)
    q = p-1+length(covarIV)

    # Find the covariance from lags 0 to p+q. Call it gamma
    gamma = fill(0.0, 1+p+q+50)
    t = 0:length(gamma)-1
    for i=1:p
        gamma += real(amplitudes[i] * (bases[i] .^ t))
    end
    gamma[1:length(covarIV)] = covarIV

    # Find the phi polynomial
    poles = 1.0 ./ bases
    phicoef = polynomial_from_roots(poles)

    # Find the nonlinear system of equations for the theta coefficients
    gammaMA = fill(0.0, 1+q)
    for t=0:q
        for i=1:p+1
            for j=1:p+1
                gammaMA[t+1] += phicoef[i]*phicoef[j]*gamma[1+abs(t+i-j)]
            end
        end
    end

    besttheta = solveGammaMA(gammaMA)
    roots_ = roots(Polynomial(besttheta))

    # Now a clever trick: any roots r of the MA polynomial that are INSIDE
    # the unit circle can be replaced by 1/r and yield the same covariance.
    @assert q == length(roots_)
    for i=1:q
        if abs2(roots_[i]) < 1
            roots_[i] = 1.0/roots_[i]
        end
    end

    # Replace thetacoef, in case any roots moved
    thetacoef = polynomial_from_roots(roots_)

    # One last problem is that theta polynomial is normalized, with theta[end]=1,
    # which we don't want. Our approach is just to see what the normalized theta
    # would give for gamma[1] and rescale.
    gammanorm,_,_ = _covar_repr(thetacoef,phicoef)
    thetacoef .*= sqrt(gamma[1]/gammanorm[1])

    ARMAModel(p,q,roots_,poles,thetacoef,phicoef,gamma[1:max(p,q+1)],bases,amplitudes)
end

"""
    hdf5save(output, model::ARMAModel)

Store the `model` to an HDF5 file.

`output` may be an HDF5 file open for writing, an HDF5 group within such a file,
or the file name of an HDF5 file which will be created.
"""
function hdf5save(output::HDF5.H5DataStore, model::ARMAModel)
    GROUPNAME = "ARMAModel"
    if haskey(output, GROUPNAME)
        o_delete(output, GROUPNAME)
    end
    grp = create_group(output, GROUPNAME)
    grp["basesR"] = real(model.expbases)
    grp["basesI"] = imag(model.expbases)
    grp["amplitudesR"] = real(model.expampls)
    grp["amplitudesI"] = imag(model.expampls)
    grp["covarIV"] = model.covarIV[1:max(model.p, model.q+1)-model.p]
    model
end

function hdf5save(filename::AbstractString, model::ARMAModel)
    h5open(filename, "w") do f
        hdf5save(f, model)
    end
end

"""
    model = hdf5load(input)

Load and return an `ARMAModel` object from `input`. The argument `input` can be
the "ARMAModel" group, the group that contains it, or the name of an HDF5 file
(group "ARMAModel" would need to be at the root of the file).
"""
function hdf5load(input::HDF5.H5DataStore)
    if haskey(input, "ARMAModel")
        return hdf5load(input["ARMAModel"])
    end
    bases = input["basesR"][:] + 1im*input["basesI"][:]
    ampls = input["amplitudesR"][:] + 1im*input["amplitudesI"][:]
    covarIV = input["covarIV"][:]

    ARMAModel(bases, ampls, covarIV)
end

function hdf5load(filename::AbstractString)
    h5open(filename, "r") do f
        return hdf5load(f)
    end
end

"""
    generate_noise(m, N)

Generate a simulated noise timeseries of length `N` from an ARMAModel `m`.
"""
function generate_noise(m::ARMAModel, N::Int)
    # eps = white N(0,1) noise; x = after MA process; z = after inverting AR
    eps = randn(N+m.q)
    eps[1:m.p] = 0
    x = fill(0.0, N)
    z = fill(0.0, N)
    for i=1:m.q+1
        x += eps[i:end+i-m.p-1] * m.thetacoef[i]
    end
    for j=1:m.p
        z[j] = x[j]
        for i = 2:j
            z[j] -= m.phicoef[i] * z[j-i+1]
        end
    end
    for j=1+m.p:N
        z[j] = x[j]
        for i = 2:m.p+1
            z[j] -= m.phicoef[i] * z[j-i+1]
        end
    end
    z
end

"""
    model_covariance(...)

Compute the ARMA model's model covariance function. Methods:

1. `model_covariance(m::ARMAModel, N::Int)`
2. `model_covariance(covarIV::Vector, phicoef::Vector, N::Int)`

The first returns covariance of model `m` from lags 0 to `N-1`. The second
returns the covariance given the initial exceptional values of covariance
`covarIV` and the coefficients `phicoef` of the Phi polynomial. (From these
coefficients, a recusion allows computation of covariance beyond the initial
values.)
"""
function model_covariance(covarIV::AbstractVector, phicoef::AbstractVector, N::Int)
    if N < length(covarIV)
        return covarIV[1:N]
    end
    @assert(length(covarIV) >= length(phicoef)-1)
    if phicoef[1] != 1.0
        phicoef = phicoef/phicoef[1]
    end
    covar = fill(0.0, N)
    covar[1:length(covarIV)] = covarIV[1:end]

    # Use Pollock eq 17.86
    for i = length(covarIV)+1:N
        for j = 1:length(phicoef)-1
            covar[i] -= phicoef[j+1] * covar[i-j]
        end
    end
    covar
end
model_covariance(m::ARMAModel, N::Int) = model_covariance(m.covarIV, m.phicoef, N)

"""
    model_psd(m,...)

The ARMA model's power spectral density function. Methods:

1. `model_psd(m::ARMAModel, freq::Vector)`
1. `model_psd(m::ARMAModel, N::Int)`

The first returns PSD at the given frequencies `freq` (frequency 0.5 is the
critical, or Nyquist, frequency). The second returns PSD at `N` equally-spaced
frequencies from 0 to 0.5.
"""
function model_psd(m::ARMAModel, freq::AbstractVector)
    z = exp.(-2im*pi *freq)
    N = length(freq)
    numer = fill(complex(m.thetacoef[1]), N)
    for i=1:m.q
        numer .+= m.thetacoef[i+1] * (z.^i)
    end
    denom = fill(complex(m.phicoef[1]), N)
    for i=1:m.p
        denom .+= m.phicoef[i+1] * (z.^i)
    end
    abs2.(numer ./ denom)
end

model_psd(m::ARMAModel, N::Int) = model_psd(m, collect(range(0, stop=0.5, length=N)))

"""
    toeplitz_whiten(m, timestream)

Approximately whiten the timestream using a Toeplitz matrix. This has the
consequence that a zero-padded delay of the input timestream is equivalent to a
zero-padded delay of the output.

Returns the whitened timestream.

If `timestream` is a matrix, then each column is whitened, and the result is
returned.

No Toeplitz matrix has the ability to make the input exactly white,
but for many purposes, the time-shift property is more valuable than
that exact whitening.
"""
function toeplitz_whiten(m::ARMAModel, timestream::AbstractVector)
    N = length(timestream)
    white = fill(0.0, N)

    # First, multiply the input by the AR matrix (a banded Toeplitz
    # matrix with the phi coefficients on the diagonal and first p
    # subdiagonals).
    # The result is the MA matrix times the whitened data.
    MAonly = m.phicoef[1] * timestream
    for i=1:m.p
        MAonly[1+i:end] .+= m.phicoef[i+1] * timestream[1:end-i]
    end

    # Second, solve the MA matrix (also a banded Toeplitz matrix with
    # q non-zero subdiagonals.)
    white[1] = MAonly[1] / m.thetacoef[1]
    if N==1
        return white
    end
    for i = 2:min(m.q, N)
        white[i] = MAonly[i]
        for j = 1:i-1
            white[i] -= white[j]*m.thetacoef[1+i-j]
        end
        white[i] /= m.thetacoef[1]
    end
    for i = m.q+1:N
        white[i] = MAonly[i]
        for j = i-m.q:i-1
            white[i] -= white[j]*m.thetacoef[1+i-j]
        end
        white[i] /= m.thetacoef[1]
    end
    white
end

function toeplitz_whiten(m::ARMAModel, M::AbstractMatrix)
    tw(v::AbstractVector) = toeplitz_whiten(m, v)
    mapslices(tw, M, dims=1)
end

include("exact_operations.jl")

end # module


text:
# Linear Discriminant Analysis

mutable struct LinearDiscriminantModel{T} <: DiscriminantModel{T}
    "Standard discriminant model parameters"
    Θ::DiscriminantParameters{T}
    "Whitening transformation"
    W::AbstractMatrix{T}
    "Canonical coordinate projection matrix with whitening"
    C::Union{Nothing,AbstractMatrix{T}}
    "Discriminant function intercept"
    δ₀::T
    function LinearDiscriminantModel{T}() where T
        new{T}(DiscriminantParameters{T}(), Matrix{T}(undef,0,0), nothing, zero(T))
    end
end

function canonical_coordinates!(LDA::LinearDiscriminantModel{T}) where T
    Θ = LDA.Θ

    m, p = check_dims(Θ.M, dims=Θ.dims)

    if p ≤ m-1
        # no dimensionality reduction is possible
        LDA.C = Matrix{T}(I, p, p)
    elseif Θ.dims == 1
        # center M by overall centroid
        M = broadcast!(-, similar(Θ.M, T), Θ.M, transpose(Θ.μ))

        # Σ_between = Mᵀdiag(π)M so need to scale by sqrt π
        broadcast!((πₖ, Mₖⱼ) -> √(πₖ)Mₖⱼ, M, Θ.π, M)
        UDVᵀ = svd!(M*LDA.W, full=false)

        # Extract m-1 components
        LDA.C = LDA.W*transpose(view(UDVᵀ.Vt, 1:m-1, :))
    else
        M = broadcast!(-, similar(Θ.M, T), Θ.M, Θ.μ)

        broadcast!((πₖ, Mⱼₖ) -> √(πₖ)Mⱼₖ, M, transpose(Θ.π), M)
        UDVᵀ = svd!(LDA.W*M, full=false)

        LDA.C = transpose(view(UDVᵀ.U, :, 1:m-1))*LDA.W
    end

    return LDA
end

"""
    fit!(LDA::LinearDiscriminantModel, Y::Vector, X::Matrix; dims::Integer=1, kwargs...)

Fit a linear discriminant model based on data matrix `X` and class indicator matrix `Y`
along dimensions `dims` and overwrite `LDA`.

# Keyword arguments
- `canonical::Bool=false`: computes canonical coordinates and performs dimensionality
  reduction if `true`
- `compute_covariance`: if true, computes `Σ`
- `centroids::Matrix`: specifies the class centroids. If `dims` is `1`, then each row
   represents a class centroid, otherwise each column represents a class centroid. If not
   specified, the centroids are computed from the data
- `priors::Vector`: specifies the class membership prior probabilties. If not specified,
  the class probabilities a computed based on the class counts
- `gamma::Real=0`: regularization parameter γ ∈ [0,1] shrinks the within-class covariance
  matrix towards the identity matrix scaled by the average eigenvalue of the covariance
  matrix
"""
function _fit!(LDA::LinearDiscriminantModel{T},
               y::Vector{<:Integer},
               X::Matrix{T},
               dims::Integer=1,
               canonical::Bool=false,
               compute_covariance::Bool=false,
               centroids::Union{Nothing,AbstractMatrix}=nothing,
               priors::Union{Nothing,AbstractVector}=nothing,
               gamma::Union{Nothing,Real}=nothing) where T
    Θ = LDA.Θ

    set_dimensions!(Θ, y, X, dims)
    set_gamma!(Θ, gamma)
    set_statistics!(Θ, y, X, centroids)
    set_priors!(Θ, priors)

    # Compute degress of freedom
    df = size(X, dims) - Θ.m

    if dims == 1
        # Overall centroid is prior-weighted average of class centroids
        Θ.μ = vec(transpose(Θ.π)*Θ.M)

        # Center the data matrix with respect to classes to compute whitening matrix
        X .-= view(Θ.M, y, :)

        # Compute covariance matrix if requested
        if compute_covariance
            Θ.Σ = lmul!(one(T)/df, transpose(X)*X)
        end
    else
        Θ.μ = Θ.M*Θ.π

        X .-= view(Θ.M, :, y)

        if compute_covariance
            Θ.Σ = lmul!(one(T)/df, X*transpose(X))
        end
    end

    # Use cholesky whitening if gamma is not specifed, otherwise svd whitening
    if Θ.Σ !== nothing
        LDA.W, LDA.δ₀ = whiten_cov_chol!(copy(Θ.Σ), Θ.γ, dims=dims)
    elseif Θ.γ === nothing
        LDA.W, LDA.δ₀ = whiten_data_chol!(X, dims=dims, df=df)
    else
        LDA.W, LDA.δ₀ = whiten_data_svd!(X, Θ.γ, dims=dims, df=df)
    end

    # Perform canonical discriminant analysis if applicable
    if canonical
        canonical_coordinates!(LDA)
    else
        LDA.C = nothing
    end

    Θ.fit = true

    return LDA
end

function lda(X::Matrix{T},
             y::Vector{<:Integer};
             dims::Integer=1,
             canonical::Bool=false,
             compute_covariance::Bool=false,
             centroids::Union{Nothing,AbstractMatrix}=nothing,
             priors::Union{Nothing,AbstractVector}=nothing,
             gamma::Union{Nothing,Real}=nothing) where T
    U = T == Float32 || T == Float64 ? T : Float64
    return _fit!(
        LinearDiscriminantModel{U}(),
        y,
        U == T ? copy(X) : convert(Matrix{U}, X),
        dims,
        canonical,
        compute_covariance,
        centroids,
        priors,
        gamma
    )
end

function discriminants!(Δ::Matrix{T}, LDA::LinearDiscriminantModel{T}, X::Matrix{T}) where T
    dims = LDA.Θ.dims

    M = LDA.Θ.M

    n, p = check_dims(X, dims=dims)
    m, p₂ = check_dims(M, dims=dims)
    n₂, m₂ = check_dims(Δ, dims=dims)

    alt_dims = dims == 1 ? 2 : 1

    p == p₂ || throw(DimensionMismatch("predictor count along dimension $(alt_dims) of " *
                                       "X must match dimension $(dims) of M (got $(p) " *
                                       "and $(p₂))"))
    n == n₂ || throw(DimensionMismatch("observation count along dimension $(dims) of X " *
                                       "must match dimension $(dims) of Δ (got $(n) " *
                                       "and $(n₂))"))
    m == m₂ || throw(DimensionMismatch("class count along dimension $(alt_dims) of Δ " *
                                       "must match dimension $(dims) of M (got $(m) " *
                                       "and $(m₂))"))

    is_row = dims == 1

    π = LDA.Θ.π
    W = LDA.C === nothing ? LDA.W : LDA.C

    X̃, M̃ = is_row ? (X*W, M*W) : (W*X, W*M)

    Z = similar(X̃)
    for k = 1:m
        if is_row
            Δₖ = view(Δ, :, k:k)
            μₖ = view(M̃, k:k, :)
        else
            Δₖ = view(Δ, k:k, :)
            μₖ = view(M̃, :, k:k)
        end

        broadcast!(-, Z, X̃, μₖ)
        sum!(abs2, Δₖ, Z)
        broadcast!((d², logπ) -> logπ - d²/2, Δₖ, Δₖ, log(π[k]))
    end

    return Δ
end


text:
using ForwardDiff
using Mimi
export grad_component

# Optimize a univariate objective, by constructing a close linear approximation
function univariatematrix(func, lb, ub, x0, tolerance)
    # Collect all derivatives
    result = DerivativeResult(x0)
    ForwardDiff.derivative!(result, func, x0);

    pointslope = ForwardDiff.derivative(result)
    pointvalue = ForwardDiff.value(result)

    pointsecond = ForwardDiff.derivative(x -> ForwardDiff.derivative(func, x), x0)

    # Determine span for which linear is good approx
    pointsafespan = sqrt(2 * tolerance / abs(pointsecond))
    leftsafespan = min(pointsafespan, tolerance / abs(func(lb) - -pointslope * (x0 - lb)))
    rightsafespan = min(pointsafespan, tolerance / abs(func(ub) - pointslope * (ub - x0)))

    if lb < x0 - leftsafespan
        # Create average slope from here to edge
        leftslope = (func(x0) - func(x0 - leftsafespan)) / leftsafespan
        # Recurse for areas left of here
        leftfs, leftmins, leftmaxs = univariatematrix(func, lb, x0 - leftsafespan,
                                                      x0 - leftsafespan, tolerance)
    else
        leftslope = (func(x0) - func(lb)) / (x0 - lb)
        leftfs = leftmins = leftmaxs = []
    end

    if ub > x0 + rightsafespan
        # Create average slope from here to edge
        rightslope = (func(x0 + rightsafespan) - func(x0)) / rightsafespan
        # Recurse for areas right of here
        rightfs, rightmins, rightmaxs = univariatematrix(func, x0 + rightsafespan, ub,
                                                         x0 + rightsafespan, tolerance)
    else
        rightslope = (func(ub) - func(x0)) / (ub - x0)
        rightfs = rightmins = rightmaxs = []
    end

    if x0 == lb
        fs = [rightslope; leftfs; rightfs]
        mins = [0.; leftmins; rightmins]
        maxs = [min(rightsafespan, ub - x0); leftmaxs; rightmaxs]
    elseif x0 == ub
        fs = [-leftslope; leftfs; rightfs]
        mins = [0.; leftmins; rightmins]
        maxs = [min(leftsafespan, x0 - lb); leftmaxs; rightmaxs]
    else
        # leftslope is negative, because if upward sloping then this would take us down
        fs = [-leftslope; rightslope; leftfs; rightfs]
        mins = [0.; 0.; leftmins; rightmins]
        maxs = [min(leftsafespan, x0 - lb); min(rightsafespan, ub - x0); leftmaxs; rightmaxs]
    end

    fs, mins, maxs
end

function grad_component(m::Model, component::Symbol, variable::Symbol, parameter::Symbol)
    fdf = unaryobjective(m, [component], [parameter], m -> vec(m[component, variable]))
    A = ForwardDiff.jacobian(fdf, convert(Vector{Float64}, vec(m.md.external_params[parameter].values.data)))
    LinearProgrammingRoom(component, variable, component, parameter, sparse(A))
end

function grad_component(indices::Dict{Symbol, Any}, add_comp!::Function, variable::Symbol, parameter::Symbol)
    m = Model(Real)
    for index in keys(indices)
        set_dimension!(m, index, indices[index])
    end

    add_comp!(m, :Singleton)
    run(m) # Needs to be run for ModelInstance

    grad_component(m, :Singleton, variable, parameter)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Information about the condition of a component.

    IoK8sApiCoreV1ComponentCondition(;
        error=nothing,
        message=nothing,
        status=nothing,
        type=nothing,
    )

    - error::String : Condition error code for a component. For example, a health check error code.
    - message::String : Message about the condition for a component. For example, information about a health check.
    - status::String : Status of the condition for a component. Valid values for \"Healthy\": \"True\", \"False\", or \"Unknown\".
    - type::String : Type of condition for a component. Valid value: \"Healthy\"
"""
mutable struct IoK8sApiCoreV1ComponentCondition <: SwaggerModel
    error::Any # spec type: Union{ Nothing, String } # spec name: error
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiCoreV1ComponentCondition(;error=nothing, message=nothing, status=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ComponentCondition, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        validate_property(IoK8sApiCoreV1ComponentCondition, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1ComponentCondition, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(IoK8sApiCoreV1ComponentCondition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiCoreV1ComponentCondition

const _property_map_IoK8sApiCoreV1ComponentCondition = Dict{Symbol,Symbol}(Symbol("error")=>Symbol("error"), Symbol("message")=>Symbol("message"), Symbol("status")=>Symbol("status"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiCoreV1ComponentCondition = Dict{Symbol,String}(Symbol("error")=>"String", Symbol("message")=>"String", Symbol("status")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ComponentCondition }) = collect(keys(_property_map_IoK8sApiCoreV1ComponentCondition))
Swagger.property_type(::Type{ IoK8sApiCoreV1ComponentCondition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ComponentCondition[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ComponentCondition }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ComponentCondition[property_name]

function check_required(o::IoK8sApiCoreV1ComponentCondition)
    (getproperty(o, Symbol("status")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ComponentCondition }, name::Symbol, val)
end


text:
using NMF: nndsvd, solve!, CoordinateDescent

"""
    NMF(;ncomps=nothing)
    NMF(ncomps)

Use non-negative matrix factorization (NMF, NNMF) to form a non-negative, low-rank, and orthonormal basis of the input. The implementation of the underlying NMF is provided by [NMF.jl](https://github.com/JuliaStats/NMF.jl). The implementation uses a non-negative SVD for initialization and a coordinate-descent solver to fit.

If `ncomps` is `nothing`, it will be set to the number of frames in the reference cube when processed.

!!! warning "Non-negativity constraint"
    NMF is not designed to fit negative values. This algorithm will warn you (but will not error) if a target or reference cube contains negative values. The output may seem reasonable, but it is not well-defined with respect to the NMF algorithm. To overcome this, rescaling the data by its minimum before processing is necessary
    ```julia
    target = cube .- minimum(cube)
    S = reconstruct(NMF(), target, angles)
    ```
    When doing full-frame reduction (e.g. `NMF()(cube, angles)`) *this is handled automatically*, so this constraint only applies to the lower-level API and methods which rely on those, such as [`GreeDS`](@ref). **In general, if you see warnings, heed them.**

# References
1. [Ren et al. 2018](http://adsabs.harvard.edu/abs/2018ApJ...852..104R) Non-negative Matrix Factorization: Robust Extraction of Extended Structures
"""
struct NMF <: ADIAlgorithm
    ncomps
end
NMF(; ncomps=nothing) = NMF(ncomps)

function subtract(alg::NMF, cube; kwargs...)
    target = normalize_nmf_input(cube)
    if :ref in keys(kwargs)
        ref_ = normalize_nmf_input(kwargs[:ref])
        S = reconstruct(alg, target; kwargs..., ref=ref_)
    else
        S = reconstruct(alg, target; kwargs...)
    end
    return target .- S
end

function normalize_nmf_input(data)
    if any(<(0), data)
        return data .- minimum(data)
    end
    return data
end

function fit(alg::NMF, data::AbstractMatrix{T}; ref=data, kwargs...) where T
    if any(<(0), data) || any(<(0), ref)
        @warn "Negative values encountered in input. Make sure to rescale your data"
    end
    if alg.ncomps === nothing
        k = size(data, 1)
    else
        k = min(alg.ncomps, size(data, 1))
    end

    X = collect(ref)
    W, H = nndsvd(X, k)
    solve!(CoordinateDescent{T}(), X, W, H)
    if ref !== data
        Y = collect(data)
        W = Y * H'
        solve!(CoordinateDescent{T}(update_H=false), Y, W, H)
    end
    return LinearDesign(H, W)
end