text:
#
# utils.jl -
#
# Implements utility methods for Julia interface to ActiveSilicon Phoenix (PHX)
# library.
#
#------------------------------------------------------------------------------
#
# This file is part of the `Phoenix.jl` package which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut (https://github.com/emmt/Phoenix.jl).
# Copyright (C) 2016, Éric Thiébaut & Jonathan Léger.
#

function histogram(arr::AbstractArray{T,N}) where {T <: Integer, N}
    amin, amax = extrema(arr)
    cmin, cmax = Int(amin), Int(amax)
    off = 1 - cmin
    n = cmax + off
    y = zeros(Int, n)
    x = cmin:cmax
    @inbounds for v in arr
        i = off + Int(v)
        y[i] += 1
    end
    return (x, y)
end

const _DST_FORMATS = Dict(PHX_DST_FORMAT_Y8   => Monochrome{8},
                          PHX_DST_FORMAT_Y10  => Monochrome{10},
                          PHX_DST_FORMAT_Y12  => Monochrome{12},
                          # FIXME: PHX_DST_FORMAT_Y12B => 12,
                          PHX_DST_FORMAT_Y14  => Monochrome{14},
                          PHX_DST_FORMAT_Y16  => Monochrome{16},
                          PHX_DST_FORMAT_Y32  => Monochrome{32},
                          PHX_DST_FORMAT_Y36  => Monochrome{36},
                          PHX_DST_FORMAT_2Y12 => Monochrome{24},
                          PHX_DST_FORMAT_BAY8  => BayerFormat{8},
                          PHX_DST_FORMAT_BAY10 => BayerFormat{10},
                          PHX_DST_FORMAT_BAY12 => BayerFormat{12},
                          PHX_DST_FORMAT_BAY14 => BayerFormat{14},
                          PHX_DST_FORMAT_BAY16 => BayerFormat{16},
                          PHX_DST_FORMAT_RGB15 => RGB{15},
                          PHX_DST_FORMAT_RGB16 => RGB{16},
                          PHX_DST_FORMAT_RGB24 => RGB{24},
                          PHX_DST_FORMAT_RGB32 => RGB{32},
                          PHX_DST_FORMAT_RGB36 => RGB{36},
                          PHX_DST_FORMAT_RGB48 => RGB{48},
                          PHX_DST_FORMAT_BGR15 => BGR{15},
                          PHX_DST_FORMAT_BGR16 => BGR{16},
                          PHX_DST_FORMAT_BGR24 => BGR{24},
                          PHX_DST_FORMAT_BGR32 => BGR{32},
                          PHX_DST_FORMAT_BGR36 => BGR{36},
                          PHX_DST_FORMAT_BGR48 => BGR{48},
                          PHX_DST_FORMAT_BGRX32 => BGRX{32},
                          PHX_DST_FORMAT_RGBX32 => RGBX{32},
                          # FIXME: PHX_DST_FORMAT_RRGGBB8 => 8,
                          PHX_DST_FORMAT_XBGR32 => XBGR{32},
                          PHX_DST_FORMAT_XRGB32 => XRGB{32},
                          PHX_DST_FORMAT_YUV422 => YUV422)

const CAPTURE_FORMATS = Union{Monochrome{8},
                              Monochrome{10},
                              Monochrome{12},
                              Monochrome{14},
                              Monochrome{16},
                              Monochrome{32},
                              Monochrome{36},
                              Monochrome{24},
                              BayerFormat{8},
                              BayerFormat{10},
                              BayerFormat{12},
                              BayerFormat{14},
                              BayerFormat{16},
                              RGB{15},
                              RGB{16},
                              RGB{24},
                              RGB{32},
                              RGB{36},
                              RGB{48},
                              BGR{15},
                              BGR{16},
                              BGR{24},
                              BGR{32},
                              BGR{36},
                              BGR{48},
                              BGRX{32},
                              RGBX{32},
                              XBGR{32},
                              XRGB{32},
                              YUV422}

"""

    capture_format(fmt) -> pix

yields the pixel format corresponding to the destination buffer pixel format
`fmt`, and conversely:

    capture_format(pix) -> fmt

For instance:

    capture_format(PHX_DST_FORMAT_Y8) -> Monochrome{8}
    capture_format(Monochrome{8}) -> PHX_DST_FORMAT_Y8

See also: [`best_capture_format`](@ref), [`capture_format_bits`](@ref).

"""
function capture_format(fmt::Integer)
    if haskey(_DST_FORMATS, fmt)
        return getindex(_DST_FORMATS, fmt)
    end
    throw(ArgumentError("unknown capture pixel format"))
end

capture_format(::Type{T}) where {T<:Monochrome{8}} = PHX_DST_FORMAT_Y8
capture_format(::Type{T}) where {T<:Monochrome{10}} = PHX_DST_FORMAT_Y10
capture_format(::Type{T}) where {T<:Monochrome{12}} = PHX_DST_FORMAT_Y12
# FIXME: PHX_DST_FORMAT_Y12B
capture_format(::Type{T}) where {T<:Monochrome{14}} = PHX_DST_FORMAT_Y14
capture_format(::Type{T}) where {T<:Monochrome{16}} = PHX_DST_FORMAT_Y16
capture_format(::Type{T}) where {T<:Monochrome{32}} = PHX_DST_FORMAT_Y32
capture_format(::Type{T}) where {T<:Monochrome{36}} = PHX_DST_FORMAT_Y36
capture_format(::Type{T}) where {T<:Monochrome{24}} = PHX_DST_FORMAT_2Y12
capture_format(::Type{T}) where {T<:BayerFormat{8}} = PHX_DST_FORMAT_BAY8
capture_format(::Type{T}) where {T<:BayerFormat{10}} = PHX_DST_FORMAT_BAY10
capture_format(::Type{T}) where {T<:BayerFormat{12}} = PHX_DST_FORMAT_BAY12
capture_format(::Type{T}) where {T<:BayerFormat{14}} = PHX_DST_FORMAT_BAY14
capture_format(::Type{T}) where {T<:BayerFormat{16}} = PHX_DST_FORMAT_BAY16
capture_format(::Type{T}) where {T<:RGB{15}} = PHX_DST_FORMAT_RGB15
capture_format(::Type{T}) where {T<:RGB{16}} = PHX_DST_FORMAT_RGB16
capture_format(::Type{T}) where {T<:RGB{24}} = PHX_DST_FORMAT_RGB24
capture_format(::Type{T}) where {T<:RGB{32}} = PHX_DST_FORMAT_RGB32
capture_format(::Type{T}) where {T<:RGB{36}} = PHX_DST_FORMAT_RGB36
capture_format(::Type{T}) where {T<:RGB{48}} = PHX_DST_FORMAT_RGB48
capture_format(::Type{T}) where {T<:BGR{15}} = PHX_DST_FORMAT_BGR15
capture_format(::Type{T}) where {T<:BGR{16}} = PHX_DST_FORMAT_BGR16
capture_format(::Type{T}) where {T<:BGR{24}} = PHX_DST_FORMAT_BGR24
capture_format(::Type{T}) where {T<:BGR{32}} = PHX_DST_FORMAT_BGR32
capture_format(::Type{T}) where {T<:BGR{36}} = PHX_DST_FORMAT_BGR36
capture_format(::Type{T}) where {T<:BGR{48}} = PHX_DST_FORMAT_BGR48
capture_format(::Type{T}) where {T<:BGRX{32}} = PHX_DST_FORMAT_BGRX32
capture_format(::Type{T}) where {T<:RGBX{32}} = PHX_DST_FORMAT_RGBX32
# FIXME: PHX_DST_FORMAT_RRGGBB8
capture_format(::Type{T}) where {T<:XBGR{32}} = PHX_DST_FORMAT_XBGR32
capture_format(::Type{T}) where {T<:XRGB{32}} = PHX_DST_FORMAT_XRGB32
capture_format(::Type{T}) where {T<:YUV422} = PHX_DST_FORMAT_YUV422
capture_format(::Type{T}) where {T<:PixelFormat} = ParamValue(0)

"""

`capture_format_bits(fmt)` yields the number of bits per pixel for the
destination buffer pixel format `fmt`, *e.g.* `PHX_DST_FORMAT_Y8`.

See also: [`best_capture_format`](@ref).

"""
capture_format_bits(fmt::Integer) :: Int =
    (haskey(_DST_FORMATS, fmt) ?
     bitsperpixel(getindex(_DST_FORMATS, fmt)) : -1)

"""

    best_capture_format(cam) -> (T, pixelformat)

yields the Julia pixel type `T` and Phoenix pixel format `pixelformat` for the
destination image buffer(s) which are as close as possible as that of the
source image buffer(s) for camera `cam`.  The value of `T` can be used to
allocate a Julia array while `pixelformat` can be used to set
`cam[PHX_DST_FORMAT]`.

The method can also be called with the color format (`cam[PHX_CAM_SRC_COL]`)
and the number of bits per pixel (`cam[PHX_CAM_SRC_DEPTH]`) of the camera:

    best_capture_format(color, depth) -> (T, pixelformat)

See also: [`capture_format_bits`](@ref).

"""
best_capture_format(cam::Camera) =
    best_capture_format(cam[PHX_CAM_SRC_COL], cam[PHX_CAM_SRC_DEPTH])

function best_capture_format(color::Integer, depth::Integer)
    local T::DataType, format::ParamValue
    if color == PHX_CAM_SRC_MONO
        if depth ≤ 8
            T, format = UInt8, PHX_DST_FORMAT_Y8
        elseif depth ≤ 16
            T, format = UInt16, PHX_DST_FORMAT_Y16
        elseif depth ≤ 32
            T, format = UInt32, PHX_DST_FORMAT_Y32
        else
            error("Unsupported monochrome camera depth $(Int(depth))")
        end
    elseif color == PHX_CAM_SRC_YUV422
        T, format = YUV422BitsType, PHX_DST_FORMAT_YUV422
    elseif color == PHX_CAM_SRC_RGB
        if depth == 8
            T, format = RGB24BitsType, PHX_DST_FORMAT_RGB24
        elseif depth == 16
            T, format = RGB48BitsType, PHX_DST_FORMAT_RGB48
        else
            error("Unsupported RGB camera depth $(Int(depth))")
        end
    elseif (color == PHX_CAM_SRC_BAY_RGGB || color == PHX_CAM_SRC_BAY_GRBG ||
            color == PHX_CAM_SRC_BAY_GBRG || color == PHX_CAM_SRC_BAY_BGGR)
        if depth ≤ 8
            T, format = UInt8, PHX_DST_FORMAT_BAY8
        elseif depth ≤ 16
            T, format = UInt16, PHX_DST_FORMAT_BAY16
        else
            error("Don't know how to interpret Bayer color format (0x$(string(color, base=16)))")
        end
    else
        error("Unknown camera color format (0x$(string(color, base=16)))")
    end
    @assert sizeof(T)*8 == capture_format_bits(format)
    return T, format
end

"""

    cstring(str, len = strlen(str))

yields a vector of bytes (`UInt8`) with the contents of the string `str` and
properly zero-terminated.  This buffer is independent from the input string and
its contents can be overwritten.  An error is thrown if `str` contains any
embedded NUL characters (which would cause the string to be silently truncated
if the C routine treats NUL as the terminator).

An alternative (without the checking of embedded NUL characters) is:

    push!(convert(Vector{UInt8}, str), convert(UInt8, 0))

"""
function cstring(str::AbstractString,
                 len::Integer = length(str)) :: Array{UInt8}
    buf = Array{UInt8}(undef, len + 1)
    i = 0
    @inbounds for c in str
        if i ≥ len
            break
        end
        c != '\0' || error("string must not have embedded NUL characters")
        i += 1
        buf[i] = c
    end
    @inbounds while i ≤ len
        i += 1
        buf[i] = zero(UInt8)
    end
    return buf
end

function Base.summary(cam::Camera)

    # Get hardware revision.
    @printf("Hardware revision:       %.2x:%.2x:%.2x\n",
            cam[PHX_REV_HW_MAJOR],
            cam[PHX_REV_HW_MINOR],
            cam[PHX_REV_HW_SUBMINOR])

    # Get software revision.
    @printf("Software revision:       %.2x:%.2x:%.2x\n",
            cam[PHX_REV_SW_MAJOR],
            cam[PHX_REV_SW_MINOR],
            cam[PHX_REV_SW_SUBMINOR])

    # Get the board properties.
    println("Board properties:")
    for str in split(cam[PHX_BOARD_PROPERTIES], '\n', keepempty = false)
        println("    ", str)
    end

    # CoaXPress camera?
    @printf("CoaXPress camera:         %s\n",
            (iscoaxpress(cam) ? "true" : "false"))

    @printf("Camera active xoffset:    %4d\n",
            Int(cam[PHX_CAM_ACTIVE_XOFFSET]))
    @printf("Camera active yoffset:    %4d\n",
            Int(cam[PHX_CAM_ACTIVE_YOFFSET]))
    @printf("Camera active xlength:    %4d\n",
            Int(cam[PHX_CAM_ACTIVE_XLENGTH]))
    @printf("Camera active ylength:    %4d\n",
            Int(cam[PHX_CAM_ACTIVE_YLENGTH]))

end

"""

`iscoaxpress(cam)` yields whether Phoenix camara `cam` is a CoaXPress
controlled camera.

See also: [`Phoenix.Camera`](@ref), [`Phoenix.assert_coaxpress`](@ref).

"""
iscoaxpress(cam::Camera) = false
iscoaxpress(cam::CoaXPressCamera) = true

"""

`assert_coaxpress(cam)` checks whether Phoenix camara `cam` is a CoaXPress
controlled camera and throw an error otherwise.

See also: [`Phoenix.iscoaxpress`](@ref).

"""
function assert_coaxpress(cam::Camera)
    if iscoaxpress(cam)
        magic = cam[CXP_STANDARD]
        if magic != CXP_MAGIC
            if magic == bswap(CXP_MAGIC)
                error("bad byte swapping flag")
            else
                error("unexpected magic number for CoaXPress camera")
            end
        end
    else
        error("not a CoaXPress camera")
    end
    return nothing
end

getvendorname(cam::Camera) =
    iscoaxpress(cam) ? cam[CXP_DEVICE_VENDOR_NAME] : ""

getmodelname(cam::Camera) =
    iscoaxpress(cam) ? cam[CXP_DEVICE_MODEL_NAME] : ""

getdevicemanufacturer(cam::Camera) =
    iscoaxpress(cam) ? cam[CXP_DEVICE_MANUFACTURER_INFO] : ""

getdeviceversion(cam::Camera) =
    iscoaxpress(cam) ? cam[CXP_DEVICE_VERSION] : ""

getdeviceserialnumber(cam::Camera) =
    iscoaxpress(cam) ? cam[CXP_DEVICE_SERIAL_NUMBER] : ""

getdeviceuserid(cam::Camera) =
    iscoaxpress(cam) ? cam[CXP_DEVICE_USER_ID] : ""

subsampling_parameter(sub::Integer) =
    (sub == 1 ? PHX_ACQ_X1 :
     sub == 2 ? PHX_ACQ_X2 :
     sub == 4 ? PHX_ACQ_X4 :
     sub == 8 ? PHX_ACQ_X8 :
     error("unsupported horizontal subsampling ratio"))

"""

    gettimeofday()

yields the current time since the epoch as a `TimeVal` structure.
To convert it in seconds, just convert the result into a float:

    float(gettimeofday())

which yields the is the same as

"""
function gettimeofday()
    tvref = Ref{TimeVal}()
    status = ccall(:gettimeofday, Cint, (Ptr{TimeVal}, Ptr{Cvoid}),
                   tvref, C_NULL)
    status == SUCCESS || error("gettimeofday failed")
    return tvref[]
end

TimeVal(sec::Real) = TimeVal(Float64(sec))
TimeVal(sec::Integer) = TimeVal(sec, 0)
TimeVal(tv::TimeVal) = tv
function TimeVal(ts::TimeSpec)
    µs, ns = divrem(ts.nsec, 1_000)
    s, µs = divrem(µs + (ns ≥ 500 ? 1 : 0), 1_000_000)
    return TimeVal(ts.sec + s, µs)
end
function TimeVal(sec::Float64)
    if isfinite(sec)
        s = trunc(Int, sec)
        µs = round(Int, (sec - s)*1_000_000)
        if µs ≥ 1_000_000
            d, µs = divrem(µs, 1_000_000)
            s += d
        end
        return TimeVal(s, µs)
    elseif isinf(sec)
        return TimeVal(sec == Inf ? typemax(_typeof_tv_sec) :
                        typemin(_typeof_tv_sec), 0)
    else
        throw(DomainError())
    end
end

"""

    TimeSpec(s [, ns])

yields a `TimeSpec` value for `s` seconds and `ns` nanoseconds.  If both `s`
and `ns` are specified, they must be integres; otherwise, `s` can be a
fractional number of seconds.  For instance:

    TimeSpec(time() + 10.3)

yields a `TimeSpec` instance for `10.3` seconds later.  Use the `float` method
to convert a `TimeSpec` value in a floating-point value in seconds.

"""
TimeSpec(sec::Real) = TimeSpec(Float64(sec))
TimeSpec(sec::Integer) = TimeSpec(sec, 0)
TimeSpec(tv::TimeVal) = TimeSpec(tv.sec, tv.usec*1_000)
TimeSpec(ts::TimeSpec) = ts
function TimeSpec(sec::Float64)
    if isfinite(sec)
        s = trunc(Int, sec)
        ns = round(Int, (sec - s)*1_000_000_000)
        if ns ≥ 1_000_000_000
            d, ns = divrem(ns, 1_000_000_000)
            s += d
        end
        return TimeSpec(s, ns)
    elseif isinf(sec)
        return TimeSpec(sec == Inf ? typemax(_typeof_tv_sec) :
                        typemin(_typeof_tv_sec), 0)
    else
        throw(DomainError())
    end
end

#Base.convert(::Type{TimeVal}, x) = TimeVal(x)
for T in (Float32, Float64)
    @eval begin
        Base.$T(tv::TimeVal)  = ($T(tv.sec) + $T(tv.usec)*$T(1E-6))
        Base.$T(ts::TimeSpec) = ($T(ts.sec) + $T(ts.nsec)*$T(1E-9))
        #Base.convert(::Type{$T}, tv::TimeVal) = $T(tv)
        #Base.convert(::Type{$T}, ts::TimeSpec) = $T(ts)
    end
end
Base.float(tv::TimeVal) = Float64(tv)
Base.float(ts::TimeSpec) = Float64(ts)

#function Base.(+)(ts::TimeSpec, tv::TimeVal)
#    sec, nsec = divrem(ts.nsec + 1_000*tv.usec, 1_000_000_000)
#    sec += ts.sec + tv.sec
#    return TimeSpec(sec, nsec)
#end
#
#Base.(+)(tv::TimeVal, ts::TimeSpec) = ts + tv

"""
    isforever(ts)

yeilds whether the time specified `ts` should be consider as being forever,
that is in a quasi-infinite future (due to the finite precision, this is
more than 292.3 billions years).

"""
isforever(ts::TimeSpec) = (ts.sec ≥ typemax(_typeof_tv_sec))

"""
    rounddown(a, b)

yields largest multiple of `b` which is smaller or equal `a`.  Arguments and
result are integers, `a` must be nonnegative, `b` must be striclty positive and
result is nonnegative.

See also: [`roundup`](@ref)

"""
rounddown(a::Integer, b::Integer) = div(a, b)*b

"""
    roundup(a, b)

yields smallest multiple of `b` which is larger or equal `a`.  Arguments and
result are integers, `a` must be nonnegative, `b` must be striclty positive and
result is nonnegative.

See also: [`rounddown`](@ref)

"""
roundup(a::Integer, b::Integer) = rounddown((b - one(b)) + a, b)

"""
     divrnd(a, b) -> q

yields `a/b` rounded to the nearest integer.  Arguments are assumed to be
positive integers.  The result is the same as `round(Int, a/b)` but (slightly)
faster.

"""
divrnd(a::Integer, b::Integer) = div(2a + b, 2b)


text:
function _SNumber(name)
    m = match(r"(?<=_S)\d+$", name)
    m != nothing ? m.match : ""
end

# Creates DataFrame with columns (SampleID, Run, Reference, S (if applicable), Rejected, Comment)
# parent - Bam folder
# runName - Name of run to create template for
# logFile - "AlignUtils.log", used to find reference used for each sample
function metadatatemplate(syn, runName::AbstractString; bamFolder::AbstractString="", logFile::AbstractString="", exclude::Pattern=r"Undetermined")
    isempty(bamFolder) && error("Bam folder must be specified.")

    # 1. List aligned samples
    # 2. Ignored excluded
    # 3. Create table with columns
    #   a. SampleID
    #   b. Run
    #   c. Reference (use Provenance?)
    #   d. S number (if applicable)
    #   e. Rejected
    #   f. Comment

    ids,names = listaligned(syn, bamFolder, runName::AbstractString)
    mask = .~matchany(exclude, names)
    ids,names = ids[mask],names[mask]
    map!(x->x[1:end-4], names, names) # remove ".bam" from end

    references = isempty(logFile) ? "" : referencefromlog(syn, logFile, names)

    # extract S numbers
    SNumbers = map(_SNumber,names)

    DataFrame(SampleID=names, Run=runName, Reference=references, S=SNumbers, Rejected="", Comment="")
end

function metadatatemplate(syn, alignmentFolder::AbstractString, runName::AbstractString; exclude::Pattern=r"Undetermined", logFileName="AlignUtils.log")
    bamFolder = getchildbyname(syn, alignmentFolder, "Bam")
    scriptsFolder = getchildbyname(syn, alignmentFolder, "Scripts")
    scriptsRunFolder = getchildbyname(syn, scriptsFolder, runName)
    logFile = getchildbyname(syn, scriptsRunFolder, logFileName)

    metadatatemplate(syn,runName,bamFolder=bamFolder,logFile=logFile,exclude=exclude)
end

# filtermetadata(metadata, filter1, filter2, ...)
# where filter is one of the following
#        :fieldname                - Check if the fieldname is a column in the metadata DataFrame. Otherwise return empty DataFrame.
#        :fieldname=>value         - Only keep samples that has the given value for the given field.
#        :fieldname=>[v1,v2...,vn] - Only keep samples that one of the given values for the given field.
#        :fieldname=>Predicate     - Only keep samples where the predicate function evaluates to true for the value of the given field name.
#        Predicate                 - Only keep samples for which Predicate returns true. Predicate acts on entire sample row.
filtermetadata(metadata::DataFrame, filter::Symbol) = filter in names(metadata) ? metadata : DataFrame()
function filtermetadata(metadata::DataFrame, filter::Pair{Symbol,T}) where {T<:AbstractArray}
    filter.first in names(metadata) || return DataFrame()
    mask = map!(x->!ismissing(x) && x in filter.second, falses(size(metadata,1)), metadata[!,filter.first])
    count(mask)==0 ? DataFrame() : metadata[mask,:]
end
function filtermetadata(metadata::DataFrame, filter::Pair{Symbol,T}) where {T}
    filter.first in names(metadata) || return DataFrame()
    mask = .~ismissing.(metadata[!,filter.first]) .& (metadata[!,filter.first] .== filter.second)
    count(mask)==0 ? DataFrame() : metadata[mask,:]
end
function filtermetadata(metadata::DataFrame, filter::Pair{Symbol,Function})
    filter.first in names(metadata) || return DataFrame()
    mask = map!(filter.second, falses(size(metadata,1)), metadata[!,filter.first])
    count(mask)==0 ? DataFrame() : metadata[mask,:]
end
function filtermetadata(metadata::DataFrame, filter::Function)
    N = size(metadata,1)
    mask = falses(N)
    for i=1:N
        mask[i] = filter(metadata[i,:])
    end
    count(mask)==0 ? DataFrame() : metadata[mask,:]
end

function filtermetadata(metadata, args...)
    for arg in args
        metadata = filtermetadata(metadata,arg)
    end
    metadata
end

function _loadmetadata(syn, fileID)
    T = CSV.read(localpath(syn, fileID); missingstrings=["","NA"])
    T[!,:MetadataID] .= fileID # make an extra column with a reference to the metadata file!
    T
end

getsamplemetadata(syn, file::File, args...) = filtermetadata(_loadmetadata(syn, file.id), args...)
function getsamplemetadata(syn, folder::Folder, args...)
    #filtermetadata(_loadmetadata(syn, metadataID), args...)

    paths, names = listfiles(syn, folder.id) # TODO: update listfiles to support arguments of type Folder
    mask = map(x->lowercase(splitext(x)[2])==".csv", names)
    paths, names = paths[mask], names[mask]

    allMetadata = Vector{DataFrame}(length(paths))
    map!(path->getsamplemetadata(syn,path,args...), allMetadata, paths)
    vcat(allMetadata...)
end

# metadataID: Synapse id, could be file or folder (should point to MyProject/Metadata).
#             If it's a folder, all .csv files in the folder will be checked.
#             filtering specified by args... (see filtermetadata())
getsamplemetadata(syn, metadataID::AbstractString, args...) = getsamplemetadata(syn, get(syn, metadataID, downloadFile=false), args...)

function appendsynapseids!(syn, metadata::DataFrame, folderID::AbstractString, fileSuffixes::AbstractArray, columnNames::AbstractArray{Symbol}; subFolder="")
    @assert length(fileSuffixes)==length(columnNames)

    for columnName in columnNames
        metadata[!,columnName] .= ""
    end

    # partition samples by run, to make as few Synapse queries as possible
    runNames = metadata[!,:Run]
    for runName in unique(runNames)
        runInd = findall(runNames.==runName) # index of samples in this run in metadata table
        sampleNames = metadata[runInd,:SampleID]

        sampleFolder = getchildbyname(syn, folderID, runName)
        isempty(subFolder) || (sampleFolder = getchildbyname(syn, sampleFolder, subFolder))

        fileIDs,filenames = listfiles(syn, sampleFolder)

        for (suffix,columnName) in zip(fileSuffixes,columnNames)
            fileMask = falses(length(filenames))
            map!( x->endswith(x,suffix), fileMask, filenames ) # find matches for suffix

            matchingIDs   = fileIDs[fileMask] # Synapse ids for files with matching suffix
            matchingNames = map(x->x[1:end-length(suffix)], filenames[fileMask]) # get rid of suffix from name

            matchingInd = indexin(sampleNames, matchingNames)

            matchingMask = matchingInd.!=nothing
            metaInd = runInd[matchingMask]
            matchingInd = matchingInd[matchingMask]

            metadata[metaInd,columnName] = matchingIDs[matchingInd]
        end

    end

    metadata
end
appendsynapseids!(syn, metadata::DataFrame, folderID::AbstractString, fileSuffix::AbstractString, columnName::Symbol; kwargs...) = appendsynapseids!(syn, metadata, folderID, [fileSuffix], [columnName]; kwargs...)

# appendswarmids!(syn, metadata::DataFrame, alignmentFolderID::AbstractString) = appendsynapseids!(syn, metadata, getchildbyname(syn,alignmentFolderID,"MutantSwarms"),".jld",:SwarmID)

# strands is one of: :both, :forward, :reverse, or an array of these
function appendswarmids!(syn, metadata::DataFrame, alignmentFolderID::AbstractString, strand::Symbol=:both)
    @assert strand in [:both, :forward, :reverse]
    if strand==:both
        appendsynapseids!(syn, metadata, getchildbyname(syn,alignmentFolderID,"MutantSwarms"), ".jld", :SwarmID)
    elseif strand==:forward
        appendsynapseids!(syn, metadata, getchildbyname(syn,alignmentFolderID,"MutantSwarms"), ".jld", :SwarmForwardID, subFolder="forward")
    else#if strand==:reverse
        appendsynapseids!(syn, metadata, getchildbyname(syn,alignmentFolderID,"MutantSwarms"), ".jld", :SwarmReverseID, subFolder="reverse")
    end
end
function appendswarmids!(syn, metadata::DataFrame, alignmentFolderID::AbstractString, strands::AbstractArray)
    for s in strands
        metadata = appendswarmids!(syn, metadata::DataFrame, alignmentFolderID::AbstractString, s)
    end
    metadata
end

function appendalignedids!(syn, metadata::DataFrame, alignmentFolderID::AbstractString)
    appendsynapseids!(syn, metadata, getchildbyname(syn,alignmentFolderID,"Bam"),
                      [".bam",".bam.bai","_consensus.fasta"],
                      [:BamID,:BaiID,:ConsensusID])
end

function downloadbymeta!(syn, metadata::DataFrame, IDColumn::Symbol, pathColumn::Symbol; cacheDir="")
    @assert hasproperty(metadata,IDColumn)

    kwargs = isempty(cacheDir) ? () :

    ids = metadata[!,IDColumn]
    mask = .~ismissing.(ids)

    paths = isempty(cacheDir) ?
            localpath(syn, ids[mask]) :
            localpath(syn, ids[mask], downloadLocation=cacheDir, ifcollision="overwrite.local")

    metadata[!,pathColumn] .= ""
    metadata[mask,pathColumn] = paths
    metadata
end

#downloadswarms!(syn, metadata::DataFrame; kwargs...) = downloadbymeta!(syn, metadata, :SwarmID, :SwarmPath; kwargs...)

# strands is one of: :both, :forward, :reverse, or an array of these
function downloadswarms!(syn, metadata::DataFrame, strand::Symbol=:both; cacheDir="", kwargs...)
    @assert strand in [:both, :forward, :reverse]
    if strand==:both
        downloadbymeta!(syn, metadata, :SwarmID, :SwarmPath; cacheDir=cacheDir, kwargs...)
    elseif strand==:forward
        isempty(cacheDir) || (cacheDir=joinpath(cacheDir,"forward"))
        downloadbymeta!(syn, metadata, :SwarmForwardID, :SwarmForwardPath; cacheDir=cacheDir, kwargs...)
    else#if strand==:reverse
        isempty(cacheDir) || (cacheDir=joinpath(cacheDir,"reverse"))
        downloadbymeta!(syn, metadata, :SwarmReverseID, :SwarmReversePath; cacheDir=cacheDir, kwargs...)
    end
end
function downloadswarms!(syn, metadata::DataFrame, strands::AbstractArray; kwargs...)
    for s in strands
        metadata = downloadswarms!(syn, metadata, s; kwargs...)
    end
    metadata
end

downloadconsensuses!(syn, metadata::DataFrame; kwargs...) = downloadbymeta!(syn, metadata, :ConsensusID, :ConsensusPath; kwargs...)
downloadbams!(syn, metadata::DataFrame; kwargs...) = downloadbymeta!(syn, metadata, :BamID, :BamPath; kwargs...)
downloadbais!(syn, metadata::DataFrame; kwargs...) = downloadbymeta!(syn, metadata, :BaiID, :BaiPath; kwargs...)
function downloadaligned!(syn, metadata::DataFrame; kwargs...)
    downloadbams!(syn, metadata; kwargs...)
    downloadbais!(syn, metadata; kwargs...)
    downloadconsensuses!(syn, metadata; kwargs...)
end

ext: jl


text:
module TestTwoInertiasAndIdealGearWithUnitsAndUncertainties

using Modia
@usingModiaPlot
using Modia.Measurements

TwoInertiasAndIdealGear = Model(
    J1 = 0.0025u"kg*m^2",
    J2 = (150.0 ± 20.0)u"kg*m^2",
    r  = 105.0,
    tau_max = 1u"N*m",

    phi2 = Var(start = (0.5 ± 0.05)u"rad"),
    w2   = Var(start = 0.0u"rad/s"),
    tau2 = Var(start = 0u"N*m"),

    equations = :[
        tau = if time < 1u"s"; tau_max elseif time < 2u"s"; 0u"N*m" elseif time < 3u"s"; -tau_max else 0u"N*m" end,

        # inertia1
        w1 = der(phi1),
        J1*der(w1) = tau - tau1,

        # gear
        phi1   = r*phi2,
        r*tau1 = tau2,

        # inertia2
        w2 = der(phi2),
        J2*der(w2) = tau2
    ]
)

twoInertiasAndIdealGear = @instantiateModel(TwoInertiasAndIdealGear, FloatType = Measurement{Float64})

simulate!(twoInertiasAndIdealGear, Tsit5(), stopTime = 4.0, logParameters=true, logStates=true)

plot(twoInertiasAndIdealGear, ["phi2", "w2", "der(w2)"])

# Linearize
println("\n... Analytic linearization")
(A1, x1) = linearize!(twoInertiasAndIdealGear, stopTime=4, analytic = true)
xNames = get_xNames(twoInertiasAndIdealGear)
@show xNames
println(IOContext(stdout, :error_digits=>15), "A1 = ", A1, ", x1 = ", x1)

println("\n... Numeric linearization with Float64")
(A2, x2) = linearize!(twoInertiasAndIdealGear, stopTime=4, analytic=false)
println(IOContext(stdout, :error_digits=>15), "A2 = ", A2, ", x2 = ", x2)

#= DoubleFloats not defined
using Test
println("\n... Numeric linearization with Double64")
using DoubleFloats
twoInertiasAndIdealGear2 = SimulationModel{Measurement{Double64}}(twoInertiasAndIdealGear)
(A3, x3) = linearize!(twoInertiasAndIdealGear2, stopTime=3, analytic=false)
println(IOContext(stdout, :error_digits=>15), "A3 = ", A3, ", x3 = ", x3)
=#

end


text:
"""
    flatten_dir(path::String)

Takes all files nested in a directory structure and takes them to the root path.
Warning: this effects your file system. Be careful.
"""
function flatten_dir(path::String)
    allfiles = [ joinpath.(root, files) for (root, _, files) in walkdir(path) ]
    file_count = 0
    for src in Iterators.flatten( allfiles )
        if dirname(src) != path
            Base.Filesystem.mv(src, Base.Filesystem.joinpath(path, basename(src) ))
            file_count += 1
        end
    end
    @warn   "ChemometricsData.jl has flattened $path containing $file_count files. \n" *
            "The directory structure has been preserved."
    return nothing
end

"""
    isnumeric(s::String)::Boolean

returns true if the String `s` contains only numeric characters.
"""
isnumeric(s::String)::Boolean = tryparse(Float64, s) isa Number

"""
    numeric_column_names( df::DataFrame )

Returns a subset of a dataframe whose column names can be interpretted as numeric values.
"""
function numeric_column_names( df::DataFrame )
    cols = String.( names( df ) )
    numer_cols = Symbol.( cols[ isnumeric.( cols ) ] )
    return df[!, numer_cols]
end

"""
    nonnumeric_column_names( df::DataFrame )

Returns a subset of a dataframe whose column names can not be interpretted as numeric values.
"""
function nonnumeric_column_names( df::DataFrame )
    cols = String.( names( df ) )
    numer_cols = Symbol.( cols[ .!isnumeric.( cols ) ] )
    return df[!, numer_cols]
end


text:
function swsetupmatrix(s1, s2)
    setupmatrix = zeros(Int, length(s1)+1, length(s2)+1)
    return setupmatrix
end

function swscorematrix(seq1, seq2; match=1, mismatch=-1, gap=-1)
    scoremat = swsetupmatrix(seq1, seq2)
    for i in 2:size(scoremat, 1) # iterate through row indices
        for j in 2:size(scoremat, 2) # iterate through column indices
            above = nwscore(seq1[i-1], nothing; match, mismatch, gap) + scoremat[i-1,j]
            left = nwscore(nothing, seq2[j-1]; match, mismatch, gap) + scoremat[i,j-1]
            diagonal = nwscore(seq1[i-1], seq2[j-1]; match, mismatch, gap) + scoremat[i-1,j-1]
            scoremat[i,j] = max(above, left, diagonal)
            # or `scoremat[i,j] = max(above, left, diagonal, 0)`

            if scoremat[i,j] < 0
                scoremat[i,j] = 0
            end
        end
    end
    return scoremat
end

"""
    swalign(seq1::String, seq2::String; match =1, mismatch=-1, gap=-1)

    This function takes the results of a matrix scored using smith waterman properties,
    then returns the proper alignment of sequences based on the matrix as a tuple of two strings.

    Ex.)
    julia> swalign("AAACCCGGG","TTTCCCAAA")
    ("CCC", "CCC")

    julia> swalign("TGTTACGG", "GGTTGACTA")
    ("GTT-AC", "GTTGAC")
"""
function swalign(seq1::String, seq2::String; match =1, mismatch=-1, gap=-1)
    if !isDNA(seq1) || !isDNA(seq2) || length(seq1) == 0 || length(seq2) == 0
        throw(ErrorException("Invalid DNA sequence provided"))
    end

    scoremat = swscorematrix(seq1, seq2; match=match, mismatch=mismatch, gap=gap)
    aligned1 = ""
    aligned2 = ""

    maxScore = findmax(scoremat)

    i = maxScore[2][1]
    j = maxScore[2][2]
    #@info i
    #@info j

    while i > 1 && j > 1
            if scoremat[i,j] == 0
                break
            #check scores from left, right, diag, and if any match current cell, then add that char to the alignments
            elseif scoremat[i,j] == nwscore(seq1[i-1], nothing; match, mismatch, gap) + scoremat[i-1,j]
                #@info "max is from above"
                aligned1 = aligned1 * seq1[i-1]
                aligned2 = aligned2 * "-"
                #update the traversal indices
                i = i-1
            elseif scoremat[i,j] == nwscore(nothing, seq2[j-1]; match, mismatch, gap) + scoremat[i,j-1]
                #@info "max is from left"
                aligned1 = aligned1 * "-"
                aligned2 = aligned2 * seq2[j-1]
                #update the traversal indices
                j = j-1
            else
                #@info "max is from diag"
                aligned1 = aligned1 * seq1[i-1]
                aligned2 = aligned2 * seq2[j-1]
                #update the traversal indices
                i = i-1
                j = j-1
            end

    end
    aligned1 = reverse(aligned1)
    aligned2 = reverse(aligned2)
    return (aligned1, aligned2)
end

# Labs 3 & 4

"""
    nwscore(base1, base2; mismatch = -1, match = 1, gap = -1 )

Your documentation here.

This function compares nucleotides for genomic sequences
to see if they are matching=1, mismatching=-1, or gaps=-1.

Ex.)
```
julia> nwscore('A','G')
-1
julia> nwscore('A','A')
1
julia> nwscore('A',nothing)
-1
julia> nwscore(nothing,nothing)
ERROR: ArgumentError: Score for two gaps is not defined
Stacktrace:
[1] nwscore(#unused#::Nothing, #unused#::Nothing)
    @ Main ~/Documents/my_repo/BISC195Labs/src/needleman_wunch.jl:25
[2] top-level scope
    @ REPL[81]:1
```

"""
function nwscore(base1::Char, base2::Char; match = 1, mismatch = -1, gap = -1)
    if base1 == base2
        return match
    else base1 != base2
        return mismatch
    end
end

function nwscore(base::Char, ::Nothing; match = 1, mismatch = -1, gap = -1)
    return gap
end

function nwscore(::Nothing, base::Char; match = 1, mismatch = -1, gap = -1)
    nwscore(base, nothing; gap=gap)
end

function nwscore(::Nothing, ::Nothing)
    throw(ArgumentError("Score for two gaps is not defined"))
end

# function nwalign(seq1::String, seq2::String; match =1, mismatch=-1, gap=-1)
#     #input sequence 1 (seq1)
#     #input sequence 2 (seq2)
#     #turn each sequence into an array/vector of characters
#     #compare each character based on position (see grid for rules) using nwscore
#     #make sure there is error if 2 nothings
#     #determine best alignment depending on score
#         #if same score, print both as touple
#     #return alignments (seq1, seq2) as an array
# end

function nwsetupmatrix(s1, s2; gap=-1)
    setupmatrix = zeros(Int, length(s1)+1, length(s2)+1)
    for j in 1:length(s2)+1
        setupmatrix[1, j] = (j-1)*gap
    end
    for i in 1:length(s1)+1
        setupmatrix[i,1] = (i-1)*gap
    end
    return setupmatrix
end

function nwscorematrix(seq1, seq2; match=1, mismatch=-1, gap=-1)
    scoremat = nwsetupmatrix(seq1, seq2; gap=gap)
    for i in 2:size(scoremat, 1) # iterate through row indices
        for j in 2:size(scoremat, 2) # iterate through column indices
            above = nwscore(seq1[i-1], nothing; match, mismatch, gap) + scoremat[i-1,j]
            left = nwscore(nothing, seq2[j-1]; match, mismatch, gap) + scoremat[i,j-1]
            diagonal = nwscore(seq1[i-1], seq2[j-1]; match, mismatch, gap) + scoremat[i-1,j-1]
            scoremat[i,j] = max(above, left, diagonal)
        end
    end
    return scoremat
end

#need to write a function that finds and returns the best alignment
"""
    nwalign(seq1::String, seq2::String; match =1, mismatch=-1, gap=-1)

    This function takes the results of a matrix scored using needleman_wunsch properties,
    then returns the proper alignment of sequences based on the matrix as a tuple of two strings.

    Ex.)
    julia> nwalign("AGGT", "ACGAT")
    ("AGG-T", "ACGAT")

    julia> nwalign("ACAT", "AGAT")
    ("ACAT", "AGAT")

    julia> nwalign("AGGTC", "ACGT")
    ("AGGTC", "ACGT-")

"""
function nwalign(seq1::String, seq2::String; match =1, mismatch=-1, gap=-1)
    if !isDNA(seq1) || !isDNA(seq2) || length(seq1) == 0 || length(seq2) == 0
        throw(ErrorException("Invalid DNA sequence provided"))
    end

    scoremat = nwscorematrix(seq1, seq2; match=match, mismatch=mismatch, gap=gap)
    aligned1 = ""
    aligned2 = ""

    i = size(scoremat,1)
    j = size(scoremat,2)

    while i > 1 && j > 1
            #check scores from left, right, diag, and if any match current cell, then add that char to the alignments
            if scoremat[i,j] == nwscore(seq1[i-1], nothing; match, mismatch, gap) + scoremat[i-1,j]
                #@info "max is from above"
                aligned1 = aligned1 * seq1[i-1]
                aligned2 = aligned2 * "-"
                #update the traversal indices
                i = i-1
            elseif scoremat[i,j] == nwscore(nothing, seq2[j-1]; match, mismatch, gap) + scoremat[i,j-1]
                #@info "max is from left"
                aligned1 = aligned1 * "-"
                aligned2 = aligned2 * seq2[j-1]
                #update the traversal indices
                j = j-1
            else
                #@info "max is from diag"
                aligned1 = aligned1 * seq1[i-1]
                aligned2 = aligned2 * seq2[j-1]
                #update the traversal indices
                i = i-1
                j = j-1
            end

    end
    aligned1 = reverse(aligned1)
    aligned2 = reverse(aligned2)
    return (aligned1, aligned2)
end


text:
using ModelParameters

function c_dam_level(volume)
    return 156.8 + 0.9463 * volume^0.2922
end

function c_dam_area(volume)
    return 0.0021 * volume^0.762
end

function c_dam_discharge(volume, max_storage)
    discharge = 0.0
    if volume > max_storage
        discharge = 0.001492 * (volume - max_storage)^1.5280
    end

    return max(0.0, discharge)
end

function c_dam_outflow(discharge, irrigation_extraction)
    return discharge + irrigation_extraction
end

Base.@kwdef mutable struct DamNode{A} <: NetworkNode{A}
    @network_node

    max_storage::Float64
    storage_coef::A = Param(0.5, bounds=(0.01, 10.0))

    # Function to use to calculate dam level from storage volume
    calc_dam_level::Function = c_dam_level

    # Function to use to calculate dam surface area from storage volume
    calc_dam_area::Function = c_dam_area

    # Function to use to calculate dam discharge from storage volume
    calc_dam_discharge::Function = c_dam_discharge

    # Function to calculate outflow from dam
    calc_dam_outflow::Function = c_dam_outflow

    storage::Array{Float64} = [0.0]

    effective_rainfall::Array{Float64} = []
    et::Array{Float64} = []
    inflow::Array{Float64} = []
    dam_area::Array{Float64} = []

    level::Array{Float64} = []
    discharge::Array{Float64} = []
    outflow::Array{Float64} = []

end

function DamNode(
    name::String,
    area::Float64,
    max_storage::Float64,
    storage_coef::Float64,
    initial_storage::Float64,
    calc_dam_level::Function,
    calc_dam_area::Function,
    calc_dam_discharge::Function,
    calc_dam_outflow::Function)
    return DamNode(name, area, max_storage, storage_coef,
                   calc_dam_level, calc_dam_area, calc_dam_discharge, calc_dam_outflow,
                   [initial_storage], [], [], [], [], [], [], [])
end

function level(node::DamNode)
    return last(node.level)
end

function area(node::DamNode)
    return last(node.dam_area)
end

function discharge(node::DamNode)
    return last(node.discharge)
end

function storage(node::DamNode)
    return last(node.storage)
end

function update_state(node::DamNode, storage, rainfall, et, area, discharge, outflow)
    push!(node.storage, storage)
    push!(node.effective_rainfall, rainfall)
    push!(node.et, et)
    push!(node.level, node.calc_dam_level(storage))
    push!(node.dam_area, area)
    push!(node.discharge, discharge)
    push!(node.outflow, outflow)

    return nothing
end

"""
Update dam volume for timestep

Parameters
----------
volume : float, current water volume in ML
node_inflow : float, inflow from previous node in ML
gamma : float, groundwater exchange (positive is gain from gw flow, negative is loss to infiltration)
rain : float, rainfall input
evap : float, evaporation loss
infiltration : float, infiltration loss
area : float, dam surface area in square kilometers
extractions : float, water extraction from dam in ML
discharge : float, discharge from dam in ML
max_store : float, maximum dam storage in ML

Returns
-------
float, volume of water stored in dam
"""
function update_volume(volume, node_inflow, gamma, rain, evap, area, extractions, discharge, max_store)::Float64

    vol = volume + (node_inflow + gamma) + (rain - evap) * area - extractions - discharge

    return max(0.0, min(max_store, vol))
end

"""
Calculate outflow for the dam node for a single time step.

Parameters
----------
node : DamNode
rain : Float64, rainfall in mm
et : Float64, evapotranspiration data in mm
irrig_ext : Float64, irrigation extractions
extractions : Float64, extraction data in ML

:returns: numeric, outflow from Dam
"""
function run_node!(node::DamNode,
                   rain::Float64,
                   et::Float64,
                   inflow::Float64,
                   extractions::Float64,
                   groundwater_gain_loss::Float64=0.0)

    volume = storage(node)
    dam_area = node.calc_dam_area(volume)
    discharge = node.calc_dam_discharge(volume, node.max_storage)

    updated_store = update_volume(volume, inflow, groundwater_gain_loss, rain, et,
                                  dam_area, extractions, discharge, node.max_storage)
    outflow = node.calc_dam_outflow(discharge, extractions)

    update_state(node, updated_store, rain, et, dam_area, discharge, outflow)

    return outflow, level(node)
end

function reset!(node::DamNode)::Nothing
    node.storage = [node.storage[1]]
    node.effective_rainfall = []
    node.et = []
    node.level = []
    node.dam_area = []

    node.discharge = []
    node.outflow = []

    return nothing
end

"""
"""
function update_params!(node::DamNode, params::Array{Float64, 1})::Nothing
    node.storage_coef = params[1]
    return nothing
end

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


text:
module Provenance

using Dates
using UUIDs
using Pkg

using DocStringExtensions

@enum ActivityStatus unfinished finished

struct SystemProvenance
end

struct SystemState
end

mutable struct Activity
    uuid::UUID
    name::String
    parent_activity::Union{UUID, Nothing}
    child_activities::Vector{UUID}
    start::Union{DateTime, Nothing}
    stop::Union{DateTime, Nothing}
    system::SystemProvenance
    input::Vector{String}
    output::Vector{String}
    samples::Vector{SystemState}
    status::ActivityStatus

    Activity(name; parent_activity=nothing) = new(
        uuid4(),
        name,
        parent_activity,
        [],
        now(),
        nothing,
        SystemProvenance(),
        [],
        [],
        [],
        unfinished
    )
end

const activities = Activity[]
push!(activities, Activity("Main activity"))

const backlog = Activity[]

function startactivity(name)
    parent_activity = activities[end]
    activity = Activity(name, parent_activity=parent_activity.uuid)
    push!(activities, activity)
    push!(parent_activity.child_activities, activity.uuid)
    nothing
end

function endactivity()
    activity = pop!(activities)
    finish(activity)
    push!(backlog, activity)
    nothing
end

function finish(activity::Activity)
    activity.status = finished
    nothing
end

end # module


text:
using Vizinanigans

n = 4
states = [randn(10,10) for i in 1:n]

visualize(states)


text:
using Revise
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using NNFEM
Random.seed!(233)

# function viscoelasticity_matrices(mu,eta,lambda,dt)
#     viscoelasticity_matrices_ = load_op_and_grad("./build/libViscoelasticityMatrices","viscoelasticity_matrices", multiple=true)
#     mu,eta,lambda,dt = convert_to_tensor(Any[mu,eta,lambda,dt], [Float64,Float64,Float64,Float64])
#     viscoelasticity_matrices_(mu,eta,lambda,dt)
# end

# TODO: specify your input parameters
domain = example_domain(10,10,0.1)
init_nnfem(domain)
N = getNGauss(domain)
mu = rand(N)
eta = rand(N)
lambda = rand(N)
dt = 0.01
S = zeros(N,3,3)
H = zeros(N,3,3)
for i = 1:N
    p = mu[i] * dt / eta[i]
    S[i,:,:] = inv(
        [
            1+2/3*p -1/3*p 0.0
            -1/3*p 1+2/3*p 0.0
            0.0 0.0 1+p
        ]
    )
    H[i,:,:] = S[i,:,:] * [
        2*mu[i] + lambda[i] lambda[i] 0.0
        lambda[i] 2mu[i] + lambda[i] 0.0
        0.0 0.0 mu[i]
    ]
end
# s, h = viscoelasticity_matrices(mu,eta,lambda,dt)
s, h = compute_maxwell_viscoelasticity_matrices(mu,lambda,eta,dt)

sess = Session(); init(sess)
@show S0, H0 = run(sess, [s,h])
@show maximum(abs.(S0-S))
@show maximum(abs.(H0-H))

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    # return sum(viscoelasticity_matrices(mu,eta,lambda,dt)[2]^2)
    return sum(compute_maxwell_viscoelasticity_matrices(mu,eta,m,dt)[2]^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(N))
v_ = rand(N)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using HttpServer
using WebSockets
import JSON

include("../ekfslam-sim.jl")

# Global object containing simulation state, updated in real time.
# It is exposed for monitoring, debugging, and flexible visualization.
simdata = ekfsim_setup(10, "../course1.txt")

function start(simdata::SimData, monitor::Function, client::WebSockets.WebSocket)
    simdata.paused = false
    @async sim!(simdata, monitor, [simdata, client])
end

"""
Callback passed to simulation providing access to simulation state during
runtime. This method writes messages to a WebSocket client for browser
visualization.
"""
function monitor(simdata::SimData, client::WebSockets.WebSocket)

    n = simdata.scene.nsteps
    tt, st = simdata.scene.true_track, simdata.scene.slam_track
    pose = simdata.state.x[1:3]

    # The 5 parameters for a rotated covariance ellipse

    # Send latest pose information
    d = Dict("ideal" => Dict("x" => tt[1, n], "y" => tt[2, n], "phi" => tt[3, n]),
             "slam"  => Dict("x" => st[1, n], "y" => st[2, n], "phi" => st[3, n]))
    send_json("tracks", d, client)

    # Send SLAM state
    d = Dict("pose" => pose, "cov" => simdata.state.cov)
    send_json("state", d, client)

    # Send line endpoints for lidar beams from vehicle to observed feature
    if simdata.state_updated && simdata.nz > 0
        lines = laser_lines(simdata.z[:,1:simdata.nz], simdata.state.x[1:3])

        nlines = size(lines, 2)
        for j = 1:nlines
            vx, vy, zx, zy = lines[1:4, j]
            @assert inbounds(vx, vy, simdata.scene) "vx, vy: $vx $vy"
            @assert inbounds(zx, zy, simdata.scene) "zx, zy: $zx $zy"
        end

        send_json("lidar", dict_array(lines, ["x1", "y1", "x2", "y2"]), client)

        # Write feature uncertainty ellipses
        if length(simdata.state.x) > 3
            ellipses = feature_ellipses(simdata.state.x, simdata.state.cov)
            ellipse_keys = ["cx", "cy", "rx", "ry", "phi"]
            send_json("feature-ellipses", dict_array(ellipses, ellipse_keys), client)
        end
    end

    # Write uncertainty ellipse for vehicle position
    let
        l,u = eig(simdata.state.cov[1:2, 1:2])
        vehicle_ellipse = [pose' sqrt(l)' atan2(u[2,1], u[1,1])]'
        ellipse_keys = ["cx", "cy", "vehicle_phi", "rx", "ry", "phi"]
        send_json("vehicle-ellipse", dict_array(vehicle_ellipse, ellipse_keys), client)
    end

    return
end

function feature_ellipses(x, cov)
    # Number of observed features/landmarks
    nf = floor(Int, (length(x) - 3)/2)

    # Allocate rows for x,y and columns for npoints/ellipse times # ellipses
    ellipses = Array{Float64}(5, nf)

    for i = 1:nf
        j = (2*i + 2):(2*i + 3)
        l,u = eig(cov[j,j])
        ellipses[:, i] = [x[j]; sqrt(l); atan2(u[2,1], u[1,1])]
    end
    ellipses
end

"""
Write JSON-formatted message to websocket client.
Recipients expect the schema defined here, so modify with care.
"""
function send_json(name::AbstractString, data::Any, client::WebSockets.WebSocket)
    msg = Dict{AbstractString, Any}()
    msg["type"] = name
    msg["data"] = data
    msg["timestamp"] = time()
    write(client, JSON.json(msg))
end

"""
Create array of Dicts from Matrix `a` using the supplied list of `keys`.
Example:
        julia> a = rand(1:10, 2, 3)
        2x3 Array{Int64,2}:
         10  5  8
          5  3  6

        julia> keys = ["a", "b"]
        2-element Array{String,1}:
         "a"
         "b"

        julia> dict_array(a, keys)
        3-element Array{Dict{String,Int64},1}:
         Dict("b"=>5,"a"=>10)
         Dict("b"=>3,"a"=>5)
         Dict("b"=>6,"a"=>8)
"""
function dict_array{T,S}(a::Matrix{T}, keys::Vector{S})
    assert(length(keys) <= size(a, 1))

    n = size(a, 2)
    da = Vector{Dict{S, T}}(n)

    # Loop over columns of a
    for j = 1:n
        d = Dict{S,T}()

        for (i, key) in enumerate(keys)
            d[key] = a[i,j]
        end
        da[j] = d
    end
    da
end

wsh = WebSocketHandler() do req, client
    println("Handling WebSocket client")
    println("    client.id: ",         client.id)
    println("    client.socket: ",     client.socket)
    println("    client.is_closed: ",  client.is_closed)
    println("    client.sent_close: ", client.sent_close)

    while true

        # Read string from client, decode, and parse to Dict
        msg = JSON.parse(String(copy(read(client))))

        if haskey(msg, "text") && msg["text"] == "ready"
            println("Received update from client: ready")
            send_json("waypoints", dict_array(simdata.scene.waypoints, ["x", "y"]), client)
            send_json("landmarks", dict_array(simdata.scene.landmarks, ["x", "y"]), client)
        end

        if haskey(msg, "text") && msg["text"] == "start"
            start(simdata, monitor, client)
        end

        if haskey(msg, "text") && msg["text"] == "reset"
            simdata.scene.nsteps = 0
            simdata.scene.true_track *= 0
            simdata.scene.slam_track *= 0

            simdata.state.x = initial_pose(simdata.scene)
            simdata.state.cov = zeros(3, 3)

            simdata.vehicle.waypoint_id = 1
            simdata.vehicle.pose = initial_pose(simdata.scene)

            simdata.state_updated = false
            simdata.paused = false
        end

        if haskey(msg, "text") && msg["text"] == "pause"
            simdata.paused = !simdata.paused
            if !simdata.paused
                start(simdata, monitor, client)
            end
        end

    end
end

"""
Serve page(s) and supporting files over HTTP.
Files will not be loaded by browser unless they appear in this list!
Paths are relative to the location of this script (or invocation location?)
"""
httph = HttpHandler() do req::Request, res::Response

    # Handle case where / means index.html
    if req.resource == "/"
        println("serving ", req.resource)
        return Response(readstring("index.html"))
    end

    files = [
    "index.html",
    "js/jquery-2.1.3.min.js",
    "js/d3.min.js",
    "js/wsclient.js"
    ]

    for file in files
        if startswith("/$file", req.resource)
            println("serving /", file)
            return Response(open(readstring, file))
        end
    end

    Response(404)
end

httph.events["error"]  = (client, err) -> println(err)
httph.events["listen"] = (port)        -> println("Listening on $port...")

# Instantiate and start a websockets/http server
server = Server(httph, wsh)
println("Starting WebSocket server.")
run(server, 8000)


text:
module SequentialRoomsDirectedModule

import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase
import ..SequentialRoomsUndirectedModule as SRUM

#####
##### game logic
#####

const NUM_OBJECTS = SRUM.NUM_OBJECTS
const AGENT = SRUM.AGENT
const WALL = SRUM.WALL
const GOAL = SRUM.GOAL
const NUM_ACTIONS = 4

mutable struct SequentialRoomsDirected{R, RNG} <: GW.AbstractGridWorld
    env::SRUM.SequentialRoomsUndirected{R, RNG}
    agent_direction::Int
end

function SequentialRoomsDirected(; R = Float32, num_rooms = 3, range_height_room = 4:6, range_width_room = 7:9, rng = Random.GLOBAL_RNG)
    env = SRUM.SequentialRoomsUndirected(R = R, num_rooms = num_rooms, range_height_room = range_height_room, range_width_room = range_width_room, rng = rng)
    agent_direction = rand(rng, 0:GW.NUM_DIRECTIONS-1)
    env = SequentialRoomsDirected(env, agent_direction)
    GW.reset!(env)
    return env
end

function GW.reset!(env::SequentialRoomsDirected)
    GW.reset!(env.env)
    env.agent_direction = rand(env.env.rng, 0:GW.NUM_DIRECTIONS-1)
    return nothing
end

function GW.act!(env::SequentialRoomsDirected, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    inner_env = env.env
    tile_map = inner_env.tile_map
    agent_position = inner_env.agent_position
    agent_direction = env.agent_direction

    if action in Base.OneTo(2)
        if action == 1
            new_agent_position = GW.move_forward(agent_position, agent_direction)
        else
            new_agent_position = GW.move_backward(agent_position, agent_direction)
        end

        if !tile_map[WALL, new_agent_position]
            tile_map[AGENT, agent_position] = false
            inner_env.agent_position = new_agent_position
            tile_map[AGENT, new_agent_position] = true
        end
    elseif action == 3
        env.agent_direction = GW.turn_left(agent_direction)
    else
        env.agent_direction = GW.turn_right(agent_direction)
    end

    if tile_map[GOAL, inner_env.agent_position]
        inner_env.reward = inner_env.terminal_reward
        inner_env.done = true
    else
        inner_env.reward = zero(inner_env.reward)
        inner_env.done = false
    end

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::SequentialRoomsDirected) = GW.get_height(env.env)
GW.get_width(env::SequentialRoomsDirected) = GW.get_width(env.env)

GW.get_action_names(env::SequentialRoomsDirected) = (:MOVE_FORWARD, :MOVE_BACKWARD, :TURN_LEFT, :TURN_RIGHT)
GW.get_object_names(env::SequentialRoomsDirected) = GW.get_object_names(env.env)

function Base.show(io::IO, ::MIME"text/plain", env::SequentialRoomsDirected)
    tile_map = env.env.tile_map
    small_tile_map = SRUM.get_small_tile_map(tile_map)

    characters = ('☻', '█', '♥', '→', '↑', '←', '↓', '⋅')

    _, height_small_tile_map, width_small_tile_map = size(small_tile_map)

    str = "tile_map:\n"

    for i in 1:height_small_tile_map
        for j in 1:width_small_tile_map
            object = findfirst(@view small_tile_map[:, i, j])
            if isnothing(object)
                char = characters[end]
            elseif object == AGENT
                char = characters[NUM_OBJECTS + 1 + env.agent_direction]
            else
                char = characters[object]
            end
            str = str * char
        end

        if i < height_small_tile_map
            str = str * "\n"
        end
    end

    str = str * "\nreward: $(env.env.reward)"
    str = str * "\ndone: $(env.env.done)"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::SequentialRoomsDirected) = ('w', 's', 'a', 'd')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: SequentialRoomsDirected} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: SequentialRoomsDirected} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: SequentialRoomsDirected} = (SRUM.get_small_tile_map(env.env.env.tile_map), env.env.agent_direction)

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: SequentialRoomsDirected} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: SequentialRoomsDirected} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: SequentialRoomsDirected} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: SequentialRoomsDirected} = env.env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: SequentialRoomsDirected} = env.env.env.done

end # module


text: github_project("QuantEcon/quantecon-notebooks-julia", precompile=true) # should see precompiling


text:
# Here we calculate the effective wavenumber and effective wave amplitudes, without any restriction to the volume fraction of particles and incident wave frequency.

"""
    wavenumbers(ω, medium, specie; kws...)

Returns all the possible effective wavenumbers with positive imaginary part when using a single type of particle called `specie`.
"""
wavenumbers(ω::T, medium::PhysicalMedium, specie::Specie{T}; kws...) where T<:Number = wavenumbers(ω, medium, [specie]; kws...)

function wavenumbers(ω::T, source::AbstractSource, material::Material{Dim,S,Sps}; kws...) where {T,Dim,S<:Shape{Dim},Sps<:Species{T,Dim}}
    return wavenumbers(ω, source.medium, material.species;
        numberofparticles = material.numberofparticles,
        # symmetry = Symmetry(source,material),
        kws...
    )
end

"""
    wavenumbers(ω, medium::PhysicalMedium, species::Species; kws...)

Returns all the possible effective wavenumbers with positive imaginary part. This function requires significantly numerical optimisation and so can be slow.
"""
function wavenumbers(ω::T, medium::PhysicalMedium, species::Species{T};
        num_wavenumbers::Int = 2, tol::T = 1e-5,
        max_Imk::T = T(2) + T(20) * imag(wavenumber_low_volumefraction(ω, medium, species; verbose = false)),
        basis_order = 3 * Int(round(maximum(outer_radius.(species)) * ω / abs(medium.c) )) + 1,
        # max_Rek::T = T(2) + T(20) * abs(real(wavenumber_low_volumefraction(ω, medium, species; verbose = false))),
        kws...) where T<:Number

    # For very low attenuation, need to search close to assymptotic root with a path method.
    k_effs::Vector{Complex{T}} = wavenumbers_path(ω, medium, species;
        num_wavenumbers = num_wavenumbers,
        max_Imk = max_Imk, tol = tol,
        basis_order = basis_order, kws...
    )

    # Take only the num_wavenumbers wavenumbers with the smallest imaginary part.
    k_effs = k_effs[1:num_wavenumbers]

    # NOTE: these search methods would significantly improve if we used the asymptotic result for multiple wavenumbers and monopole scatterers. This would give a reasonable length scale for the mesh, and good estimates on where to start searching.
    if num_wavenumbers > 2
        # box_k = box_keff(ω, medium, species; tol = tol)
        # max_imag = max(3.0 * maximum(imag.(k_effs)), max_Imk)
        # max_imag = max(max_imag, box_k[2][2])
        # max_real = max(2.0 * maximum(real.(k_effs)), max_Rek)

        max_imag = maximum(imag.(k_effs))
        max_real = maximum(real.(k_effs))
        box_k = [[-max_real,max_real], [0.0,max_imag]]

        k_effs2 = wavenumbers_bisection(ω, medium, species;
            # num_wavenumbers=num_wavenumbers,
            tol = tol, box_k = box_k,
            basis_order = basis_order,
            kws...)
        k_effs = [k_effs; k_effs2]
        k_effs = reduce_kvecs(k_effs, tol)
        k_effs = sort(k_effs, by = imag)
    end

    return Complex{Float64}.(k_effs)
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
using Zygote, SparseArrays

Zygote.@adjoint Base.:(-)(a, b) = a-b, Δ -> (Δ, -Δ)

_, back = Zygote._forward(sprand(2, 2, 0.1), sprand(2, 2, 0.2)) do x, y
    x - y
end


text:
## Functions to make running and analyzing persistent homology a little easier.

function createSimplexLists(incidenceMatrix,simplex_named_list)

    # Will make lists of edges, 2-simplices, and 3-simplices
    edge_list = Any[]
    simp2_list = Any[]
    simp3_list = Any[]

    for (index,simplex_name) in enumerate(simplex_named_list)

        indices_of_simplex_nodes = findall(incidenceMatrix[index,:].>0)

        edges_from_simplex = collect(combinations(indices_of_simplex_nodes,2))
        simp2_from_simplex = collect(combinations(indices_of_simplex_nodes,3))
        simp3_from_simplex = collect(combinations(indices_of_simplex_nodes,4))

        for (index2,edge) in enumerate(edges_from_simplex)
            push!(edge_list,edge)
        end
        for (index2,simp2) in enumerate(simp2_from_simplex)
            push!(simp2_list,simp2)
        end
        for (index2,simp3) in enumerate(simp3_from_simplex)
            push!(simp3_list,simp3)
        end
    end
    edge_list = unique(edge_list)
    simp2_list = unique(simp2_list)
    simp3_list = unique(simp3_list)

    return edge_list, simp2_list, simp3_list
end

function createCellList(nNodes,edgeList,keepTriangles)
    n2Simps = length(keepTriangles)
    nEdges = length(edgeList)
    #E = Array{Array{Float64}}((Int(nNodes+nEdges+n2Simps)),1)
    E = Array{Array{Float64,1},1}(undef,Int(nNodes+nEdges+n2Simps))
    nCells = length(E)
    for a0 in collect(1:nCells)

        if a0 <= nNodes
            E[a0] = [a0]
        elseif a0 > (nNodes) && a0<=(nNodes+nEdges)
            E[a0] = edgeList[a0-nNodes]
        else
            E[a0] = keepTriangles[a0-(nNodes+nEdges)]

        end

    end
    E
end

function createDMatrix(E,nNodes,nEdges)
    D_01 = Int.(issubset.(E[1:nNodes],E[(nNodes+1):(nNodes+nEdges)]'))
    D_12 = Int.(issubset.(E[(nNodes+1):(nNodes+nEdges)],E[(nNodes+nEdges+1):length(E)]'))
    D = zeros(length(E),length(E))
    D[1:nNodes,(nNodes+1):(nNodes+nEdges)] = D_01
    D[(nNodes+1):(nNodes+nEdges),(nNodes+nEdges+1):length(E)] = D_12
    D = sparse(D)
end

function determineCellWeight(cell)
    cellWeight = maximum(node_weights[Int.(cell)])
end

function bettiCurveFromBarcode(barcode,nNodes,node_weights)

    nNodes = Int(nNodes)
    bettiCurve = zeros(nNodes+1,1)
    ordered_weights = sort(unique(node_weights))

    for bar in collect(1:size(barcode,1))

        birth_index = Int(findall(ordered_weights.==barcode[bar,1])[1])

        if barcode[bar,2]>nNodes

            bettiCurve[birth_index:Int(nNodes+1),1] = bettiCurve[birth_index:Int(nNodes+1),1] .+1
        else

            death_index = Int(findall(ordered_weights.==barcode[bar,2])[1])
            bettiCurve[birth_index:death_index,1] = bettiCurve[birth_index:death_index,1].+1

        end
    end

    return bettiCurve
end

function plotBarcode(barcode,nNodes)

    nNodes = Int(nNodes)
    pbar = plot(1:6,zeros(6,1),c=:black)

    barcode_sorted = barcode[sortperm(barcode[:,1]),:]
    nbars = size(barcode)[1]

    for bar_index in collect(1:nbars)
        birth = barcode[bar_index,1]
        death = barcode[bar_index,2]
        if death>nNodes
            death = nNodes+1
        end

        plot!([birth, death],[bar_index, bar_index], legend = false,
                            xlim = (0,nNodes))
    end

    return pbar
end


text:
using DexApiClients
using DexApiClients.api
using gRPCClient
using Test

function test_version(dex, api_ver)
    inp = VersionReq()
    out, status = GetVersion(dex, inp)
    gRPCCheck(status)

    @test isa(out, VersionResp)
    @test out.api == api_ver
    nothing
end

function test_list_passwords(dex, min_npwd, find_user)
    inp = ListPasswordReq()
    out, status = ListPasswords(dex, inp)
    gRPCCheck(status)

    @test isa(out, ListPasswordResp)
    @test length(out.passwords) >= min_npwd
    if find_user !== nothing
        filtered_passwords = filter(find_user, out.passwords)
        @test !isempty(filtered_passwords)
    end
    nothing
end

function test_update_password(dex, email, new_hash, new_username, expect_not_found)
    inp = UpdatePasswordReq(;email=email, new_hash=new_hash, new_username=new_username)
    out, status = UpdatePassword(dex, inp)
    gRPCCheck(status)

    @test isa(out, UpdatePasswordResp)
    @test out.not_found == expect_not_found
    nothing
end

function test_create_password(dex, email, hash, username, user_id)
    pwd = Password(;email=email, hash=hash, username=username, user_id=user_id)
    inp = CreatePasswordReq(;password=pwd)
    out, status = CreatePassword(dex, inp)
    gRPCCheck(status)

    @test isa(out, CreatePasswordResp)
    @test !out.already_exists
    nothing
end

function test_delete_password(dex, email, expect_not_found)
    inp = DeletePasswordReq(;email=email)
    out, status = DeletePassword(dex, inp)
    gRPCCheck(status)

    @test isa(out, DeletePasswordResp)
    @test out.not_found == expect_not_found
    nothing
end

function test_all()
    dex = DexBlockingClient("http://127.0.0.1:5557"; verbose=false)

    test_version(dex, 2)
    test_list_passwords(dex, 1, nothing)

    test_update_password(dex, "abc@def.com", Vector{UInt8}("\$2a\$10\$2b2cU8CPhOTaGrs1HRQuAueS7JTT5ZHsHSzYiFPm1leZck7Mc8T4W"), "new abc", true)
    test_delete_password(dex, "abc@def.com", true)

    test_create_password(dex, "abc@def.com", Vector{UInt8}("\$2a\$10\$2b2cU8CPhOTaGrs1HRQuAueS7JTT5ZHsHSzYiFPm1leZck7Mc8T4W"), "abc", "100")
    test_list_passwords(dex, 2, x->((x.email == "abc@def.com") && (x.username == "abc") && (x.user_id == "100")))
    test_update_password(dex, "abc@def.com", Vector{UInt8}("\$2a\$10\$2b2cU8CPhOTaGrs1HRQuAueS7JTT5ZHsHSzYiFPm1leZck7Mc8T4W"), "new abc", false)
    test_list_passwords(dex, 2, x->((x.email == "abc@def.com") && (x.username == "new abc") && (x.user_id == "100")))
    test_delete_password(dex, "abc@def.com", false)
end


text:
export iauA2tf
"""
Decompose radians into hours, minutes, seconds, fraction.

 This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

 Status:  vector/matrix support function.

 Given:
   ndp     int     resolution (Note 1)
   angle   double  angle in radians

 Returned:
   sign    char    '+' or '-'
   ihmsf   int[4]  hours, minutes, seconds, fraction

 Called:
   iauD2tf      decompose days to hms

 Notes:

    1. The argument ndp is interpreted as follows:
| ndp | resolution |
| :-: | :--------: |
| : | ...00,000,000 |
| -7 | 10,000,000 |
| -6 | 1,000,000 |
| -5 | 100,000 |
| -4 | 10,000 |
| -3 | 1,000 |
| -2 | 100 |
| -1 | 10 |
| 0 | 1 |
| 1 | 0.1 |
| 2 | 0.01 |
| 3 | 0.001 |
| : | 0.000... |

    2. The largest positive useful value for ndp is determined by the
    size of angle, the format of doubles on the target platform, and
    the risk of overflowing ihmsf[4].  On a typical platform, for
    angle up to 2pi, the available floating-point precision might
    correspond to ndp=12.  However, the practical limit is typically
    ndp=9, set by the capacity of a 32-bit int, or ndp=4 if int is
    only 16 bits.

    3. The absolute value of angle may exceed 2pi.  In cases where it
    does not, it is up to the caller to test for and handle the
    case where angle is very nearly 2pi and rounds up to 24 hours,
    by testing for ihmsf[1]=24 and setting ihmsf[1:4] to zero.

This revision:  2013 July 31

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
function iauA2tf(ndp::Int, angle::Real)
    sign  = Array{UInt8}(undef, 1)
    ihmsf = Array{Int32}(undef, 4)

    ccall((:iauA2tf, libsofa_c), Cvoid,
        (Cint, Cdouble, Ptr{UInt8}, Ptr{Cint}),
        convert(Int32, ndp), convert(Float64, angle), sign, ihmsf)

    return Char(sign[1]), ihmsf
end


text:
# this recursively flattens all dict's in a JSON structure
function flatten(d::Dict{String, Any}; prefix::String = "")
    d_out = Dict{String, Any}()
    for k in keys(d)
        if isa(d[k], Dict)
            merge!(d_out, flatten(d[k], prefix = "$(prefix)$(k)_"))
        else
            push!(d_out, "$(prefix)$k" => d[k])
        end
    end
    d_out
end

# flatten a vector of Dict's to a DataFrame
function flatten_to_df(dicts::Array{Dict{String,Any},1})
    flat = flatten.(dicts)
    dicts_to_df(flat)
end
function flatten_to_df(dicts::Array{Any,1})
    if any(typeof.(dicts) .!= Dict{String, Any})
        error("Vector contains Dict with non-string keys")
    end
    flat = flatten.(dicts)
    dicts_to_df(flat)
end
# flatten a vector of Dict's to a CSV file
function flatten_to_csv(dicts::Array{Dict{String,Any},1}, csvfilename::String)
    flat = flatten.(dicts)
    dicts_to_csv(flat,csvfilename)
end
function flatten_to_csv(dicts::Array{Any,1}, csvfilename::String)
    if any(typeof.(dicts) .!= Dict{String, Any})
        error("Vector contains Dict with non-string keys")
    end
    flat = flatten.(dicts)
    dicts_to_csv(flat,csvfilename)
end

# The Dict's are assumed to NOT contain any further dicts
function dicts_to_df(dicts::Array{Dict{String,Any},1})
    # get variable names
    variables = Vector{String}()
    for i = 1:length(dicts)
        for k in keys(dicts[i])
            if !any(variables .== k)
                push!(variables,k)
            end
        end
    end
    # construct DataFrame
    df = DataFrame(_row = collect(1:length(dicts)))
    for v in variables
        df[:,Symbol(v)] = [(haskey(dicts[i], v) ? dicts[i][v] : missing) for i = 1:length(dicts)]
    end
    return df
end

# write to csv
function dicts_to_csv(dicts::Array{Dict{String,Any},1}, csvfilename::String)
    df = dicts_to_df(dicts)
    CSV.write(csvfilename, df, missingstring="",quotestrings=true)
end

# Additional cleaning code

# cleanvariables = ["content", "highlight", "sgcool_label_text"]
# for v in cleanvariables
#     if (typeof(df[:,Symbol(v)]) == Array{String,1})
#         x = df[!,Symbol(v)]
#         x .= replace.(x, "\r\n" .=> " ");
#         x .= replace.(x, "\n" .=> " ");
#         x .= replace.(x, "\r" .=> " ");
#         x .= replace.(x, "\\r\\n" .=> " ");
#         x .= replace.(x, "\\n" .=> " ");
#         x .= replace.(x, "\\r" .=> " ");
#     end
#     if (typeof(df[:,Symbol(v)]) == Array{Union{Missing, String},1})
#         x = skipmissing(df[!,Symbol(v)])
#         df[.!(ismissing.(df[:,Symbol(v)])),Symbol(v)] = replace.(x, "\r\n" .=> " ");
#         df[.!(ismissing.(df[:,Symbol(v)])),Symbol(v)] = replace.(x, "\n" .=> " ");
#         df[.!(ismissing.(df[:,Symbol(v)])),Symbol(v)] = replace.(x, "\r" .=> " ");
#         df[.!(ismissing.(df[:,Symbol(v)])),Symbol(v)] = replace.(x, "\\r\\n" .=> " ");
#         df[.!(ismissing.(df[:,Symbol(v)])),Symbol(v)] = replace.(x, "\\n" .=> " ");
#         df[.!(ismissing.(df[:,Symbol(v)])),Symbol(v)] = replace.(x, "\\r" .=> " ");
#     end
# end

# CSV.write(outfilename, df, missingstring="",quotestrings=true)


text:
# using Pkg
# Pkg.activate("c:/git/CuCountmap")
using SortingAlgorithms
using CuCountMap

v = rand(UInt64, 100_000_000);

@benchmark sorted_v = gpuradixsort($v)
# BenchmarkTools.Trial:
#   memory estimate:  1.49 GiB
#   allocs estimate:  2550
#   --------------
#   minimum time:     3.351 s (0.10% GC)
#   median time:      3.404 s (1.82% GC)
#   mean time:        3.404 s (1.82% GC)
#   maximum time:     3.457 s (3.48% GC)
#   --------------
#   samples:          2
#   evals/sample:     1

@benchmark sorted_v_cpu = sort($v, alg=RadixSort)
# BenchmarkTools.Trial:
#   memory estimate:  1.49 GiB
#   allocs estimate:  18
#   --------------
#   minimum time:     3.648 s (0.08% GC)
#   median time:      3.708 s (1.75% GC)
#   mean time:        3.708 s (1.75% GC)
#   maximum time:     3.767 s (3.36% GC)
#   --------------
#   samples:          2
#   evals/sample:     1


text:
module ServoWithRampAndRevolute

using ModiaLang
using Unitful
import Modia3D

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

Controller = Model(
    # Interface
    refLoadAngle     = input,
    actualMotorAngle = input,
    actualMotorSpeed = input,
    refTorque        = output,
    # Components
    gainOuter     = Gain,
    angleFeedback = Feedback,
    gainInner     = Gain,
    speedFeedback = Feedback,
    PI = PI,

    connect = :[
        (refLoadAngle, gainOuter.u)
        (gainOuter.y, angleFeedback.u1)
        (actualMotorAngle, angleFeedback.u2)
        (angleFeedback.y, gainInner.u)
        (gainInner.y, speedFeedback.u1)
        (actualMotorSpeed, speedFeedback.u2)
        (speedFeedback.y, PI.u)
        (PI.y, refTorque)]
)

Drive = Model(
    # Interface
    refTorque   = input,
    actualSpeed = output,
    actualAngle = output,
    flange      = Flange,
    # Components
    torque       = UnitlessTorque,
    speedSensor  = UnitlessSpeedSensor,
    angleSensor  = UnitlessAngleSensor,
    motorInertia = Inertia,
    idealGear    = IdealGear,

    connect = :[
        (refTorque, torque.tau)
        (torque.flange, speedSensor.flange, angleSensor.flange, motorInertia.flange_a)
        (motorInertia.flange_b, idealGear.flange_a)
        (idealGear.flange_b, flange)
        (speedSensor.w, actualSpeed)
        (angleSensor.phi, actualAngle)
        ]
)

Servo = Model(
    refLoadAngle = input,
    flange       = Flange,
    ppi          = Controller,
    drive        = Drive,

    connect = :[
        # inputs of ppi
        (refLoadAngle, ppi.refLoadAngle)
        (ppi.actualMotorSpeed, drive.actualSpeed)
        (ppi.actualMotorAngle, drive.actualAngle)
        # output of ppi --> input of motor
        (ppi.refTorque, drive.refTorque)
        # Flange of drive
        (drive.flange, flange)
        ]
)

arm_joint_1_obj = joinpath(Modia3D.path, "objects/robot_KUKA_YouBot/arm_joint_1.obj")

m = 2
axis = 3
vmat1 = VisualMaterial(color="LightBlue", transparency=0.5)
m1=1.390
translation1 =[0.033,0,0]
rotation1 =Modia3D.rot1(180u"°")

k1=50.0
k2=0.1
T2=1.0
gearRatio=156.0
motorInertia = (0.0000135 + 0.000000409)*u"kg*m^2"
J1 = 0.32

servoParameters = Map(
                      ppi = Map(
                          gainOuter = Map(k=gearRatio),
                          gainInner = Map(k=k1),
                          PI        = Map(k=k2, T=T2)
                      ),
                      drive = Map(
                          motorInertia = Map(J=motorInertia),
                          idealGear = Map(ratio=gearRatio)
                      )
                  )

TestServo = Model(
    world = Object3D(feature=Scene()),
    body  = Object3D(feature = Solid(
            shape = FileMesh(filename = arm_joint_1_obj), massProperties = MassPropertiesFromShapeAndMass(mass=m1), visualMaterial=vmat1)),
    obj2  = Object3D(parent = :body, translation = translation1, rotation = rotation1),

    rev = RevoluteWithFlange(obj1 = :world, obj2 = :obj2, axis=axis, phi = Var(init = 0.0), w=Var(init=0.0)),

    ramp  = Ramp | Map(duration=1.18u"s", height=2.95),
    servo = Servo | servoParameters,
    connect = :[
        (ramp.y, servo.refLoadAngle)
        (servo.flange, rev.flange)
        ]
)

servo = @instantiateModel(buildModia3D(TestServo), unitless=true, logCode=true, log=true)

stopTime = 4.0
tolerance = 1e-6
requiredFinalStates = [2.949999509178657, 4.908194611512363e-7, 0.02159751818019237]
simulate!(servo, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plotVariables = [("ramp.y", "rev.phi"); "rev.w"; "servo.ppi.PI.x"]
plot(servo, plotVariables, figure=1)

end


text:
"""
Module to provide file inut and output functionality for `Rimu`.
Provides convenience functions:

* [`RimuIO.save_df(filename, df::DataFrame)`](@ref) Save dataframe in Arrow format.
* [`RimuIO.load_df(filename)`](@ref) Load Arrow file into dataframe.
* [`RimuIO.save_dvec(filename, dv)`](@ref) Save dict vector in BSON format.
* [`RimuIO.load_dvec(filename)`](@ref) Load BSON file into dict vector.
"""
module RimuIO

using ..DictVectors
using ..Interfaces

using Arrow, DataFrames, BSON

export save_df, load_df, save_dvec, load_dvec

const COMPLEX = Symbol("JuliaLang.Complex")
ArrowTypes.arrowname(::Type{<:Complex}) = COMPLEX
ArrowTypes.ArrowType(::Type{Complex{T}}) where {T} = Tuple{T,T}
ArrowTypes.JuliaType(::Val{COMPLEX}, ::Type{Tuple{T,T}}) where {T} = Complex{T}
ArrowTypes.toarrow(a::Complex) = (a.re, a.im)
ArrowTypes.fromarrow(::Type{Complex{T}}, t::Tuple{T,T}) where {T} = Complex{T}(t...)

"""
    RimuIO.save_df(filename, df::DataFrame)
Save dataframe in Arrow format.
"""
save_df(filename, df::DataFrame) = Arrow.write(filename, df)

"""
    RimuIO.load_df(filename) -> DataFrame
Load Arrow file into dataframe.
"""
load_df(filename) = DataFrame(Arrow.Table(filename))

"""
    RimuIO.save_dvec(filename, dvec)
Save `dvec` in [BSON](https://github.com/JuliaIO/BSON.jl) format.

## Notes

* Only the [`localpart`](@ref) is saved. You may need to re-wrap the result in
  [`MPIData`](@ref) if using MPI.

* When using this function with MPI, make sure to save the vectors from different ranks to
  different files, e.g. by saving as `RimuIO.save_dvec("filename-\$(mpi_rank()).bson", dvec)`.

"""
save_dvec(filename, dvec) = bson(filename, Dict(:dvec => localpart(dvec)))

"""
    RimuIO.load_dvec(filename) -> AbstractDVec
Load `AbstractDVec` stored in [BSON](https://github.com/JuliaIO/BSON.jl).
"""
load_dvec(filename) = BSON.load(filename)[:dvec]

end # module RimuIO


text:
##for debugging
# cd(dirname(@__FILE__))
run(`make clean -C apriltag-0.10.0/`)

# Make c library
run(`make -C apriltag-0.10.0/`)

# set linker path to library
libpath = joinpath(dirname(@__FILE__),"apriltag-0.10.0")

#create a file to push to the DL_LOAD_PATH when module loads,
f = open("loadpath.jl", "w")

write(f,"# This is an automatically generated file, do not edit.\n")
write(f,"using Libdl\n")
write(f,"push!(Libdl.DL_LOAD_PATH, \"$(libpath)\")")

close(f)


text:
# Pairwise Function References

base_functions = (
    SquaredEuclidean,
    SineSquared,
    ChiSquared,
    ScalarProduct
)

base_functions_initiate = Dict(
    SquaredEuclidean => 0,
    SineSquared      => 0,
    ChiSquared       => 0,
    ScalarProduct    => 0
)

base_functions_aggregate = Dict(
    SquaredEuclidean => (s,x,y) -> s + (x-y)^2,
    SineSquared      => (s,x,y) -> s + sin((x-y))^2,
    ChiSquared       => (s,x,y) -> x == y == 0 ? s : s + ((x-y)^2/(x+y)),
    ScalarProduct    => (s,x,y) -> s + x*y
)

base_functions_return = Dict{DataType,Any}()

# Kernel Function References

kernel_functions = (
    ExponentialKernel,
    SquaredExponentialKernel,
    GammaExponentialKernel,
    RationalQuadraticKernel,
    GammaRationalKernel,
    MaternKernel,
    LinearKernel,
    PolynomialKernel,
    ExponentiatedKernel,
    PeriodicKernel,
    PowerKernel,
    LogKernel,
    SigmoidKernel
)

kernel_functions_arguments = Dict(
    ExponentialKernel        => ((1.0,),        (2.0,)),
    SquaredExponentialKernel => ((1.0,),        (2.0,)),
    GammaExponentialKernel   => ((1.0,1.0),     (2.0,0.5)),
    RationalQuadraticKernel  => ((1.0,1.0),     (2.0,2.0)),
    GammaRationalKernel      => ((1.0,1.0,1.0), (2.0,2.0,0.5)),
    MaternKernel             => ((1.0,1.0),     (2.0,2.0)),
    LinearKernel             => ((1.0,1.0),     (2.0,2.0)),
    PolynomialKernel         => ((1.0,1.0,3),   (2.0,2.0,2)),
    ExponentiatedKernel      => ((1.0,),        (2.0,)),
    PeriodicKernel           => ((1.0,),        (2.0,)),
    PowerKernel              => ((1.0,),        (0.5,)),
    LogKernel                => ((1.0,1.0),     (2.0,0.5)),
    SigmoidKernel            => ((1.0,1.0),     (2.0,2.0))
)

kernel_functions_kappa = Dict(
    ExponentialKernel        => (z,α)     -> exp(-α*sqrt(z)),
    SquaredExponentialKernel => (z,α)     -> exp(-α*z),
    GammaExponentialKernel   => (z,α,γ)   -> exp(-α*z^γ),
    RationalQuadraticKernel  => (z,α,β)   -> (1 + α*z)^(-β),
    GammaRationalKernel      => (z,α,β,γ) -> (1 + α*z^γ)^(-β),
    MaternKernel             => (z,ν,θ)   -> begin
                                                v1 = sqrt(2*ν) * z / θ
                                                v1 = v1 < eps(typeof(z)) ? eps(typeof(z)) : v1
                                                2*(v1/2)^ν * besselk(ν,v1)/gamma(ν)
                                             end,
    LinearKernel             => (z,a,c)   -> (a*z+c),
    PolynomialKernel         => (z,a,c,d) -> (a*z+c)^d,
    ExponentiatedKernel      => (z,a)     -> exp(a*z),
    PeriodicKernel           => (z,α)     -> exp(-α*z),
    PowerKernel              => (z,γ)     -> z^γ,
    LogKernel                => (z,α,γ)   -> log(α*z^γ+1),
    SigmoidKernel            => (z,a,c)   -> tanh(a*z+c)
)

kernel_functions_base = Dict(
    LinearKernel             => ScalarProduct,
    PolynomialKernel         => ScalarProduct,
    ExponentiatedKernel      => ScalarProduct,
    PeriodicKernel           => SineSquared,
    SigmoidKernel            => ScalarProduct
)

base_functions_properties = Dict(
                       #|stnry |isotrop
    SineSquared      => (true,  false),
    ScalarProduct    => (false, false),
    SquaredEuclidean => (true,  true),
    ChiSquared       => (false, false)
 )


text:
"""
#### function ```get_area(x1::Float64, x2::Float64, a::Float64, b::Float64)```

    Description of ```get_area```
    ----------------------------------------
Returns the area under trapezoid, where "roof" is defined by a linear function y = ax +b.

"""
function get_area(x1::Float64, x2::Float64, a::Float64, b::Float64)
    # area under the trapezoid is A = Δx * ⟨y⟩
    # where Δx is the width of the base of trapezoid
    # and ⟨y⟩ is mean height

    Δx = abs(x2-x1)
    mean_y = abs( get_line_point(x2,a,b) + get_line_point(x1,a,b)) / 2
    return Δx * mean_y
end

"""
#### function ```get_row_volume(df_row::DataFrame, thickness = 0.001)```

    Description of ```get_row_volume```
    ----------------------------------------
Returns the volume of trapezoid, where "roof" is defined by a linear function y = ax +b and thickness is the step between rows - ΔE.
```df_row``` is the df subset where df.E1 .== e1, e1 is the constant energy for which the linear approximation was made.

"""
function get_row_volume(df_row::DataFrame, thickness = 0.001)
    totalArea = sum(get_area.(df_row.minE, df_row.maxE, df_row.a, df_row.b))
    volume = totalArea*thickness
    return volume
end

"""
#### function ```get_total_volume(df_rebinned::DataFrame, thickness = 0.001)```

    Description of ```get_total_volume```
    ----------------------------------------
Returns the total volume as a sum of all rows. Thickness is the step between energies, ΔE.

"""
function get_total_volume(df_rebinned::DataFrame, thickness = 0.001)
    volume = 0.0
    for e in unique(df_rebinned.E1)
        volume += get_row_volume(df_rebinned[df_rebinned.E1 .== e, :], thickness)
    end
    return volume
end

"""
#### function ```get_segment_volume(df::DataFrame, X1::Float64, X2::Float64, thickness = 0.001)```

    Description of ```get_segment_volume```
    ----------------------------------------
Returns the volume of a segment enclosing x1, x2.

"""
function get_segment_volume(df::DataFrame, X1::Float64, X2::Float64, thickness = 0.001)
    X1 < minimum(df.minE) && (X1 = minimum(df.minE))  #if x1 is less than the lowest minE value in df, set X1 to the lowest value
    X2 > maximum(df.maxE) && (X2 = maximum(df.maxE))  #if x2 is more than the highest minE value in df, set X2 to the highest value

    df_segment = subset(df, :minE => e -> e .<= X1, :maxE => e -> e .> X1, view = :true) # return df subset where minE < X1 < maxE

    if nrow(df_segment) == 0
        return 0.0
    end

    A = 0.0

    if(X1 > df_segment.minE[1]) && (X2 < df_segment.maxE[1] )                     # if the bounds are within one line region
        A += get_area(X1, X2, df_segment.a[1], df_segment.b[1])
        return A*thickness
    end

    while df_segment.maxE[1] - X2 < 0
        A += get_area(X1, df_segment.maxE[1], df_segment.a[1], df_segment.b[1])
        X1 = df_segment.maxE[1]
        df_segment = subset(df, :minE => e -> e .<= X1, :maxE => e -> e .> X1, view = :true)

    end
    A += get_area(df_segment.minE[1], X2, df_segment.a[1], df_segment.b[1])
    return A*thickness
end

"""
#### function ```get_volume_matrix(df::DataFrame, eMin, eMax, dE, thickness = 0.001)```

    Description of ```get_volume_matrix```
    ----------------------------------------
Returns a square matrix of volumes. Each cell defines volume of a segment defined by dE1*dE2 = dE^2.

"""
function get_volume_matrix(df::DataFrame, eMin, eMax, dE, thickness = 0.001) #returns square matrix with volumes
    dims = length(minimum(df.E1):dE:maximum(df.E1))  # dimension of the square matrix

    volumesTruncated = zeros(dims,dims) # final volume to output with dimensions n x n
    temp_volumes = zeros(dims, length(unique(df.E1)))  # 30x2997

    for (idxE1, uniqueE1) in enumerate(unique(df.E1))  #1:2997
        df_e1 = subset(df, :E1 => e -> e .== uniqueE1)

        for (idxEnergy, energy) in enumerate(eMin:dE:eMax) # 1:30
            temp_volumes[idxEnergy, idxE1] = get_segment_volume(df_e1, energy, energy+dE, thickness)
        end
    end

    n = floor(Int,length(unique(df.E1))/ dims )  # how many cols to truncate
    for r in 1:length(temp_volumes[:,1])              # truncating n cols into 1
        volumesTruncated[r,:] = [sum(temp_volumes[r,i:i+n]) for i in 1:n:length(temp_volumes[r,:])-n]
    end
    return volumesTruncated
end


text:
using Oceananigans
using Oceananigans.Units
using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation: Transformation
using Distributions
using DataDeps
using GLMakie
using Printf

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

θ_constant_Ri = (;
                  CᵂwΔ  = 8.38,
                  Cᵂu★  = 9.06,
                  Cᴰ    = 1.37,
                  Cᴸᵇ   = 1.13,
                  Cᴷu⁻  = 0.0887,
                  Cᴷc⁻  = 0.164,
                  Cᴷe⁻  = 2.13,
                 )

function perturbation_prior(θ★, ϵ=0.7)
    L = θ★ / 4
    U = 2θ★
    return ScaledLogitNormal(bounds=(L, U))
end

θ★ = θ_constant_Ri
priors = NamedTuple(name => perturbation_prior(θ★[name]) for name in keys(θ★))
free_parameters = FreeParameters(priors)

case_path(case) = @datadep_str("two_day_suite_4m/$(case)_instantaneous_statistics.jld2")

case = "weak_wind_strong_cooling"

times = range(2hours, step=10minutes, stop=48hours)
Nt = length(times)
transformation = Transformation(time=TimeIndices([2, round(Int, Nt/2), Nt]))
field_names = (:u, :v, :T, :e)
regrid_size = nothing
observations = SyntheticObservations(case_path(case); field_names, regrid_size, times, transformation)

α = observations.metadata.parameters.thermal_expansion_coefficient
equation_of_state = LinearEquationOfState(; α)

mixing_length = MixingLength(Cᴬu   = 0.0,
                             Cᴬc   = 0.0,
                             Cᴬe   = 0.0,
                             CᴷRiᶜ = 2.0,
                             Cᴷuʳ  = 0.0,
                             Cᴷcʳ  = 0.0,
                             Cᴷeʳ  = 0.0)

catke = CATKEVerticalDiffusivity(; mixing_length)

Nensemble = 40

simulation = ensemble_column_model_simulation(observations;
                                              Nensemble,
                                              architecture = GPU(),
                                              buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=true),
                                              tracers = (:T, :e),
                                              closure = catke)

simulation.Δt = 2.0

Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
Qᵀ = simulation.model.tracers.T.boundary_conditions.top.condition
dTdz = simulation.model.tracers.T.boundary_conditions.bottom.condition

Qᵘ .= observations.metadata.parameters.momentum_flux
Qᵀ .= observations.metadata.parameters.temperature_flux
dTdz .= observations.metadata.parameters.dθdz_deep

calibration = InverseProblem(observations, simulation, free_parameters)

output_paths = []

model = simulation.model

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, model.diffusivity_fields),
                     prefix = "catke_calibration_0",
                     schedule = SpecifiedTimes(times...),
                     force = true)

push!(output_paths, simulation.output_writers[:fields].filepath)

eki = EnsembleKalmanInversion(calibration; noise_covariance=1e-2)

for i = 1:5
    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, merge(model.velocities, model.tracers, model.diffusivity_fields),
                         prefix = string("catke_calibration_", i),
                         schedule = SpecifiedTimes(times...),
                         force = true)

    push!(output_paths, simulation.output_writers[:fields].filepath)

    iterate!(eki)
end

u = []
v = []
T = []
e = []

for i = 0:5
    @show i
    path = string("catke_calibration_", i, ".jld2")
    push!(u, FieldTimeSeries(path, "u", boundary_conditions=nothing))
    push!(v, FieldTimeSeries(path, "v", boundary_conditions=nothing))
    push!(T, FieldTimeSeries(path, "T", boundary_conditions=nothing))
    push!(e, FieldTimeSeries(path, "e", boundary_conditions=nothing))
end

z = znodes(first(u))

fig = Figure(resolution=(1200, 600))

ax_T = Axis(fig[1, 1], xlabel = "Temperature (ᵒC)", ylabel = "z (m)")
ax_u = Axis(fig[1, 2], xlabel = "Velocity (m s⁻¹)", ylabel = "z (m)")
ax_e = Axis(fig[1, 3], xlabel = "Turbulent Kinetic Energy (m² s⁻²)", ylabel = "z (m)")

Nt = length(observations.times)
time_slider = Slider(fig[2, :], range=1:Nt, horizontal=true, startvalue=1)
n = time_slider.value

iter_slider = Slider(fig[3, :], range=1:length(u), horizontal=true, startvalue=1)
iter = iter_slider.value

profile_max(u, i, n) = Point2f.(maximum(interior(u[n])[:, i, :], dims=1)[:], z)
profile_min(u, i, n) = Point2f.(minimum(interior(u[n])[:, i, :], dims=1)[:], z)

#=
u_max = @lift profile_max(u[$iter], 1, $n)
v_max = @lift profile_max(v[$iter], 1, $n)
T_max = @lift profile_max(T[$iter], 1, $n)
e_max = @lift profile_max(e[$iter], 1, $n)

u_min = @lift profile_min(u[$iter], 1, $n)
v_min = @lift profile_min(v[$iter], 1, $n)
T_min = @lift profile_min(T[$iter], 1, $n)
e_min = @lift profile_min(e[$iter], 1, $n)

color = (:pink, 0.4)
band!(ax_T, T_max, T_min; color)
band!(ax_e, e_max, e_min; color)
band!(ax_u, u_max, u_min; color)
band!(ax_u, v_max, v_min; color)
=#

for k = 1:Nensemble
    uk = @lift interior(u[$iter][$n])[k, 1, :]
    vk = @lift interior(v[$iter][$n])[k, 1, :]
    Tk = @lift interior(T[$iter][$n])[k, 1, :]
    ek = @lift interior(e[$iter][$n])[k, 1, :]

    lines!(ax_T, Tk, z, color=(:royalblue1, 0.8), linewidth=1,   label="CATKE")
    lines!(ax_u, uk, z, color=(:royalblue1, 0.8), linewidth=1,   label="u, CATKE")
    lines!(ax_u, vk, z, color=(:orange,     0.8), linewidth=1, label="v, CATKE")
    lines!(ax_e, ek, z, color=(:royalblue1, 0.8), linewidth=1,   label="CATKE")
end

u_truth = @lift interior(observations.field_time_serieses.u[$n])[1, 1, :]
v_truth = @lift interior(observations.field_time_serieses.v[$n])[1, 1, :]
T_truth = @lift interior(observations.field_time_serieses.T[$n])[1, 1, :]
e_truth = @lift interior(observations.field_time_serieses.e[$n])[1, 1, :]

lines!(ax_T, T_truth, z, color=(:gray23,  0.4),  linewidth=6, label="LES")
lines!(ax_u, u_truth, z, color=(:gray23,  0.4),  linewidth=6, label="u, LES")
lines!(ax_u, v_truth, z, color=(:darkred, 0.4),  linewidth=6, label="v, LES")
lines!(ax_e, e_truth, z, color=(:gray23,  0.4),  linewidth=6, label="LES")

xlims!(ax_u, -0.2, 0.3)
xlims!(ax_e, -0.0002, 0.0016)
xlims!(ax_T, 19.6, 20.0)

axislegend(ax_T, position=:rb, merge=true)
axislegend(ax_u, position=:rb, merge=true)
axislegend(ax_e, position=:rb, merge=true)

title = @lift @sprintf("LES and CATKE solutions at %.1f hours with iteration %d parameters", u[1].times[$n] / hour, $iter - 1)
Label(fig[0, :], title)

display(fig)

nswitch = floor(Int, Nt/length(output_paths))

#record(fig, "calibration_demo.mp4", 1:Nt; framerate=16) do nn
#    @info "Drawing frame $nn of $Nt..."
#    if nn % nswitch == 0 && iter.val < 2
#        iter[] += iter.val
#    end
#
#    n[] = nn
#end


text:
using XGBoost
using FactCheck

include("utils.jl")

facts("Sparse matrices") do
    X = sparse(randn(100,10) .* bitrand(100,10))
    y = randn(100)
    DMatrix(X, label = y)

    X = sparse(convert(Matrix{Float32}, randn(10,100) .* bitrand(10,100)))
    y = randn(100)
    DMatrix(X, true)

    X = sparse(randn(100,10) .* bitrand(100,10))
    y = randn(100)
    DMatrix(X)
end

facts("DMatrix loading") do
    dtrain = DMatrix("../data/agaricus.txt.train")
    train_X, train_Y = readlibsvm("../data/agaricus.txt.train", (6513, 126))
    @fact dtrain --> not(nothing)

    labels = get_info(dtrain, "label")

    @fact train_Y --> labels
end

facts("Agaricus training") do
    dtrain = DMatrix("../data/agaricus.txt.train")
    dtest = DMatrix("../data/agaricus.txt.test")
    watchlist = [(dtest, "eval"), (dtrain, "train")]

    bst = xgboost(dtrain, 2, watchlist=watchlist, eta = 1, max_depth = 2,
                  objective = "binary:logistic", silent = 1)
    @fact bst --> not(nothing)

    preds = XGBoost.predict(bst, dtest)

    labels = get_info(dtest, "label")
    @fact size(preds) --> size(labels)

    err = countnz((preds .> 0.5) .!= labels) / length(preds)
    @fact err --> less_than(0.1)
end

facts("Cross validation") do
    dtrain = DMatrix("../data/agaricus.txt.train")
    dtest = DMatrix("../data/agaricus.txt.test")
    watchlist = [(dtest, "eval"), (dtrain, "train")]

    bst = nfold_cv(dtrain, 5, 3, eta = 1, max_depth = 2, objective = "binary:logistic", silent = 1,
                   seed = 12345)
    # important_features = importance(bst)
    #
    # @fact startswith(important_features[1].fname, "f28") --> true
    # @pending important_features[1].fname --> "f28"
end

facts("Feature importance") do
    dtrain = DMatrix("../data/agaricus.txt.train")
    dtest = DMatrix("../data/agaricus.txt.test")
    watchlist = [(dtest, "eval"), (dtrain, "train")]

    bst = xgboost(dtrain, 5, watchlist = watchlist, eta = 1, max_depth = 2,
                  objective = "binary:logistic", silent = 1, seed = 12345)
    important_features = importance(bst)

    @fact startswith(important_features[1].fname, "f28") --> true
    @pending important_features[1].fname --> "f28"
end

facts("Example is running") do
    include("example.jl")
end