text:
# TODO: Remove some T's from output type signatures

#' @description
#'
#' A type that extends the normal Julia Array{T} type by allowing
#' the entries to be missing (i.e. NULL or NA) in addition to
#' taking on values of type `T`.
#'
#' @field data::Array{T, N} The elements of the DataArray are stored here,
#'        unless the given element is missing. In this case, the contents
#'        of the corresponding entry in `data` can be any arbitrary value.
#' @field na::BitArray{N} A 1-bit per entry Boolean mask specifying,
#'        for each entry of the DataArray, whether that entry is `NA`
#'        or not.
#'
#' @examples
#'
#' dv = DataArray([1, 2, 3], [false, false, true])
#'
#' dm = DataArray([1 2; 3 4], [false false; true false])
type DataArray{T, N} <: AbstractDataArray{T, N}
    data::Array{T, N}
    na::BitArray{N}

    # Ensure data values and missingness metadata match
    function DataArray(d::Array{T, N}, m::BitArray{N})
        if size(d) != size(m)
            msg = "Data and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        new(d, m)
    end
end

#' @description
#'
#' An DataVector is an DataArray of order 1.
typealias DataVector{T} DataArray{T, 1}

#' @description
#'
#' An DataMatrix is an DataArray of order 2.
typealias DataMatrix{T} DataArray{T, 2}

#' @description
#'
#' Create a DataArray from a set of data values and missingness
#' mask as an `Array{Bool}`.
#'
#' NB: This definition exists because Julia requires that we redefine
#' inner DataArray constructor as an outer constuctor.
#'
#' @param d::Array The non-NA values of the DataArray being constructed.
#' @param m::BitArray A Boolean mask indicating whether each element of
#'        the DataArray is missing or not.
#'
#' @returns out::DataArray A DataArray based on the input data and
#'          missingness mask.
#'
#' @examples
#'
#' da = DataArray([1, 2, 3], falses(3))
function DataArray{T, N}(d::Array{T, N},
                         m::BitArray{N} = falses(size(d))) # -> DataArray{T}
    return DataArray{T, N}(d, m)
end

#' @description
#'
#' Create a DataArray from a set of data values and missingness
#' mask as an `Array{Bool}`.
#'
#' @param d::Array The non-NA values of the DataArray being constructed.
#' @param m::Array{Bool} A Boolean mask indicating whether each element of
#'        the DataArray is missing or not.
#'
#' @returns out::DataArray A DataArray based on the input data and
#'          missingness mask.
#'
#' @examples
#'
#' da = DataArray([1, 2, 3], [false, false, true])
function DataArray(d::Array, m::Array{Bool}) # -> DataArray{T}
    return DataArray(d, bitpack(m))
end

#' @description
#'
#' Create a DataArray of a given type and dimensionality. All
#' entries will be set to `NA` by default.
#'
#' @param T::Type The type of the output DataArray.
#' @param dims::Integer... The size of the output DataArray.
#'
#' @returns out::DataArray A DataArray of the desired type and size.
#'
#' @examples
#'
#' da = DataArray(Int, 2, 2)
function DataArray(T::Type, dims::Integer...) # -> DataArray{T}
    return DataArray(Array(T, dims...), trues(dims...))
end

#' @description
#'
#' Create a DataArray of a given type and dimensionality. All
#' entries will be set to `NA` by default.
#'
#' @param T::Type The type of the output DataArray.
#' @param dims::NTuple{N, Int} The size of the output DataArray.
#'
#' @returns out::DataArray A DataArray of the desired type and size.
#'
#' @examples
#'
#' da = DataArray(Int, (2, 2))
function DataArray{N}(T::Type, dims::NTuple{N, Int}) # -> DataArray{T}
    return DataArray(Array(T, dims...), trues(dims...))
end

#' @description
#'
#' Create a copy of a DataArray.
#'
#' @param da::DataArray The DataArray that will be copied.
#'
#' @returns out::DataArray A copy of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' dv_new = copy(dv)
function Base.copy(d::DataArray) # -> DataArray{T}
    return DataArray(copy(d.data), copy(d.na))
end

function Base.copy!(dest::DataArray, src::DataArray) # -> DataArray{T}
    copy!(dest.data, src.data)
    copy!(dest.na, src.na)
    dest
end

#' @description
#'
#' Create a deep copy of a DataArray.
#'
#' @param da::DataArray The DataArray that will be deep copied.
#'
#' @returns out::DataArray A deep-copy of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' dv_new = deepcopy(dv)
function Base.deepcopy(d::DataArray) # -> DataArray{T}
    return DataArray(deepcopy(d.data), deepcopy(d.na))
end

#' @description
#'
#' Append the elements of items to the end of d.
#'
#' @param da::DataArray{T,1} The DataArray to append elements to.
#' @param items::DataArray{T,1} Elements to append to da.
#'
#' @returns: out::DataArray{T,1} The modified data array.
#'
#' @examples
#'
#' x = @data [1, 2, 3]
#' y = @data [4, 5, 6]
#' append!(x, y)
function Base.append!{T}(da::DataArray{T,1}, items::DataArray{T,1})
    append!(da.data, items.data)
    append!(da.na, items.na)
    da
end

#' @description
#'
#' Append the elements of items to the end of d.
#'
#' @param da::DataArray{T,1} The DataArray to append elements to.
#' @param items::AbstractArray{T,1} Elements to append to da.
#'
#' @returns: out::DataArray{T,1} The modified data array.
#'
#' @examples
#'
#' x = @data [1, 2, 3]
#' y = [4, 5, 6]
#' append!(x, y)
function Base.append!{T}(da::DataArray{T,1}, items::AbstractArray{T,1})
    append!(da.data, items)
    oldn = length(da.na)
    resize!(da.na, oldn + length(items))
    da.na[oldn+1:end] = false
    da
end

#' @description
#'
#' Create a new DataArray{T} that is similar to an existing DataArray.
#'
#' @param da::DataArray DataArray based on which a new DataArray
#'        will be created.
#' @param T::Type The element type of the output DataArray.
#' @param dims::Dims The dimensionality of the output DataArray.
#'
#' @returns out::DataArray{T} A new DataArray with the desired type and
#'          dimensionality.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' dv_new = similar(dv, Float64, 2, 2, 2)
function Base.similar(da::DataArray, T::Type, dims::Dims) #-> DataArray{T}
    return DataArray(Array(T, dims), trues(dims))
end

#' @description
#'
#' Find the sizes along each dimension of the DataArray.
#'
#' @param da::DataArray{Bool} DataArray whose size is desired.
#'
#' @returns a::Array{Int} Array containing the indices of all `true` values
#'          in `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' inds = find(dv)
Base.size(d::DataArray) = size(d.data) # -> (Int...)

#' @description
#'
#' Determine the number of dimensions (i.e. order) of a DataArray.
#'
#' @param da::DataArray{Bool} DataArray whose dimensionality is desired.
#'
#' @returns d::Int The number of dimensions of `da`.
#'
#' @examples
#'
#' dm = @data [false false; true false]
#' inds = ndims(dm)
Base.ndims(da::DataArray) = ndims(da.data) # -> Int

#' @description
#'
#' Determine the length (e.g. number of elements) of a DataArray.
#'
#' @param da::DataArray The DataArray whose length is desired.
#'
#' @returns n::Int The length of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' i = length(dv)
Base.length(d::DataArray) = length(d.data) # -> Int

#' @description
#'
#' Find the index of the last element of `da`.
#'
#' @param da::DataArray The DataArray whose final element's index is desired.
#'
#' @returns i::Int The index of the final element of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' i = endof(dv)
Base.endof(da::DataArray) = endof(da.data) # -> Int

#' @description
#'
#' Find the indices of all `true` values in a `DataArray{Bool}`.
#'
#' @param da::DataArray{Bool} DataArray whose `true` values will be found.
#'
#' @returns a::Array{Int} Array containing the indices of all `true` values
#'          in `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' inds = find(dv)
function Base.find(da::DataArray{Bool}) # -> Array{Int}
    data = da.data
    ntrue = 0
    @inbounds @bitenumerate da.na i na begin
        ntrue += !na && data[i]
    end
    res = Array(Int, ntrue)
    count = 1
    @inbounds @bitenumerate da.na i na begin
        if !na && data[i]
            res[count] = i
            count += 1
        end
    end
    return res
end

#' @description
#'
#' Turn a DataArray into an Array. Raises an error if NA's are encountered.
#'
#' @param da::DataArray{T} DataArray that will be converted to an Array.
#'
#' @returns a::Array{T} Array containing values of `da`.
#'
#' @examples
#'
#' dv = @data [1, 2, 3, 4]
#' v = array(dv)
#'
#' dm = @data [1 2; 3 4]
#' m = array(dm)
function array{T}(da::DataArray{T}) # -> Array{T}
    res = Array(T, size(da))
    for i in 1:length(da)
        if da.na[i]
            throw(NAException())
        else
            res[i] = da.data[i]
        end
    end
    return res
end

#' @description
#'
#' Turn a DataArray into an Array. Replace any NA's with the value
#' of second argument, `replacement`.
#'
#' @param da::DataArray{T} DataArray that will be converted to an Array.
#' @param replacement::T Value that will replace NA's in `da`.
#'
#' @returns a::Array{T} Array containing values of `da` plus replacements.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = array(dv, 3)
#'
#' dm = @data [1 2; NA 4]
#' m = array(dm, 3)
function array{T}(da::DataArray{T}, replacement::T) # -> Array{T}
    res = Array(T, size(da))
    for i in 1:length(da)
        if da.na[i]
            res[i] = replacement
        else
            res[i] = da.data[i]
        end
    end
    return res
end

#' @description
#'
#' Turn a DataArray into an Array. Replace any NA's with the value
#' of second argument, `replacement`.
#'
#' @param da::DataArray{T} DataArray that will be converted to an Array.
#' @param replacement::Any Value that will replace NA's in `da`.
#'        Converted to the `eltype`, `T`, of `da`.
#'
#' @returns a::Array{T} Array containing values of `da` plus replacements.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = array(dv, 3)
#'
#' dm = @data [1 2; NA 4]
#' m = array(dm, 3)
function array{T}(da::DataArray{T}, replacement::Any) # -> Array{T}
    return array(da, convert(T, replacement))
end

#' @description
#'
#' NO-OP: Turn a Vector into a Vector. See dropna(dv::DataVector) for
#'        rationale.
#'
#' @param v::Vector{T} Vector that will be converted to a Vector.
#'
#' @returns v::Vector{T} Vector containing all of the values of `v`.
#'
#' @examples
#'
#' v = [1, 2, 3, 4]
#' v = dropna(v)
dropna(v::AbstractVector) = copy(v) # -> AbstractVector

#' @description
#'
#' Turn a DataVector into a Vector. Drop any NA's.
#'
#' NB: Because NA's are dropped instead of replaced, this function only
#'     works on DataVector's and will not work on DataArray's of higher
#'     order.
#'
#' @param dv::DataVector{T} DataArray that will be converted to an Array.
#'
#' @returns v::Array{T} Array containing only the non-NA values of `dv`.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = dropna(dv)
dropna(dv::DataVector) = copy(dv.data[!dv.na]) # -> Vector

# Iterators
# TODO: Use values()
#       Use DataValueIterator type?

type EachFailNA{T}
    da::AbstractDataArray{T}
end
each_failNA{T}(da::AbstractDataArray{T}) = EachFailNA(da)
Base.start(itr::EachFailNA) = 1
function Base.done(itr::EachFailNA, ind::Integer)
    return ind > length(itr.da)
end
function Base.next(itr::EachFailNA, ind::Integer)
    if isna(itr.da[ind])
        throw(NAException())
    else
        (itr.da[ind], ind + 1)
    end
end

type EachDropNA{T}
    da::AbstractDataArray{T}
end
each_dropna{T}(da::AbstractDataArray{T}) = EachDropNA(da)
Base.start(itr::EachDropNA) = 1
function Base.done(itr::EachDropNA, ind::Integer)
    return ind > length(itr.da)
end
function Base.next(itr::EachDropNA, ind::Integer)
    while ind <= length(itr.da) && isna(itr.da[ind])
        ind += 1
    end
    (itr.da[ind], ind + 1)
end

type EachReplaceNA{S, T}
    da::AbstractDataArray{S}
    replacement_val::T
end
function each_replaceNA(da::AbstractDataArray, val::Any)
    EachReplaceNA(da, convert(eltype(da), val))
end
function each_replaceNA(val::Any)
    x -> each_replaceNA(x, val)
end
Base.start(itr::EachReplaceNA) = 1
function Base.done(itr::EachReplaceNA, ind::Integer)
    return ind > length(itr.da)
end
function Base.next(itr::EachReplaceNA, ind::Integer)
    if isna(itr.da[ind])
        (itr.replacement_val, ind + 1)
    else
        (itr.da[ind], ind + 1)
    end
end

# Indexing

typealias SingleIndex Real
# TODO: Remove BitVector here
typealias MultiIndex Union(Vector, BitVector, Ranges, Range1)
typealias BooleanIndex Union(BitVector, Vector{Bool})

# TODO: Solve ambiguity warnings here without
#       ridiculous accumulation of methods

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param v::Vector A Vector whose elements will be retrieved.
#' @param inds::AbstractDataVector{Bool} An AbstractDataVector
#'        containing a Boolean mask specifying, for each element,
#'        whether that element should be retrieved or not. `NA`
#'        values in the mask are treated as `false`.
#'
#' @returns out::Vector The values of `v` at the requested indices.
#'
#' @examples
#'
#' v = [1, 2, 3]
#' inds = @data([false, true])
#' v[inds]
function Base.getindex(v::Vector,
                       inds::AbstractDataVector{Bool})
    return v[find(inds)]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param v::Vector A Vector whose elements will be retrieved.
#' @param inds::AbstractDataArray{Bool} A DataArray containing
#'        a Boolean mask specifying, for each element, whether
#'        that element should be retrieved or not. `NA` values
#'        in the mask are treated as `false`.
#'
#' @returns out::Vector The values of `v` at the requested indices.
#'
#' @examples
#'
#' v = [1, 2, 3]
#' inds = @data([false, true])
#' v[inds]
function Base.getindex(v::Vector,
                       inds::AbstractDataArray{Bool})
    return v[find(inds)]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param a::Array An Array whose elements will be retrieved.
#' @param inds::AbstractDataArray{Bool} A DataArray containing
#'        a Boolean mask specifying, for each element, whether
#'        that element should be retrieved or not. `NA` values
#'        in the mask are treated as `false`.
#'
#' @returns out::Vector{S} The values of `a` at the requested indices.
#'
#' @examples
#'
#' a = [1, 2, 3]
#' inds = @data([false, true])
#' a[inds]
function Base.getindex(a::Array,
                       inds::AbstractDataArray{Bool})
    return a[find(inds)]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param a::Vector{S} A Vector whose elements will be retrieved.
#' @param inds::AbstractDataArray{T} A DataArray containing the indices
#'        of the requested elements.
#'
#' @returns out::Vector{S} The values of `a` at the requested indices.
#'
#' @examples
#'
#' a = [1, 2, 3]
#' inds = @data([1, 2])
#' da[1]
#'
#' da = @data([NA, 2, 3])
#' da[1]
function Base.getindex{S, T}(a::Vector{S},
                             inds::AbstractDataArray{T})
    return a[dropna(inds)]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param a::Array{T, N} An Array whose elements will be retrieved.
#' @param inds::AbstractDataArray{T} A DataArray containing the indices
#'        of the requested elements.
#'
#' @returns out::Array{T} The values of `a` at the requested indices.
#'
#' @examples
#'
#' a = [1, 2, 3]
#' inds = @data([1, 2])
#' da[1]
#'
#' da = @data([NA, 2, 3])
#' da[1]
function Base.getindex{S, T}(a::Array{S},
                             inds::AbstractDataArray{T})
    return a[dropna(inds)]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be retrieved.
#' @param ind::Real The index of the element to be retrieved.
#'
#' @returns out::Union(T, NAtype) The value of `da` at the requested
#'          index.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[1]
#'
#' da = @data([NA, 2, 3])
#' da[1]
function Base.getindex(da::DataArray, ind::SingleIndex)
	if da.na[ind]
		return NA
	else
		return da.data[ind]
	end
end

# d[MultiItemIndex]
# TODO: Return SubDataArray
#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be retrieved.
#' @param inds::AbstractDataVector{Bool} A Boolean mask specifying for
#'        each element whether that element will be retrieved or not.
#'
#' @returns out::DataArray{T} A DataArray containing the elements of
#'          `da` indexed.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' inds = @data([true, false])
#' da[inds]
function Base.getindex(d::DataArray,
                       inds::AbstractDataVector{Bool})
    inds = find(inds)
    return d[inds]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be retrieved.
#' @param inds::AbstractDataVector The indices of the elements that
#'        will be retrieved.
#'
#' @returns out::DataArray{T} A DataArray containing the elements of
#'          `da` indexed.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' inds = @data([1, 2])
#' da[inds]
function Base.getindex(d::DataArray,
                       inds::AbstractDataVector)
    inds = dropna(inds)
    return d[inds]
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be retrieved.
#' @param inds::BooleanIndex A Boolean mask specifying whether each
#'        element of `da` should be retrieved.
#'
#' @returns out::DataArray{T} A DataArray containing the elements of
#'          `da` indexed.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[true, false]]
function Base.getindex(d::DataArray, inds::BooleanIndex)
    res = similar(d, sum(inds))
    j = 1
    for i in 1:length(inds)
        if inds[i]
            if !d.na[i]
                res[j] = d.data[i]
            end
            j += 1
        end
    end
    return res
end

#' @description
#'
#' Get a set of elements of a DataArray.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be retrieved.
#' @param inds::MultiIndex The indices of the elements of `da` to be
#'        retrieved.
#'
#' @returns out::DataArray{T} A DataArray containing the elements of
#'          `da` indexed.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[1, 2]]
function Base.getindex(d::DataArray, inds::MultiIndex)
    res = similar(d, length(inds))
    for i in 1:length(inds)
        ix = inds[i]
        if !d.na[ix]
            res[i] = d.data[ix]
        end
    end
    return res
end

#' @description
#'
#' Set one element of a DataArray to `NA`.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be modified.
#' @param val::NAtype The `NA` value being assigned to an element of `da`.
#' @param ind::Real A real value specifying the index of the element
#'        of `da` being modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[1] = NA
function Base.setindex!(da::DataArray,
                        val::NAtype,
                        i::SingleIndex) # -> NAtype
	da.na[i] = true
    return da
end

#' @description
#'
#' Set one element of a DataArray to a new value, `val`.
#'
#' @param da::DataArray{T, N} A DataArray whose element will be modified.
#' @param val::Any The value being assigned to an element of `da`.
#' @param ind::Real A real value specifying the index of the element
#'        of `da` being modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[1] = 4
function Base.setindex!(da::DataArray,
                        val::Any,
                        ind::SingleIndex) # -> Any
	da.data[ind] = val
	da.na[ind] = false
    return da
end

#' @description
#'
#' Set a specified set of elements of a DataArray to `NA`.
#'
#' NB: The indices in `inds` must be exhaustive.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param val::NAtype The NA value being assigned to elements of `da`.
#' @param inds::AbstractVector{Bool} A Boolean vector specifying for every
#'        element of `da` whether it will be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[true, false, true]] = NA
function Base.setindex!(da::DataArray,
                        val::NAtype,
                        inds::AbstractVector{Bool}) # -> NAtype
    da.na[find(inds)] = true
    return da
end

#' @description
#'
#' Set a specified set of elements of a DataArray to `NA`.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param val::NAtype The NA value being assigned to elements of `da`.
#' @param inds::AbstractVector A vector of indices of `da` to
#'        be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[1:3] = NA
function Base.setindex!(da::DataArray,
                        val::NAtype,
                        inds::AbstractVector) # -> NAtype
    da.na[inds] = true
    return da
end

#' @description
#'
#' Determine if the entries of an DataArray are `NA`.
#'
#' @param a::DataArray{T, N} The DataArray whose missingness will
#'        be assessed.
#'
#' @returns na::BitArray{N} Elementwise Boolean whether entry is missing.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' isna(da)
isna(da::DataArray) = copy(da.na) # -> BitArray

#' @description
#'
#' Safe and type-stable way to determine if element `i` of an
#' DataArray is `NA`.
#'
#' @param a::DataArray The DataArray whose missingness will
#'        be assessed.
#' @param inds::Any The indices of the elements to be checked for `NA`.
#'
#' @returns na::Any Are the indexed elements `NA` or not?
#'
#' @examples
#'
#' a = @data([1, 2, 3])
#' isna(a, 1)
#' isna(a, 1:2)
isna(da::DataArray, inds::Any...) = getindex(da.na, inds...)

#' @description
#'
#' Determine if the entries of an DataArray are `NaN`.
#'
#' @param a::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns na::DataArray{Bool} Elementwise Boolean whether entry is `NaN`.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' isnan(da)
function Base.isnan(da::DataArray) # -> DataArray{Bool}
    return DataArray(isnan(da.data), copy(da.na))
end

#' @description
#'
#' Determine if the entries of an DataArray are finite, which means
#' neither `+/-NaN` nor `+/-Inf`.
#'
#' @param a::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns na::DataArray{Bool} Elementwise Boolean whether entry is finite.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' isfinite(da)
function Base.isfinite(da::DataArray) # -> DataArray{Bool}
    n = length(da)
    res = Array(Bool, size(da))
    for i in 1:n
        if !da.na[i]
            res[i] = isfinite(da.data[i])
        end
    end
    return DataArray(res, copy(da.na))
end

# Promotion rules

# promote_rule{T, T}(::Type{AbstractDataArray{T}},
#                    ::Type{T}) = promote_rule(T, T)
# promote_rule{S, T}(::Type{AbstractDataArray{S}},
#                    ::Type{T}) = promote_rule(S, T)
# promote_rule{T}(::Type{AbstractDataArray{T}}, ::Type{T}) = T

#' @description
#'
#' Convert a DataArray{T} to an Array{S}. Throws an `NAException`
#' if the input contains `NA` values that prohibit conversion.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns a::Array{S} The (possibly type-converted) elements of
#'         `da` if none were `NA`.
#'
#' @examples
#'
#' da = @data [1 2; 3 NA]
#' a = convert(Array{Float64}, da)
#'
#' da = @data [1 2; 3 4]
#' a = convert(Array{Float64}, da)
function Base.convert{S, T, N}(::Type{Array{S, N}},
                               x::DataArray{T, N}) # -> Array{S, N}
    if anyna(x)
        err = "Cannot convert DataArray with NA's to desired type"
        throw(NAException(err))
    else
        return convert(Array{S, N}, x.data)
    end
end

#' @description
#'
#' Convert a DataArray{T} to an Array{T}. Throws an `NAException`
#' if the input contains `NA` values that prohibit conversion.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns a::Array{T} The elements of `da` if none were `NA`.
#'
#' @examples
#'
#' da = @data [1 2; 3 NA]
#' a = convert(Array, da)
#'
#' da = @data [1 2; 3 4]
#' a = convert(Array, da)
# TODO: Consider making a copy here
function Base.convert{T, N}(::Type{Array},
                            x::DataArray{T, N}) # -> Array{T}
    if anyna(x)
        err = "Cannot convert DataArray with NA's to base type"
        throw(NAException(err))
    else
        return x.data
    end
end

#' @description
#'
#' Convert an Array{T} to a DataArray{S}.
#'
#' @param a::Array{T} The Array that will be converted.
#'
#' @returns da::DataArray{S} The converted DataArray with potential
#'          type-conversion of the elements of `a`.
#'
#' @examples
#'
#' a = [1 2; 3 4]
#' da = convert(DataArray{Float64}, a)
function Base.convert{S, T, N}(::Type{DataArray{S, N}},
                               a::AbstractArray{T, N}) # -> DataArray{S, N}
    return DataArray(convert(Array{S, N}, a), falses(size(a)))
end

#' @description
#'
#' Convert an Array{T} to a DataArray{T}.
#'
#' @param a::Array{T} The Array that will be converted.
#'
#' @returns da::DataArray{T} The converted DataArray.
#'
#' @examples
#'
#' a = [1 2; 3 4]
#' da = convert(DataArray, a)
function Base.convert{T, N}(::Type{DataArray},
                            a::AbstractArray{T, N}) # -> DataArray{T, N}
    return DataArray(convert(Array{T, N}, a), falses(size(a)))
end

#' @description
#'
#' Convert a DataArray{T} to a DataArray{S}.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns out::DataArray{S} The converted DataArray.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' dv_alt = convert(DataVector{Float64}, dv)
function Base.convert{S, T, N}(::Type{DataArray{S, N}},
                               x::DataArray{T, N}) # -> DataArray{S, N}
    return DataArray(convert(Array{S}, x.data), x.na)
end

#' @description
#'
#' NO-OP: See convert(DataArray{S}, DataArray{T}) for rationale.
#' TODO: Make operative by doing copy?
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns out::DataArray{T} The converted DataArray.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' dv_alt = convert(DataVector, dv)
function Base.convert{T, N}(::Type{DataArray},
                            x::DataArray{T, N}) # -> DataArray{T, N}
    return DataArray(x.data, x.na)
end

#' @description
#'
#' Convert an AbstractArray to a DataArray.
#'
#' @param a::AbstractArray{T} The AbstractArray that will be converted to
#'        a DataArray.
#'
#' @returns da::DataArray{T} `a` after conversion to a DataArray.
#'
#' @examples
#'
#' dv = data([1, 2, 3, 4])
#' dm = data([1 2; 3 4])
#' dv = data(1:10)
#' dv = data(falses(3))

data(a::AbstractArray) = convert(DataArray, a)

#' @description
#'
#' Convert a DataArray to an Array of int, float or bool type.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns a::Array{Union(Int, Float64, Bool)} An Array containing the
#'          type-converted values of `da`.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = int(dv)
#' v = float(dv)
#' v = bool(dv)
#
# TODO: Make sure these handle copying correctly
# TODO: Remove these? They have odd behavior, because they convert to Array's.
# TODO: Rethink multi-item documentation approach
for f in (:(Base.int), :(Base.float), :(Base.bool))
    @eval begin
        function ($f)(da::DataArray) # -> DataArray
            if anyna(da)
                err = "Cannot convert DataArray with NA's to desired type"
                throw(NAException(err))
            else
                ($f)(da.data)
            end
        end
    end
end

#' @description
#'
#' Compute the hash of an AbstractDataArray.
#'
#' @param da::DataArray{T} DataArray whose hash is desired.
#'
#' @returns h::Uint An unsigned integer hash value for `da`.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' k = hash(dv)
#
# TODO: Make sure this agrees with is_equals()
function Base.hash(a::AbstractDataArray) # -> Uint
    h = hash(size(a)) + 1
    for i in 1:length(a)
        h = bitmix(h, int(hash(a[i])))
    end
    return uint(h)
end

#' @internal
#' @description
#'
#' Find the unique values in a DataArray, noting if `NA` occurs in the
#' DataArray.
#'
#' @param da::DataArray{T} DataArray whose unique values are desired.
#'
#' @returns v::Vector{T} Vector containing the unique values from `da`.
#' @returns hasna::Bool Did `da` contain any `NA` entries?
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' distinct_values, hasna = finduniques(dv)
function finduniques{T}(da::DataArray{T}) # -> Vector{T}, Bool
    unique_values = Dict{T, Bool}()
    n = length(da)
    hasna = false
    for i in 1:n
        if da.na[i]
            hasna = true
        else
            unique_values[da.data[i]] = true
        end
    end
    return unique_values, hasna
end

#' @description
#'
#' Return a DataVector containing the unique values of a DataArray,
#' including NA if it is encountered.
#'
#' @param da::DataArray{T} DataArray whose unique values are desired.
#'
#' @returns dv::DataVector{T} DataVector containing the unique values
#'          from `da`, including NA if there are any missing entries
#'          in `da`.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' distinct_values = unique(dv)
function Base.unique{T}(da::DataArray{T}) # -> DataVector{T}
    unique_values, hasna = finduniques(da)
    n = length(unique_values)
    if hasna
        res = DataArray(Array(T, n + 1))
        i = 0
        for val in keys(unique_values)
            i += 1
            res.data[i] = val
        end
        res.na[n + 1] = true
        return res
    else
        return DataArray(collect(keys(unique_values)))
    end
end

#' @description
#'
#' Return a DataVector containing the unique values of a DataArray,
#' excluding any NA's.
#'
#' @param da::DataArray{T} DataArray whose unique values are desired.
#'
#' @returns dv::DataVector{T} DataVector containing the unique values
#'          from `da`, excluding any NA's.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' distinct_values = levels(dv)
function levels(da::DataArray) # -> DataVector{T}
    unique_values, hasna = finduniques(da)
    return DataArray(collect(keys(unique_values)))
end


text:
Base.:*(x::Union{Num, Symbolic}, y::Unitful.AbstractQuantity) = x * y
Base.:/(x::Union{Num, Symbolic}, y::Unitful.AbstractQuantity) = x / y

struct ValidationError <: Exception
    message::String
end

"Throw exception on invalid unit types, otherwise return argument."
function screen_unit(result)
    result isa Unitful.Unitlike ||
        throw(ValidationError("Unit must be a subtype of Unitful.Unitlike, not $(typeof(result))."))
    result isa Unitful.ScalarUnits ||
        throw(ValidationError("Non-scalar units such as $result are not supported. Use a scalar unit instead."))
    result == u"°" &&
        throw(ValidationError("Degrees are not supported. Use radians instead."))
    result
end

"""Test unit equivalence.

Example of implemented behavior:
```julia
using ModelingToolkit, Unitful
MT = ModelingToolkit
@parameters γ P [unit = u"MW"] E [unit = u"kJ"] τ [unit = u"ms"]
@test MT.equivalent(u"MW" ,u"kJ/ms") # Understands prefixes
@test !MT.equivalent(u"m", u"cm") # Units must be same magnitude
@test MT.equivalent(MT.get_unit(P^γ), MT.get_unit((E/τ)^γ)) # Handles symbolic exponents
```
"""
equivalent(x, y) = isequal(1 * x, 1 * y)
unitless = Unitful.unit(1)

#For dispatching get_unit
Literal = Union{Sym, Symbolics.ArrayOp, Symbolics.Arr, Symbolics.CallWithMetadata}
Conditional = Union{typeof(ifelse), typeof(IfElse.ifelse)}
Comparison = Union{typeof.([==, !=, ≠, <, <=, ≤, >, >=, ≥])...}

"Find the unit of a symbolic item."
get_unit(x::Real) = unitless
get_unit(x::Unitful.Quantity) = screen_unit(Unitful.unit(x))
get_unit(x::AbstractArray) = map(get_unit, x)
get_unit(x::Num) = get_unit(value(x))
get_unit(x::Literal) = screen_unit(getmetadata(x, VariableUnit, unitless))
get_unit(op::Differential, args) = get_unit(args[1]) / get_unit(op.x)
get_unit(op::Difference, args) = get_unit(args[1]) / get_unit(op.t)
get_unit(op::typeof(getindex), args) = get_unit(args[1])
get_unit(x::SciMLBase.NullParameters) = unitless

function get_unit(op, args) # Fallback
    result = op(1 .* get_unit.(args)...)
    try
        unit(result)
    catch
        throw(ValidationError("Unable to get unit for operation $op with arguments $args."))
    end
end

function get_unit(op::Integral, args)
    unit = 1
    if op.domain.variables isa Vector
        for u in op.domain.variables
            unit *= get_unit(u)
        end
    else
        unit *= get_unit(op.domain.variables)
    end
    return get_unit(args[1]) * unit
end

function get_unit(x::Pow)
    pargs = arguments(x)
    base, expon = get_unit.(pargs)
    @assert expon isa Unitful.DimensionlessUnits
    if base == unitless
        unitless
    else
        pargs[2] isa Number ? base^pargs[2] : (1 * base)^pargs[2]
    end
end

function get_unit(x::Add)
    terms = get_unit.(arguments(x))
    firstunit = terms[1]
    for other in terms[2:end]
        termlist = join(map(repr, terms), ", ")
        equivalent(other, firstunit) ||
            throw(ValidationError(", in sum $x, units [$termlist] do not match."))
    end
    return firstunit
end

function get_unit(op::Conditional, args)
    terms = get_unit.(args)
    terms[1] == unitless ||
        throw(ValidationError(", in $x, [$(terms[1])] is not dimensionless."))
    equivalent(terms[2], terms[3]) ||
        throw(ValidationError(", in $x, units [$(terms[2])] and [$(terms[3])] do not match."))
    return terms[2]
end

function get_unit(op::typeof(Symbolics._mapreduce), args)
    if args[2] == +
        get_unit(args[3])
    else
        throw(ValidationError("Unsupported array operation $op"))
    end
end

function get_unit(op::Comparison, args)
    terms = get_unit.(args)
    equivalent(terms[1], terms[2]) ||
        throw(ValidationError(", in comparison $op, units [$(terms[1])] and [$(terms[2])] do not match."))
    return unitless
end

function get_unit(x::Symbolic)
    if SymbolicUtils.istree(x)
        op = operation(x)
        if op isa Sym || (op isa Term && operation(op) isa Term) # Dependent variables, not function calls
            return screen_unit(getmetadata(x, VariableUnit, unitless)) # Like x(t) or x[i]
        elseif op isa Term && !(operation(op) isa Term)
            gp = getmetadata(x, Symbolics.GetindexParent, nothing) # Like x[1](t)
            return screen_unit(getmetadata(gp, VariableUnit, unitless))
        end  # Actual function calls:
        args = arguments(x)
        return get_unit(op, args)
    else # This function should only be reached by Terms, for which `istree` is true
        throw(ArgumentError("Unsupported value $x."))
    end
end

"Get unit of term, returning nothing & showing warning instead of throwing errors."
function safe_get_unit(term, info)
    side = nothing
    try
        side = get_unit(term)
    catch err
        if err isa Unitful.DimensionError
            @warn("$info: $(err.x) and $(err.y) are not dimensionally compatible.")
        elseif err isa ValidationError
            @warn(info*err.message)
        elseif err isa MethodError
            @warn("$info: no method matching $(err.f) for arguments $(typeof.(err.args)).")
        else
            rethrow()
        end
    end
    side
end

function _validate(terms::Vector, labels::Vector{String}; info::String = "")
    valid = true
    first_unit = nothing
    first_label = nothing
    for (term, label) in zip(terms, labels)
        equnit = safe_get_unit(term, info * label)
        if equnit === nothing
            valid = false
        elseif !isequal(term, 0)
            if first_unit === nothing
                first_unit = equnit
                first_label = label
            elseif !equivalent(first_unit, equnit)
                valid = false
                @warn("$info: units [$(first_unit)] for $(first_label) and [$(equnit)] for $(label) do not match.")
            end
        end
    end
    valid
end

function validate(jump::Union{ModelingToolkit.VariableRateJump,
                              ModelingToolkit.ConstantRateJump}, t::Symbolic;
                  info::String = "")
    newinfo = replace(info, "eq." => "jump")
    _validate([jump.rate, 1 / t], ["rate", "1/t"], info = newinfo) && # Assuming the rate is per time units
        validate(jump.affect!, info = newinfo)
end

function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::String = "")
    left_symbols = [x[1] for x in jump.reactant_stoch] #vector of pairs of symbol,int -> vector symbols
    net_symbols = [x[1] for x in jump.net_stoch]
    all_symbols = vcat(left_symbols, net_symbols)
    allgood = _validate(all_symbols, string.(all_symbols); info)
    n = sum(x -> x[2], jump.reactant_stoch, init = 0)
    base_unitful = all_symbols[1] #all same, get first
    allgood && _validate([jump.scaled_rates, 1 / (t * base_unitful^n)],
              ["scaled_rates", "1/(t*reactants^$n))"]; info)
end

function validate(jumps::ArrayPartition{<:Union{Any, Vector{<:JumpType}}}, t::Symbolic)
    labels = ["in Mass Action Jumps,", "in Constant Rate Jumps,", "in Variable Rate Jumps,"]
    all([validate(jumps.x[idx], t, info = labels[idx]) for idx in 1:3])
end

function validate(eq::ModelingToolkit.Equation; info::String = "")
    _validate([eq.lhs, eq.rhs], ["left", "right"]; info)
end
function validate(eq::ModelingToolkit.Equation,
                  term::Union{Symbolic, Unitful.Quantity, Num}; info::String = "")
    _validate([eq.lhs, eq.rhs, term], ["left", "right", "noise"]; info)
end
function validate(eq::ModelingToolkit.Equation, terms::Vector; info::String = "")
    _validate(vcat([eq.lhs, eq.rhs], terms),
              vcat(["left", "right"], "noise  #" .* string.(1:length(terms))); info)
end

"Returns true iff units of equations are valid."
function validate(eqs::Vector; info::String = "")
    all([validate(eqs[idx], info = info * " in eq. #$idx") for idx in 1:length(eqs)])
end
function validate(eqs::Vector, noise::Vector; info::String = "")
    all([validate(eqs[idx], noise[idx], info = info * " in eq. #$idx")
         for idx in 1:length(eqs)])
end
function validate(eqs::Vector, noise::Matrix; info::String = "")
    all([validate(eqs[idx], noise[idx, :], info = info * " in eq. #$idx")
         for idx in 1:length(eqs)])
end
function validate(eqs::Vector, term::Symbolic; info::String = "")
    all([validate(eqs[idx], term, info = info * " in eq. #$idx") for idx in 1:length(eqs)])
end
validate(term::Symbolics.SymbolicUtils.Symbolic) = safe_get_unit(term, "") !== nothing

"Throws error if units of equations are invalid."
function check_units(eqs...)
    validate(eqs...) ||
        throw(ValidationError("Some equations had invalid units. See warnings for details."))
end
all_dimensionless(states) = all(x -> safe_get_unit(x, "") in (unitless, nothing), states)

ext: jl


text:
include("./util.jl")
include("./parsers.jl")
include("./Algorithms/TEAVAR_Star.jl")
using Dates;

env = Gurobi.Env()
setparam!(env, "Method", 2) # choose only barrier; since crossover is needed for concurrent and that takes too long
setparam!(env, "Crossover", 0) # disable barrier crossover

print("method= ", getparam(env, "Method"), "; crossover= ", getparam(env, "Crossover"))

topology = ARGS[1]
# failure probabilities must always come from topology
weibull = false
shape = 0.8
scale = 0.01

# we will change demand numerals and beta
demand_num = parse(Int, ARGS[2]) # used to be 1
beta=parse(Float64, ARGS[3]) # used to be 0.99
max_cf=ARGS[4] == "mcf"
paths = ARGS[5] # used to be "SMORE" should be "SMORE4" or "SMORE8"
what_to_read = ARGS[6]
cutoff_downscale = parse(Float64, ARGS[7]) # used to 10

# output file for result analysis and debugging
outputfile=string(topology, "_d", demand_num, "_beta", beta, "_mcf", max_cf, "_paths", paths, "_topo", what_to_read, "_cutoffDS", cutoff_downscale)

links, capacity, link_probs, nodes = readTopology(topology, what_to_read=what_to_read)
demand, flows = readDemand("$(topology)/demand", length(nodes), demand_num, matrix=true)
T, Tf, k = parsePaths("$(topology)/paths/$(paths)", links, flows)

if weibull
  probabilities = weibullProbs(length(links), shape=shape, scale=scale)
else
  probabilities = link_probs
end
println("FailureProbs= ", probabilities)

beginning_cutoff = parse(Float64, ARGS[8])

let cutoff = beginning_cutoff
	while true
		println(Dates.format(now(), "HH:MM:SS"), " going to subScenariosRecursion with cutoff=", cutoff)
		flush(stdout)

		task_cutoff_finder = @task global scenarios, scenario_probs = subScenariosRecursion(probabilities, cutoff)
		t = Timer(60)
		# run for 60s or until task cutoff finder finishes
		schedule(task_cutoff_finder)
		while (!@isdefined(scenarios) || length(scenarios) == 0 || isopen(t)) && !istaskdone(task_cutoff_finder)
			# println(Dates.format(now(), "HH::MM::SS"), ": testing")
			# flush(stdout)
			yield()
		end

		nscenarios = length(scenarios)
		total_scenario_prob = sum(scenario_probs)

		println(Dates.format(now(), "HH:MM:SS"), " cutoff =", cutoff, " #scenarios=", nscenarios, " total_scenar_prob=", total_scenario_prob)
		flush(stdout)

		if total_scenario_prob >= 1 - (1-beta)/ cutoff_downscale
			break
		end

		if !isopen(t)
			println(Dates.format(now(), "HH:MM:SS"), ": timed out!")
			flush(stdout)
			break
		end

		cutoff = cutoff / cutoff_downscale
	end
end

TEAVAR_Star(env, links, capacity, flows, demand, beta, k,
    T, Tf, scenarios, scenario_probs, outputfile, explain=true, verbose=true,
        utilization=true, max_concurrent_flow=max_cf)


text:
using GeometricIntegrators
using GeometricProblems.HarmonicOscillator
using SimpleSolvers
using Test

SimpleSolvers.set_config(:nls_atol_break, Inf)
SimpleSolvers.set_config(:nls_rtol_break, Inf)
SimpleSolvers.set_config(:nls_stol_break, Inf)

Δt = 0.1
nt = 10000
ns = 10

nsave1 = 1
nwrte1 = 100

nsave2 = 10
nwrte2 = 100

h5file = "test.hdf5"

tab = TableauImplicitMidpoint()
ode = harmonic_oscillator_ode([rand(2) for i in 1:ns])

# @testset "$(rpad("Serial Simulation",80))" begin

    sim1 = Simulation(ode, tab, Δt, "Harmonic Oscillator Test 1 (Serial)", h5file, nt; nsave=nsave1, nwrite=nwrte1)
    run!(sim1)
    @test isfile(h5file)
    ssol1 = SolutionODE(h5file)
    @test ssol1.ntime == nt
    @test ssol1.nsave == nsave1
    rm(h5file)

    sim2 = Simulation(ode, tab, Δt, "Harmonic Oscillator Test 2 (Serial)", h5file, nt; nsave=nsave2, nwrite=nwrte2)
    run!(sim2)
    @test isfile(h5file)
    ssol2 = SolutionODE(h5file)
    @test ssol2.ntime == nt
    @test ssol2.nsave == nsave2
    rm(h5file)

# end

# @testset "$(rpad("Parallel Simulation (" * string(Threads.nthreads()) * " Threads)",80))" begin

    sim1 = ParallelSimulation(ode, tab, Δt, "Harmonic Oscillator Test 1 (" * string(Threads.nthreads()) * " Threads)", h5file, nt; nsave=nsave1, nwrite=nwrte1)
    run!(sim1)
    @test isfile(h5file)
    psol1 = SolutionODE(h5file)
    @test psol1.ntime == nt
    @test psol1.nsave == nsave1
    rm(h5file)

    sim2 = ParallelSimulation(ode, tab, Δt, "Harmonic Oscillator Test 2 (" * string(Threads.nthreads()) * " Threads)", h5file, nt; nsave=nsave2, nwrite=nwrte2)
    run!(sim2)
    @test isfile(h5file)
    psol2 = SolutionODE(h5file)
    @test psol2.ntime == nt
    @test psol2.nsave == nsave2
    rm(h5file)

# end

@test psol1.q == ssol1.q
@test psol2.q == ssol2.q


text:
# TODO: Support Symbol ordering

"""
    Atomix.get(ref, order) -> x
    Atomix.get(ref) -> x

Atomically load the value `x` stored in `ref` with ordering `order`.  The default ordering
`Atomix.sequentially_consistent` is used when not specified.

# Examples
```jldoctest
julia> using Atomix

julia> a = [111, 222, 333];

julia> ref = Atomix.IndexableRef(a, (1,));

julia> Atomix.get(ref)
111
```
"""
Atomix.get

"""
    Atomix.set!(ref, x, order)
    Atomix.set!(ref, x)

Atomically store the value `x` in `ref` with ordering `order`.  The default ordering
`Atomix.sequentially_consistent` is used when not specified.

# Examples
```jldoctest
julia> using Atomix

julia> a = [111, 222, 333];

julia> ref = Atomix.IndexableRef(a, (1,));

julia> Atomix.set!(ref, 123);

julia> a[1]
123
```
"""
Atomix.set!

"""
    Atomix.modify!(ref, op, x, order) -> (old => new)
    Atomix.modify!(ref, op, x) -> (old => new)

Atomically update `ref` from stored value `old` to `new = op(old, x)` with ordering `order`
(default: `Atomix.sequentially_consistent`).  Return a pair `old => new`.

# Examples
```jldoctest
julia> using Atomix

julia> a = [111, 222, 333];

julia> ref = Atomix.IndexableRef(a, (1,));

julia> Atomix.modify!(ref, +, 123)
111 => 234
```
"""
Atomix.modify!

"""
    Atomix.swap!(ref, new, order) -> old
    Atomix.swap!(ref, new) -> old

Swap the `old` stored in `ref` with the `new` value and establish the memory ordering
`order` (default: `Atomix.sequentially_consistent`).

Notes for implementers: `Atomix.swap!(ref, new, order)` is defined as `Atomix.modify!(ref,
Atomix.right, x, order)`.  Thus, only `Atomix.modify!` has to be defined.

# Examples
```jldoctest
julia> using Atomix

julia> a = [111, 222, 333];

julia> ref = Atomix.IndexableRef(a, (1,));

julia> Atomix.swap!(ref, 123)
111

julia> a[1]
123
```
"""
Atomix.swap!

"""
    Atomix.replace!(ref, expected, desired, success_order, fail_order) -> (; old, success)
    Atomix.replace!(ref, expected, desired, order) -> (; old, success)
    Atomix.replace!(ref, expected, desired) -> (; old, success)

Atomically replace the value stored in `ref` to `desired` if `expected` is stored.  A named
tuple `(; old::eltype(ref), success::Bool)` is returned.

# Examples
```jldoctest
julia> using Atomix

julia> a = [111, 222, 333];

julia> ref = Atomix.IndexableRef(a, (1,));

julia> Atomix.replace!(ref, 111, 123)
(old = 111, success = true)
```
"""
Atomix.replace!

@inline Atomix.get(ref) = Atomix.get(ref, seq_cst)
@inline Atomix.set!(ref, x) = Atomix.set!(ref, x, seq_cst)

@inline Atomix.modify!(ref, op::OP, x) where {OP} = Atomix.modify!(ref, op, x, seq_cst)

@inline Atomix.replace!(ref, expected, desired, order::Ordering = seq_cst) =
    Atomix.replace!(ref, expected, desired, order, order)

@inline Atomix.swap!(ref, x, order::Ordering = seq_cst) =
    first(Atomix.modify!(ref, right, x, order))


text:
@inline function getVelocityDelta(joint::Rotational2, body1::AbstractBody, body2::Body{T}, ω::Union{T,SVector{1,T}}) where T
    ω = joint.V3' * ω
    Δω = vrotate(ω, joint.qoff) # in body1 frame
    return Δω
end

@inline function getPositionDelta(joint::Rotational2, body1::AbstractBody, body2::Body{T}, θ::Union{T,SVector{1,T}}) where T
    q = Quaternion(cos(θ/2),(joint.V3*sin(θ/2))...)
    Δq = joint.qoff * q # in body1 frame
    return Δq
end

@inline function setForce!(joint::Rotational2, body1::Body, body2::Body{T}, τ::Union{T,SVector{1,T}}, No) where T
    τ1 = vrotate(joint.V3' * -τ, body1.q[No] * joint.qoff)
    τ2 = -τ1

    body1.τ[No] = τ1
    body2.τ[No] = τ2
    return
end

@inline function setForce!(joint::Rotational2, body1::Origin, body2::Body{T}, τ::Union{T,SVector{1,T}}, No) where T
    body2.τ[No] = vrotate(joint.V3' * τ, joint.qoff)
    return
end

@inline function minimalCoordinates(joint::Rotational2, body1::Origin, body2::Body, No)
    q2 = joint.qoff \ body2.q[No]
    joint.V3 * axis(q2) * angle(q2)
end

@inline function minimalCoordinates(joint::Rotational2, body1::Body, body2::Body, No)
    q2 = joint.qoff \ (body1.q[No] \ body2.q[No])
    joint.V3 * axis(q2) * angle(q2)
end

@inline function g(joint::Rotational2, body1::Origin, body2::Body, Δt, No)
    joint.V12 * VLᵀmat(joint.qoff) * getq3(body2, Δt)
end

@inline function g(joint::Rotational2, body1::Body, body2::Body, Δt, No)
    joint.V12 * (VLᵀmat(joint.qoff) * Lᵀmat(getq3(body1, Δt)) * getq3(body2, Δt))
end

@inline function ∂g∂posa(joint::Rotational2{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        X = @SMatrix zeros(T, 2, 3)
        R = -joint.V12 * VLᵀmat(joint.qoff) * Rmat(body2.q[No]) * RᵀVᵀmat(body1.q[No])

        return [X R]
    else
        return ∂g∂posa(joint)
    end
end

@inline function ∂g∂posb(joint::Rotational2{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        X = @SMatrix zeros(T, 2, 3)
        R = joint.V12 * VLᵀmat(joint.qoff) * Lᵀmat(body1.q[No]) * LVᵀmat(body2.q[No])

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂vela(joint::Rotational2{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        V = @SMatrix zeros(T, 2, 3)
        Ω = joint.V12 * VLᵀmat(joint.qoff) * Rmat(ωbar(body2, Δt)) * Rmat(body2.q[No]) * Rᵀmat(body1.q[No]) * Tmat(T) * derivωbar(body1, Δt)

        return [V Ω]
    else
        return ∂g∂vela(joint)
    end
end

@inline function ∂g∂velb(joint::Rotational2{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        V = @SMatrix zeros(T, 2, 3)
        Ω = joint.V12 * VLᵀmat(joint.qoff) * Lᵀmat(ωbar(body1, Δt)) * Lᵀmat(body1.q[No]) * Lmat(body2.q[No]) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end

@inline function ∂g∂posb(joint::Rotational2{T}, body1::Origin, body2::Body, No) where T
    if body2.id == joint.cid
        X = @SMatrix zeros(T, 2, 3)
        R = joint.V12 * VLᵀmat(joint.qoff) * LVᵀmat(body2.q[No])

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂velb(joint::Rotational2{T}, body1::Origin, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        V = @SMatrix zeros(T, 2, 3)
        Ω = joint.V12 * VLᵀmat(joint.qoff) * Lmat(body2.q[No]) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end


text:
function fetch_stypes(xroot)
    stype_vals = findall("//member[@values]", xroot)
    res = OrderedDict()
    for stype ∈ stype_vals
        type = stype.parentnode["name"]
        stype_value = stype["values"]
        res[type] = stype_value
    end
    res
end

stypes = fetch_stypes(xroot)


text:
module JBmad

function track(particles)
  (dummy, NP) = size(particles)
  ccall((:track, "../lib/libtest.so"), Cvoid, (Ref{Cdouble},Ref{Int64}), particles, NP)
end

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ e81954be-8854-46eb-9f26-0e2d5c8896e5
begin
    import Pkg
    Pkg.activate(Base.current_project())
    Pkg.instantiate()
end

# ╔═╡ 8b4426b8-588b-49ad-83c5-a79ed698d704
using PlutoUI

# ╔═╡ 18fbfdc6-9f0b-451f-a6df-fb602725a90c
using DSP

# ╔═╡ 8382586d-170f-45fd-82eb-e1b526b2cf1c
using JLD2

# ╔═╡ f07ab63f-dfb4-493c-a40c-6b3363a858c3
using Plots

# ╔═╡ 58ba75d2-ab31-43b5-a266-df944fd79530
using StatsBase

# ╔═╡ 746daaba-22c5-4a47-822a-3a13e687cd37
md"# Data Preparation"

# ╔═╡ cc3e29fc-b063-4421-a2a2-4fe29cb46abe
@bind dirpath TextField(default="../data/2021-01-12_20-25-30")

# ╔═╡ d07a8c82-7341-4629-b772-b3d8b5324e66
filepaths = filter(path -> endswith(path, ".bin"), readdir(dirpath, join=true))

# ╔═╡ 719cdc2e-72db-4790-b606-12c797dbdf76
function convertdata(rawdata, eltype::Type{T}) where T <: AbstractFloat
	data = (convert(Vector{eltype}, ntoh.(reinterpret(Int16, rawdata))))
	data[1:2:end], data[2:2:end]
end

# ╔═╡ 775f53db-5c7c-481a-b1c7-08c6fc49e08d
function readbins(filepaths, nb, eltype::Type{T}) where T <: AbstractFloat
	re = r".+_ch(?P<first_channel_num>[0-9]+)&(?P<second_channel_num>[0-9]+)\.bin"
	data = Dict{Int, Vector{eltype}}()
	for filepath in filepaths
		m = match(re, filepath)
		rawdata = read(filepath, nb)
		first_channel_data, second_channel_data = convertdata(rawdata, eltype)
		data[parse(Int, m[:first_channel_num])] = first_channel_data
		data[parse(Int, m[:second_channel_num])] = second_channel_data
	end
	sort(data)
end

# ╔═╡ ba1ac3d7-8376-4955-aeee-1d9daf83964e
md"""

## Bandpass

Enable? $(@bind tofilter CheckBox(default=true))

Lowpass frequency (kHz) $(@bind lopassfreq Slider(1:2000, default=50, show_value=true))

Highpass frequency (kHz) $(@bind hipassfreq Slider(1:2000, default=600, show_value=true))

Number of poles in Butterworth filter $(@bind bwpolesnum NumberField(1:10, default=4))
"""

# ╔═╡ 92061d02-57f6-4073-8fea-726a3d82a153
md"""## Resample

Enable? $(@bind toresample CheckBox(default=true))

Decimation factor $(@bind resamplefactor Slider(2:25, default=5, show_value=true))
"""

# ╔═╡ 88625f5f-6e3d-4688-81ae-d6feb3646f8b
samplefreq = 10_000 # in KHz

# ╔═╡ 885e5df3-aa89-4009-88da-23f8e05ca86a
begin
	data = readbins(filepaths, typemax(Int), Float64)
	if (tofilter && hipassfreq > lopassfreq)
		responsetype = Bandpass(lopassfreq, hipassfreq, fs=samplefreq)
		designmethod = Butterworth(bwpolesnum)
		Threads.@threads for n in collect(keys(data))
			ys = data[n]
			xs = axes(data[n], 1)
			β = (mean(xs .* ys) - mean(xs) * mean(ys)) / std(xs, corrected=false)
			α = mean(ys) - β * mean(xs)
			δs = @. ys - α - β * xs
			data[n] = filtfilt(digitalfilter(responsetype, designmethod), δs)
		end
	end
	if toresample
		Threads.@threads for n in collect(keys(data))
			data[n] = resample(data[n], 1 //resamplefactor)
		end
	end
	data
end

# ╔═╡ 67b94e20-6332-4993-bdef-8cc29e51b1ea
md"""
Plot x-range $(@bind npts Slider(100:100:5000, default=1000, show_value=true))

Plot window position $(@bind center_percent Slider(0:0.01:100, default=50, show_value=true))
"""

# ╔═╡ dd071ca8-ff4e-4af2-9d44-48e7947c3431
let
	center = round(Int, (center_percent / 100) * minimum(size(series, 1) for (_, series) in data))
	plots = []
	numchannels = length(data)
	for (nch, series) in data
		a = max(firstindex(series, 1), center - div(npts, 2))
		b = min(lastindex(series, 1), center + div(npts, 2))
		ylim = extrema(view(series, a:b))
		push!(plots, plot(a:b, view(series, a:b),
			  ylim=ylim,
			  title="Channel $nch",
			  titlealign=:left,
			  showaxis=false,
			  ticks=nothing,
			  label=nothing))
	end
	plot(plots..., layout=(numchannels, 1), size=(800, 100numchannels))
end

# ╔═╡ 60930b7e-c200-4219-929b-58e33cd8b939
new_samplefreq = samplefreq // resamplefactor

# ╔═╡ a4728746-a8a0-4bbb-9e9e-228eee132fdf
pgrams = Dict(channel => periodogram(series, fs=new_samplefreq) for (channel, series) in data)

# ╔═╡ 246153f0-c3ed-4cec-b88c-3d8bb5790c17
md"""
Which channel? $(@bind pgramplotch Select(sort(collect(keys(pgrams)))))

Downsampling $(@bind pgram_downsample Slider(10:10:100, default=50, show_value=true))
"""

# ╔═╡ db92c8ac-9fce-4e36-bb2f-de917654af00
plot(resample(pgrams[pgramplotch].freq, 1//pgram_downsample),
	 resample(pgrams[pgramplotch].power, 1//pgram_downsample),
     title="Channel $pgramplotch", label=nothing, dpi=300)

# ╔═╡ 01c58f23-0da5-4d5f-9530-c96d43c655a0
bad_channels = [1, 2, 3]

# ╔═╡ a86b6c8a-9b90-4388-8208-d257ec6951d2
@bind outputpath TextField(default="../data/2021-01-12_20-25-30/data.jld2")

# ╔═╡ 908feaef-7304-48b4-b8a5-bea3ab7d8e43
save(outputpath, "data", filter(!(x -> x.first in bad_channels), data))

# ╔═╡ Cell order:
# ╟─746daaba-22c5-4a47-822a-3a13e687cd37
# ╠═e81954be-8854-46eb-9f26-0e2d5c8896e5
# ╠═8b4426b8-588b-49ad-83c5-a79ed698d704
# ╠═18fbfdc6-9f0b-451f-a6df-fb602725a90c
# ╠═8382586d-170f-45fd-82eb-e1b526b2cf1c
# ╠═f07ab63f-dfb4-493c-a40c-6b3363a858c3
# ╠═cc3e29fc-b063-4421-a2a2-4fe29cb46abe
# ╠═d07a8c82-7341-4629-b772-b3d8b5324e66
# ╠═719cdc2e-72db-4790-b606-12c797dbdf76
# ╠═775f53db-5c7c-481a-b1c7-08c6fc49e08d
# ╟─ba1ac3d7-8376-4955-aeee-1d9daf83964e
# ╠═92061d02-57f6-4073-8fea-726a3d82a153
# ╠═88625f5f-6e3d-4688-81ae-d6feb3646f8b
# ╠═58ba75d2-ab31-43b5-a266-df944fd79530
# ╠═885e5df3-aa89-4009-88da-23f8e05ca86a
# ╟─67b94e20-6332-4993-bdef-8cc29e51b1ea
# ╟─dd071ca8-ff4e-4af2-9d44-48e7947c3431
# ╟─60930b7e-c200-4219-929b-58e33cd8b939
# ╟─a4728746-a8a0-4bbb-9e9e-228eee132fdf
# ╟─246153f0-c3ed-4cec-b88c-3d8bb5790c17
# ╟─db92c8ac-9fce-4e36-bb2f-de917654af00
# ╠═01c58f23-0da5-4d5f-9530-c96d43c655a0
# ╟─a86b6c8a-9b90-4388-8208-d257ec6951d2
# ╠═908feaef-7304-48b4-b8a5-bea3ab7d8e43


text:
using FEASTSolver
using LinearAlgebra
# using MatrixMarket
using SparseArrays

function mmwrite(filename, matrix :: SparseMatrixCSC)
  open(filename, "w") do file
    elem = eltype(matrix) <: Bool ? "pattern" :
           eltype(matrix) <: Integer ?  "integer" :
           eltype(matrix) <: AbstractFloat ? "real" :
           eltype(matrix) <: Complex ? "complex" :
           error("Invalid matrix type")
      sym = "general"

      # write mm header
      write(file, "%%MatrixMarket matrix coordinate $elem $sym\n")

      # only use lower triangular part of symmetric and Hermitian matrices
      #if issymmetric(matrix) || ishermitian(matrix)
      #    matrix = tril(matrix)
      #end

      # write matrix size and number of nonzeros
      write(file, "$(size(matrix, 1)) $(size(matrix, 2)) $(nnz(matrix))\n")

      rows = rowvals(matrix)
      vals = nonzeros(matrix)
      for i in 1:size(matrix, 2)
          for j in nzrange(matrix, i)
              write(file, "$(rows[j]) $i")
              if elem == "pattern" # omit values on pattern matrices
              elseif elem == "complex"
                  write(file, " $(real(vals[j])) $(imag(vals[j]))")
              else
                  write(file, " $(vals[j])")
              end
              write(file, "\n")
          end
      end
  end
end

    N = diagm(-1 => [1,1,1,1,1,1,1])
    Mh0 = 1/6*(4*I + N + N')
    Mh1 = N - N'
    Mh2 = -1*(2*I - N - N')
    Mh3 = Mh1
    Mh4 = -Mh2
    c = [0.6 1.3; 1.3 0.1; 0.1 1.2; 1.0 1.0; 1.2 1.0]
    I8 = Matrix(I, 8, 8)
    M0 = c[1,1] * kron(I8, Mh0) + c[1,2] * kron(Mh0, I8)
    M1 = c[2,1] * kron(I8, Mh1) + c[2,2] * kron(Mh1, I8)
    M2 = c[3,1] * kron(I8, Mh2) + c[3,2] * kron(Mh2, I8)
    M3 = c[4,1] * kron(I8, Mh3) + c[4,2] * kron(Mh3, I8)
    M4 = c[5,1] * kron(I8, Mh4) + c[5,2] * kron(Mh4, I8)
    function butterfly(A0, A1, A2, A3, A4)
        f(z) =  z^4*A4 + z^3*A3 + z^2*A2 + z*A1 + A0
    end

    mmwrite("data/butterflyM0.mtx", sparse(M0))
    mmwrite("data/butterflyM1.mtx", sparse(M1))
    mmwrite("data/butterflyM2.mtx", sparse(M2))
    mmwrite("data/butterflyM3.mtx", sparse(M3))
    mmwrite("data/butterflyM4.mtx", sparse(M4))


text:
using DataFrames
using ProgressMeter
using Sockets

# ============================================================================ #
# [Helpers]                                                                    #
# ============================================================================ #

function write_query(sock::ClickHouseSock, query::AbstractString)::Nothing
    query = ClientQuery("", sock.client_info, "", 2, 0, query)
    write_packet(sock, query)
    write_packet(sock, make_block())
    nothing
end

function dict2columns(
    dict::Dict{Symbol, T} where T,
    valid_columns::Dict{Symbol, String},
)::Vector{Column}
    @assert begin
        diff = symdiff(dict |> keys |> Set, valid_columns |> keys |> Set)
        isempty(diff)
    end "Mismatched columns: $(diff)"

    # TODO: Check if column types match.

    [
        Column(string(name), valid_columns[name], column)
        for (name, column) ∈ dict
    ]
end

function columns2dict(cols::Vector{Column})::Dict{Symbol, Any}
    Dict(Symbol(x.name) => x.data for x ∈ cols)
end

function make_block(columns::Vector{Column} = Column[])::Block
    num_columns = length(columns)
    num_rows = num_columns == 0 ? 0 : length(columns[1].data)
    @assert all(length(col.data) == num_rows for col ∈ columns)
    Block("", BlockInfo(), num_columns, num_rows, columns)
end

"""
Check dirty state & mark as dirty.

This is intententionally not a `guard(x) do y .. end` function because the
closure would result in all captured variables to be boxed, resulting in
runtime overhead. It is also intentional that the dirty flag isn't cleared
in case of exceptions between `enter_dirty` and `exit_dirty`.
"""
function enter_dirty(sock::ClickHouseSock)::Nothing
    !sock.dirty || error(
        "The socket is dirty. This means that another query is either " *
        "currently in progress or an error left the socket in a dirty state. " *
        "Please create a new connection."
    )

    sock.dirty = true
    nothing
end

function exit_dirty(sock::ClickHouseSock)::Nothing
    @assert sock.dirty
    sock.dirty = false
    nothing
end

# ============================================================================ #
# [Queries]                                                                    #
# ============================================================================ #

"Establish a connection to a given ClickHouse instance."
function connect(
    host::AbstractString = "localhost",
    port::Integer = 9000;
    database::AbstractString = "",
    username::AbstractString = "default",
    password::AbstractString = "",
)::ClickHouseSock
    sock = ClickHouseSock(Sockets.connect(host, port))

    # Say hello to the server!
    write_packet(sock, ClickHouse.ClientHello(
        CLIENT_NAME,
        DBMS_VER_MAJOR,
        DBMS_VER_MINOR,
        DBMS_VER_REV,
        database,
        username,
        password,
    ))

    # Read server info.
    server_info = read_server_packet(sock)::ServerInfo
    sock.server_tz = server_info.server_timezone

    sock
end

"Send a ping request and wait for the response."
function ping(sock::ClickHouseSock)::Nothing
    enter_dirty(sock)
    write_packet(sock, ClientPing())
    read_server_packet(sock)::ServerPong
    exit_dirty(sock)
    nothing
end

"Execute a DDL query."
function execute(
    sock::ClickHouseSock,
    ddl_query::AbstractString,
)::Nothing
    enter_dirty(sock)
    write_query(sock, ddl_query)
    read_server_packet(sock)::ServerEndOfStream
    exit_dirty(sock)
    nothing
end

"""
Insert blocks into a table, reading from an iterable.
The iterable is expected to yield values of type `Dict{Symbol, Any}`.
"""
function insert(
    sock::ClickHouseSock,
    table::AbstractString,
    iter,
)::Nothing where T
    enter_dirty(sock)

    # TODO: We might want to escape the table name here...
    write_query(sock, "INSERT INTO $(table) VALUES")

    packet = read_server_packet(sock)
    sample_block = if packet isa ServerTableColumns
        packet.sample_block
    else
        packet::Block
    end

    valid_columns = Dict(
        Symbol(x.name) => x.type
        for x ∈ sample_block.columns
    )

    for block_dict ∈ iter
        columns = dict2columns(block_dict, valid_columns)
        block = make_block(columns)
        write_packet(sock, block)
    end

    # Empty block = end of data.
    write_packet(sock, make_block())

    read_server_packet(sock)::ServerEndOfStream

    exit_dirty(sock)
    nothing
end

"Execute a query, invoking a callback for each block."
function select_callback(
    callback::Function,
    sock::ClickHouseSock,
    query::AbstractString;
    show_progress::Bool = false,
    progress_kwargs::Dict = Dict(),
)::Nothing
    enter_dirty(sock)
    write_query(sock, query)

    sample_block = read_server_packet(sock)
    sample_block::Block
    @assert UInt64(sample_block.num_rows) == 0

    progress_bar = nothing
    progress_rows::Int64 = 0
    while true
        packet = read_server_packet(sock)

        handle(x::ServerProfileInfo) = true
        handle(x::ServerEndOfStream) = false

        function handle(x::ServerProgress)
            if show_progress
                if progress_bar === nothing
                    progress_bar = Progress(
                        x.total_rows |> UInt64 |> Int64;
                        progress_kwargs...
                    )
                end
                rows = x.rows |> UInt64 |> Int64
                progress_rows += rows
                update!(progress_bar, progress_rows)
            end

            true
        end

        function handle(block::Block)
            if UInt64(block.num_rows) != 0
                block.columns |> columns2dict |> callback
            end
            true
        end

        if !handle(packet)
            break
        end
    end

    exit_dirty(sock)
end

"Execute a query, streaming the resulting blocks through a channel."
function select_df_channel(
    query::AbstractString,
    sock::ClickHouseSock = connect();
    csize = 0,
    kwargs...,
)::Channel{DataFrame}
    Channel(ctype = DataFrame, csize = csize) do ch
        select_callback(sock, query; kwargs...) do row
            put!(ch, DataFrame!(row))
        end
    end
end

function select_channel(
    query::AbstractString,
    sock::ClickHouseSock = connect();
    csize = 0,
    kwargs...,
)::Channel{Dict{Symbol, Any}}
    Channel(ctype = Dict{Symbol, Any}, csize = csize) do ch
        select_callback(sock, query; kwargs...) do row
            put!(ch, row)
        end
    end
end

"Execute a query, flattening blocks into a single dict of column arrays."
function select(
    query::AbstractString,
    sock::ClickHouseSock = connect();
    kwargs...
)::Dict{Symbol, Any}
    result = Dict{Symbol, Any}()
    select_callback(sock, query; kwargs...) do block
        for (col_name, col_data) ∈ block
            arr = get(result, col_name, nothing)
            if arr === nothing
                result[col_name] = col_data
            else
                append!(arr, col_data)
            end
        end
    end
    result
end

"Execute a query, flattening blocks into a dataframe."
function select_df(
    query::AbstractString,
    sock::ClickHouseSock = connect();
    kwargs...
)::DataFrame
    DataFrame(select(query, sock; kwargs...))
end

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


text:
#sim_res = run_sim(FT(22),FT(1.5),FT(5),FT(10),FT(0.025*0.4),FT(32));

#sim_res1 = run_sim(FT(1.8e1), FT(8.2e-1), FT(1.5e1), FT(9.5e-4), FT(2.5), FT(2.0), FT(1.1e3), FT(4.2e-1), 48*365*1 - 52*48, 48*365, 0.35);

using StatsBase
using DataFrames
using CSV
using Parameters
using Land
using Land.Photosynthesis
using Land.CanopyLayers
using Land.PlantHydraulics
using Land.SoilPlantAirContinuum
using Land.StomataModels
using Thermodynamics

using CLIMAParameters
using CLIMAParameters:AbstractEarthParameterSet
using CLIMAParameters.Planet: grav

const FT = Float32

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

KG_H_2_MOL_S = FT(55.55 / 3600);
mpa2mm = FT(10^6/9.8);
K_STEFAN = FT(Stefan());

#deltaT = FT(60*60)

include(string(PROJECT_HOME,"/simulation_code/create_spac_setsoil.jl"))
#include(string(PROJECT_HOME,"/simulation_code/create_spac_multi.jl"))

include(string(PROJECT_HOME,"/simulation_code/land_utils4.jl"))
include(string(PROJECT_HOME,"/simulation_code/cap_funs_inv.jl"))
include(string(PROJECT_HOME,"/simulation_code/create_node_newvol_varyPV.jl"))

function run_sim_vary(vcmax_par, k_rel_crit, k_weibB, k_weibC, k_plant, k_soil, z_soil, istart, N, smc0, storage_mult, buffrate,scheme_number,df_raw,g1,deltaT,alpha,n)

df = deepcopy(df_raw[istart:(istart+N-1),:]);

df[!,"leafpot"] = zeros(N)
df[!,"leafpotStore"] = zeros(N)

df[!,"VPD"] = zeros(N)

#df[!,"leafvol"] = zeros(N)
df[!,"stempot"] = zeros(N)
df[!,"soilpot"] = zeros(N)
df[!,"ETmod"] = zeros(N)
df[!,"trunkflow"] = zeros(N)

#df[!,"APARsun"] = zeros(N)
#df[!,"Zenith"] = zeros(N)

#df[!,"APARshade"] = zeros(N)

#df[!,"APAR_ps"] = zeros(N)

#df[!,"Ecrit"] = zeros(N)

df[!,"leafStore"] = zeros(N)

df[!,"ColumnSMC"] = zeros(N)
df[!,"ColumnSWP"] = zeros(N)
#df[!,"fvis"] = zeros(N)
df[!,"glw"] = zeros(N)
#df[!,"GPP"] = zeros(N)
df[!,"NPP"] = zeros(N)
#df[!,"Pcrit"] = zeros(N)
#df[!,"Runoff"] = zeros(N)

#df[!,"pl1"] = zeros(N)
#df[!,"pl2"] = zeros(N)

#if scheme_number == 3
node = create_moflux_node(vcmax_par, k_plant, z_soil, smc0, storage_mult,1,deltaT,alpha,n);
#else
#	node = create_spac(OSMWang{FT}(),vcmax_par,k_plant,z_soil, FT(40),alpha,n);
#end

beta_curve = WeibullSingle(FT(k_weibB*(1-k_rel_crit)), FT(k_weibC));
sm1 = Land.StomataModels.ESMMedlyn(FT(0),FT(g1));
#sm1 = Land.StomataModels.ESMBallBerry(FT(0),FT(g1));

node.stomata_model = sm1;

#k_weibB = 3;
#k_weibC = 2;
#k_rel_crit = FT(0.33);
update_Weibull!(node, FT(k_weibB), FT(k_weibC));

@unpack angles, can_opt, can_rad, canopy_rt, envirs, f_SL, ga, in_rad,
		latitude, leaves_rt, n_canopy, photo_set, plant_hs, plant_ps,
		rt_con, soil_opt, stomata_model, wl_set, n_root = node;
@unpack lidf, nAzi, nIncl = canopy_rt;
@unpack dWL, iPAR = wl_set;
in_rad_bak = deepcopy(in_rad);
nSL = nAzi * nIncl;
in_Erad = in_rad_bak.E_direct .+ in_rad_bak.E_diffuse;
#in_PPFD = sum( e2phot(dWL, in_Erad)[iPAR] ) * FT(1e6);
in_PPFD = sum( Land.CanopyLayers.e2phot(wl_set.WL, in_Erad/1000)[iPAR] .* dWL[iPAR] ) * FT(1e6);

maxvol = get_v_max(plant_hs);

if scheme_number == 2

	#node.plant_hs.branch[1].v_maximum = node.plant_hs.branch[2].v_maximum
	#node.plant_hs.branch[1].v_storage = node.plant_hs.branch[2].v_storage

	mypv = PVCurveLinear{Float32}(0.2f0, FT(buffrate));
	for i_can in 1:n_canopy
		plant_hs.leaves[i_can].pv = mypv
		plant_hs.branch[i_can].pv = mypv
	end

	plant_hs.trunk.pv = mypv;

	for i_root in 1:n_root
		plant_hs.roots[i_root].pv = mypv
	end

	kstore = storage_mult
	for i in 1:n_canopy
		plant_hs.leaves[i].v_maximum = FT(kstore) * plant_hs.leaves[i].v_maximum
		plant_hs.branch[i].v_maximum = FT(kstore) * plant_hs.branch[i].v_maximum
	end

	plant_hs.trunk.v_maximum = FT(kstore) * plant_hs.trunk.v_maximum

	for i in 1:n_root
		plant_hs.roots[i].v_maximum = FT(kstore) * plant_hs.roots[i].v_maximum
	end

	kinit = 0.95
	for i in 1:n_canopy
		plant_hs.leaves[i].v_storage = FT(kinit) * plant_hs.leaves[i].v_maximum
		plant_hs.branch[i].v_storage = FT(kinit) * plant_hs.branch[i].v_maximum
	end

	plant_hs.trunk.v_storage = FT(kinit) * plant_hs.trunk.v_maximum

	for i in 1:n_root
		plant_hs.roots[i].v_storage = FT(kinit) * plant_hs.roots[i].v_maximum
	end

	w_soil = FT(smc0);
	node.swc = ones(FT, length(node.swc))*w_soil;
	# update soil water matrices
	# need to adjust SWC to avoid problem in residual SWC at Niwot Ridge
	ψ_soil = soil_p_25_swc(plant_hs.roots[1].sh, w_soil);
	for root in plant_hs.roots
		root.p_ups = ψ_soil;
	end;

	pot_init = 1+0.2*ψ_soil

	for i in 1:n_canopy
		plant_hs.leaves[i].v_storage = FT(pot_init) * plant_hs.leaves[i].v_maximum
		plant_hs.branch[i].v_storage = FT(pot_init) * plant_hs.branch[i].v_maximum
	end

	plant_hs.trunk.v_storage = FT(pot_init) * plant_hs.trunk.v_maximum

	for i in 1:n_root
		plant_hs.roots[i].v_storage = FT(pot_init) * plant_hs.roots[i].v_maximum
	end

end

if scheme_number == 1
	w_soil = FT(smc0);
	node.swc = ones(FT, length(node.swc))*w_soil;
	# update soil water matrices
	# need to adjust SWC to avoid problem in residual SWC at Niwot Ridge
	ψ_soil = soil_p_25_swc(plant_hs.roots[1].sh, w_soil);
	for root in plant_hs.roots
		root.p_ups = ψ_soil;
	end;
end

smc_record = zeros(FT,N,length(node.swc));
psi_record_leaf = zeros(FT,N,length(node.plant_hs.leaves));
psi_record_branch = zeros(FT,N,length(node.plant_hs.leaves));
qin_record = zeros(FT,N,length(node.plant_hs.roots));
psi_record_trunk = zeros(FT,N);

rbase =  Q10TD{FT}(4, 298.15, 1.7)

# iterate through the weather data
for i in eachindex(df.Day)

#=
	if mod(i, (24*365)) == 1
		println(Int((i-1)/(24*365)))
		#println(node.plant_hs.leaves[length(node.plant_hs.leaves)].p_leaf)
	end
=#

	t_soil = FT(df.T_SOIL[1] + 273.15);
	update_LAI!(node, FT(df.LAI_modis[i]));
	#update_LAI!(node, FT(df.LAI_modis[i]*0.5 + 1));
	#update_LAI!(node, FT(mean(df.LAI_modis*0.5 .+ 1)));

	# update PAR related information
	zenith = zenith_angle(latitude, FT(df.Day[i]), FT(df.Hour[i]),
						  FT(df.Minu[i]));
	zenith = min(88, zenith);
	angles.sza = zenith;
	in_rad.E_direct  .= in_rad_bak.E_direct  .* df.PPFD_in[i] ./ in_PPFD;
	in_rad.E_diffuse .= in_rad_bak.E_diffuse .* df.PPFD_in[i] ./ in_PPFD;
	canopy_geometry!(canopy_rt, angles, can_opt, rt_con);
	canopy_matrices!(leaves_rt, can_opt);
	short_wave!(canopy_rt, can_opt, can_rad, in_rad, soil_opt, rt_con);
	canopy_fluxes!(canopy_rt, can_opt, can_rad, in_rad, soil_opt,
				   leaves_rt, wl_set, rt_con);

	#update_Weibull!(node, FT(weib_b), FT(3.0));
	#update_Kmax!(node, FT(k_plant));

	#df.pl1[i] = node.plant_hs.leaves[1].p_leaf*1

	# update fluxes
	f_H₂O = 0;
	f_CO₂ = 0;
	for i_can in 1:n_canopy
		iEN = envirs[i_can];
		iHS = plant_hs.leaves[i_can];
		iPS = plant_ps[i_can];
		iRT = n_canopy + 1 - i_can;

		# update environmental conditions
		iEN.t_air = df.T_AIR[i] + 273.15;
		iEN.p_atm = df.P_ATM[i] * 1000;
		iEN.p_a   = iEN.p_atm * 4e-4;
		iEN.p_O₂  = iEN.p_atm * 0.209;
		iEN.p_sat = saturation_vapor_pressure( param_set, iEN.t_air , Liquid());

		iEN.RH    = FT(df.RelHum[i]);
		iEN.p_H₂O = iEN.RH * iEN.p_sat;
		iEN.vpd   = (FT(1)-iEN.RH) * iEN.p_sat;

		#iEN.vpd   = df.VPD[i] * 100;
		#iEN.p_H₂O = iEN.p_sat - iEN.vpd;
		#iEN.RH    = iEN.p_H₂O / iEN.p_sat;

		iEN.wind  = df.WIND[i];

		# prescribe leaf temperature
		_tl = (df.LW_OUT[i] / 0.97 / K_STEFAN ) ^ 0.25;
		iPS.T = _tl;
		update_leaf_TP!(photo_set, iPS, iHS, iEN);
		temperature_effects!(iHS, FT(_tl));

		# calculate the fraction of sunlit and shaded leaves
		f_view = (can_opt.Ps[iRT] + can_opt.Ps[iRT+1]) / 2;
		for iLF in 1:nSL
			iPS.APAR[iLF] = can_rad.absPAR_sun[(iRT-1)*nSL+iLF] *
							FT(1e6);
			iPS.LAIx[iLF] = f_view * f_SL[iLF];
		end;
		iPS.APAR[end] = can_rad.absPAR_shade[iRT] * FT(1e6);
		iPS.LAIx[end] = 1 - f_view;

	end

	#df.pl2[i] = node.plant_hs.leaves[1].p_leaf*1

	glw_mean = 0

	subIter = 15

	for i_can in 1:n_canopy

		iEN = envirs[i_can];
		iHS = plant_hs.leaves[i_can];
		iPS = plant_ps[i_can];
		iRT = n_canopy + 1 - i_can;

		iHS.p_crt = FT(-iHS.vc.b * log(FT(1/k_rel_crit))^(1/iHS.vc.c));
		iPS.ec    = critical_flow(iHS, iPS.ec);
        iPS.ec    = max(FT(0), iPS.ec);

		# iterate for 15 times to find steady state solution
		for subI in 1:subIter
		# calculate the photosynthetic rates
			gas_exchange_new!(photo_set, iPS, iEN, GswDrive());
			#prognostic_gsw!(iPS, iEN, stomata_model, FT(1), deltaT/subIter);
			update_gsw!(iPS, stomata_model, photo_set, iEN, FT(deltaT/subIter),FT(1));
			#update_gsw!(iPS, stomata_model, photo_set, iEN, FT(deltaT/subIter));
			gsw_control!(photo_set, iPS, iEN);
		end

		iHS.p_crt = -iHS.vc.b * log(FT(1000)) ^ (1/iHS.vc.c);
		iPS.ec    = critical_flow(iHS, iPS.ec);
        iPS.ec    = max(FT(0), iPS.ec);

		# update the flow rates
		for iLF in 1:(nSL+1)
			f_CO₂ += iPS.An[iLF] * iPS.LAIx[iLF] * iPS.LA;
			f_H₂O += iPS.g_lw[iLF] * (iPS.p_sat - iEN.p_H₂O) / iEN.p_atm *
					 iPS.LAIx[iLF] * iPS.LA;
		end;

		glw_mean += sum(iPS.g_lw)/(nSL+1)/n_canopy;
	end;

	# update flow profile and pressure history along the tree
	if scheme_number==1
		branch_total_flow = 0;
		for i_can in 1:n_canopy
			iEN = envirs[i_can];
			iLF = plant_hs.leaves[i_can];
			iPS = plant_ps[i_can];
			iST = plant_hs.branch[i_can];
			iLF.flow = sum(iPS.g_lw .* iPS.LAIx) *
					   (iPS.p_sat - iEN.p_H₂O) / iEN.p_atm;
			#iLF.flow = FT(max(0,df.LE[i])/44200)*node.ga/node.la
			iST.flow = iLF.flow * iPS.LA;
			branch_total_flow += iST.flow;
		end;
		plant_hs.trunk.flow = branch_total_flow; #sum([iST.flow for iST in plant_hs.branch]);
		for iRT in plant_hs.roots
			iRT.flow = plant_hs.trunk.flow / length(plant_hs.roots);
		end;
	end

	if (scheme_number==2) | (scheme_number==3)
		for i_can in 1:n_canopy
			iEN = envirs[i_can];
			iLF = plant_hs.leaves[i_can];
			iPS = plant_ps[i_can];
			iST = plant_hs.branch[i_can];
			iLF.q_out = sum(iPS.g_lw .* iPS.LAIx) *
					   (iPS.p_sat - iEN.p_H₂O) / iEN.p_atm #/ iPS.LAI ;
			#iLF.q_out = FT(max(0,df.LE[i])/44200)*node.ga/node.la;

		end;
	end

	subIter2 = 8

	if scheme_number==3
		dd1, dd2 = update_cap_mat!(plant_hs,deltaT);
		#=
		try
			dd1, dd2 = update_cap_mat!(plant_hs,deltaT);
		catch err
			return df, smc_record, psi_record_leaf, psi_record_branch, psi_record_trunk, node
		end
		=#
		#=
		try
			dd1, dd2 = find_new_cap(plant_hs,deltaT);
			if mean(dd1 .<= maxvol) == 1
				push_cap!(plant_hs,dd1,dd2)
			else
				cap_iter = 100
				for icap in 1:cap_iter
					dd1, dd2 = find_new_cap(plant_hs,deltaT/cap_iter);
					push_cap!(plant_hs,dd1,dd2)
				end
			end
		catch err
			return df, smc_record, psi_record_leaf, psi_record_branch, psi_record_trunk, node
			#=
			cap_iter = 4
			for icap in 1:cap_iter
				dd1, dd2 = find_new_cap(plant_hs,deltaT/cap_iter);
				push_cap!(plant_hs,dd1,dd2)
			end
			=#
		end
		=#
	end

	for subI2 in 1:subIter2

		if scheme_number==1
			pressure_profile!(node.plant_hs, SteadyStateMode(); update=false);
			do_soil_drain!(node, FT(df.RAIN[i])/subIter2, deltaT/subIter2, k_soil,FT(0))
		end
		if scheme_number==2
			do_soil_nss_drain!(node, FT(df.RAIN[i])/subIter2, deltaT/subIter2, k_soil,FT(0))
			update_PVF!(node.plant_hs,deltaT/subIter2);
		end
		if scheme_number==3
			do_soil_nss_drain!(node, FT(df.RAIN[i])/subIter2, deltaT/subIter2, k_soil, FT(0.0))
		end

		for i_root in eachindex(node.plant_hs.roots)
			rootI = node.plant_hs.root_index_in_soil[i_root]
			node.plant_hs.roots[i_root].p_ups = soil_p_25_swc(node.plant_hs.roots[1].sh, node.swc[rootI])
		end

		# update canopy layer p_ups, which will be passed to each leaf
		for _i_can in eachindex(node.plant_hs.leaves)
			_iHS = node.plant_hs.leaves[_i_can];
			_iPS = node.plant_ps[_i_can];
			_iPS.p_ups = _iHS.p_ups;
		end
	end

	# calculate respiration from other components
	_r = photo_TD_from_set(rbase, t_soil);

	for i_can in 1:n_canopy
		ratioI = xylem_k_ratio(beta_curve, plant_hs.leaves[i_can].p_leaf);
		update_VJR_leaf!(node.plant_ps[i_can], ratioI)
	end

	# update the data frame
	df.glw[i] = glw_mean
	df.NPP[i] = f_CO₂ / ga - _r;
	df.ETmod[i] = f_H₂O / ga;
	df.trunkflow[i] = plant_hs.trunk.q_in
	df.stempot[i] = mean(node.plant_hs.trunk.p_element);
	df.leafpot[i] = node.plant_hs.leaves[length(node.plant_hs.leaves)].p_leaf;
	smc_record[i,:] = node.swc
    df.ColumnSMC[i] = sum(node.swc .* -diff(node.soil_bounds)) / abs(node.soil_bounds[end])
    swpI = [soil_p_25_swc(node.plant_hs.roots[1].sh, x) for x in node.swc];
    df.ColumnSWP[i] = sum(swpI .* -diff(node.soil_bounds)) / abs(node.soil_bounds[end])
	#df.Ecrit[i] = node.plant_ps[length(node.plant_hs.leaves)].ec;
	#df.Pcrit[i] = node.plant_hs.leaves[length(node.plant_hs.leaves)].p_crt;
	df.leafStore[i] = node.plant_hs.leaves[length(node.plant_hs.leaves)].v_storage;
	df.leafpotStore[i] = node.plant_hs.leaves[length(node.plant_hs.leaves)].p_storage;
	df.VPD[i] = node.envirs[1].vpd

	for ican in eachindex(node.plant_hs.leaves)
		psi_record_leaf[i,ican] = node.plant_hs.leaves[ican].p_leaf;
		psi_record_branch[i,ican] = mean(node.plant_hs.branch[ican].p_element);
	end
	psi_record_trunk[i] = mean(node.plant_hs.trunk.p_element);
	qin_record[i,:] = [x.q_in for x in node.plant_hs.roots];

end;

return df, smc_record, psi_record_leaf, psi_record_branch, psi_record_trunk, qin_record#, node

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

using NodeNumbering

adjacency = Dict(
            1 => [3, 8, 9],
            2 => [3, 8, 7, 6],
            3 => [1, 2],
            4 => [8, 9],
            5 => [7, 8],
            6 => [2, 7],
            7 => [5, 2, 6],
            8 => [1, 2, 4, 5],
            9 => [1, 4]);

finalorder =  Dict(
              7 => 2,
              9 => 9,
              4 => 8,
              2 => 3,
              3 => 5,
              5 => 4,
              8 => 6,
              6 => 1,
              1 => 7);

@testset "Create RCM adjacency" begin
    expected =  Dict(
                7 => [5, 6, 9],
                9 => [7, 8],
                4 => [2, 6],
                2 => [4, 3, 1],
                3 => [5, 6, 2, 1],
                8 => [6, 9],
                5 => [7, 3],
                6 => [7, 3, 8, 4],
                1 => [3, 2])
    result = NodeNumbering.create_RCM_adjacency(adjacency, finalorder)
    @test result == expected
end


text:
using MakieLayout
using Makie

scene = Scene(resolution = (600, 600));
screen = display(scene)
campixel!(scene);

la1 = LayoutedAxis(scene)
la2 = LayoutedAxis(scene)
la3 = LayoutedAxis(scene)
la4 = LayoutedAxis(scene)
la5 = LayoutedAxis(scene)

linkxaxes!(la3, la4)
linkyaxes!(la3, la5)

# lines!(la1.scene, rand(200, 2) .* 100, color=:black, show_axis=false)
img = rand(100, 100)
image!(la1.scene, img, show_axis=false)
lines!(la2.scene, rand(200, 2) .* 100, color=:blue, show_axis=false)

linkeddata = rand(200, 2) .* 100
scatter!(la3.scene, linkeddata, markersize=3, color=:red, show_axis=false)
scatter!(la4.scene, linkeddata, markersize=3, color=:orange, show_axis=false)
scatter!(la5.scene, linkeddata, markersize=3, color=:pink, show_axis=false)
update!(scene)

gl = GridLayout([], 2, 2, [1, 1], [1, 1], 0.01, 0.01)

gl[2, 1:2] = AxisLayout(BBox(75, 0, 0, 75), la1)
gl[1, 2] = AxisLayout(BBox(75, 0, 0, 75), la2)

gl2 = GridLayout([], 2, 2, [0.8, 0.2], [0.2, 0.8], 0.01, 0.01)
gl2[2, 1] = AxisLayout(BBox(75, 0, 0, 75), la3)
gl2[1, 1] = AxisLayout(BBox(75, 0, 0, 75), la4)
gl2[2, 2] = AxisLayout(BBox(75, 0, 0, 75), la5)

gl[1, 1] = gl2

sg = outersolve(gl, BBox(shrinkbymargin(pixelarea(scene)[], 30)))
applylayout(sg)
# when the scene is resized, apply the outersolve'd outermost grid layout
# this recursively updates all layout objects that are contained in the grid
on(scene.events.window_area) do area
    sg = outersolve(gl, BBox(shrinkbymargin(pixelarea(scene)[], 30)))
    applylayout(sg)
end


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

using Test
include("StudentData.jl")
using .StudentData

function insert!(heap, elm;
    min = true,
    compareFn = min ? (x, y) -> x < y : (x, y) -> x > y)
    push!(heap, elm)
    val = percolateup!(heap, min = min, compareFn = compareFn)
    return val
end

function delete!(heap;
    min = true,
    compareFn = min ? (x, y) -> x < y : (x, y) -> x > y)
    isempty(heap) && return
    swap!(heap, 1, lastindex(heap))
    elem = pop!(heap)
    percoloatedown!(heap, min = min, compareFn = compareFn)
    return elem
end

function percolateup!(heap;
    curr = lastindex(heap),
    min = true,
    compareFn = min ? (x, y) -> x < y : (x, y) -> x > y)
    parent = curr ÷ 2
    if parent != 0 && compareFn(heap[curr], heap[parent])
        swap!(heap, curr, parent)
        percolateup!(heap, curr = parent, min = min, compareFn = compareFn)
    else
        heap
    end
end

function findtargetchild(heap, lchild;
    min = true,
    compareFn = min ? (x, y) -> x < y : (x, y) -> x > y)
    rchild = lchild + 1
    if rchild > lastindex(heap)
        lchild
    else
        compareFn(heap[lchild], heap[rchild]) ? lchild : rchild
    end
end

function percoloatedown!(heap;
    curr = 1,
    min = true,
    compareFn = min ? (x, y) -> x < y : (x, y) -> x > y)
    lchild = curr * 2
    lchild > lastindex(heap) && return heap
    targetchild = findtargetchild(heap, lchild, min = min, compareFn = compareFn)
    if compareFn(heap[targetchild], heap[curr])
        swap!(heap, curr, targetchild)
        percoloatedown!(heap, curr = targetchild, min = min, compareFn = compareFn)
    else
        return heap
    end
end

function swap!(heap, index1, index2)
    temp = heap[index1]
    heap[index1] = heap[index2]
    heap[index2] = temp
end

function heapsort(arr::Vector;
    desc = true,
    min = !desc,
    compareFn = min ? (x, y) -> x < y : (x, y) -> x > y)
    insertFn! = (_heap, _elm) -> insert!(_heap, _elm, min = min, compareFn = compareFn)
    deleteFn! = (_heap) -> delete!(_heap, min = min, compareFn = compareFn)
    heap = reduce(insertFn!, arr, init = [])
    res = []
    while !isempty(heap)
        push!(res, deleteFn!(heap))
    end
    return res
end

# tests for heapsort
@testset "Heapsort tests" begin
    @test heapsort([11, 9, 7, 10, 8, 6]) == [11, 10, 9, 8, 7, 6]
    @test heapsort([11, 9, 7, 10, 8, 6], desc = false) == [6, 7, 8, 9, 10, 11]
    randArray = rand(Int8, 10)
    @test heapsort(randArray) == sort(randArray, rev = true) # Note: our heapsort implementation by default prefers descending order
    @test heapsort(randArray, desc = false) == sort(randArray)
    @test heapsort(students, compareFn = studentCompareFn) == sortedStudents
end

max_stars_repo_name: codekomali/JuliaPractice

max_issues_repo_name: codekomali/JuliaPractice

max_forks_repo_name: codekomali/JuliaPractice


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

using Markdown
using InteractiveUtils

# ╔═╡ 0da325a8-d6f6-4963-8562-9b7847878c07
begin
	using Pkg
	Pkg.add("CairoMakie")
	Pkg.add("LaTeXStrings")
	# Pkg.develop(url="/Users/lyx/SynologyDrive/Develop/MakiePublication")
	Pkg.add(url="https://github.com/liuyxpp/MakiePublication.jl")
	using CairoMakie
	using LaTeXStrings
end

# ╔═╡ 9e2e4de7-bf62-4bf9-b263-cb4781a5d46c
using Revise

# ╔═╡ 44af131a-0904-4d36-b098-01823d6fb8a8
using MakiePublication

# ╔═╡ 30b0601e-63b9-4986-a311-9685560b55ef
function myline()
	x = 0:0.001:1.0
	fig = Figure(figure_padding=(0, 6, 0, 4))
	ax = Axis(fig, xlabel=L"x", ylabel=L"f(x)")
	lines!(ax, x, x, label=L"x")
	lines!(ax, x, x.^2, label=L"x^2")
	lines!(ax, x, x.^3, label=L"x^3")
	lines!(ax, x, x.^4, label=L"x^4")
	lines!(ax, x, log.(x.+1), label=L"\ln(1+x)")
	lines!(ax, x, .√x, label=L"x^{1/2}")
	lines!(ax, x, x.^(1/3), label=L"x^{1/3}")
	xlims!(0, 1)
	ylims!(0, 1)

	fig[1,1] = ax

	return fig
end

# ╔═╡ 5b508e70-8921-407f-8e19-9cd2a3c22e80
myline()

# ╔═╡ 57ece2db-0068-4d82-8391-ddffefbde28b
with_theme(myline, theme_web())

# ╔═╡ 2d15b715-6441-4b7c-b751-599c0f5894a3
fig = with_theme(myline, theme_acs())

# ╔═╡ b78fde28-81ed-4c31-a161-43ceef0d7c42
savefig("test_line.pdf", fig)

# ╔═╡ 6c1907e1-5298-438d-b6ef-9adb92d849f1
function myscatter()
	x = 0:0.001:1.0
	xs = 0:0.05:1.0
	fig = Figure(figure_padding=(0, 6, 0, 4))
	ax = Axis(fig, xlabel=L"x", ylabel=L"f(x)")
	l1 = lines!(ax, x, x)
	s1 = scatter!(ax, xs, xs, label=L"x")
	l2 = lines!(ax, x, x.^2)
	s2 = scatter!(ax, xs, xs.^2, label=L"x^2")
	l3 = lines!(ax, x, x.^3)
	s3 = scatter!(ax, xs, xs.^3, label=L"x^3")
	l4 = lines!(ax, x, x.^4)
	s4 = scatter!(ax, xs, xs.^4, label=L"x^4")
	l5 = lines!(ax, x, log.(x.+1))
	s5 = scatter!(ax, xs, log.(xs.+1), label=L"\ln(1+x)")
	l6 = lines!(ax, x, .√x)
	s6 = scatter!(ax, xs, .√xs, label=L"x^{1/2}")
	l7 = lines!(ax, x, x.^(1/3))
	s7 = scatter!(ax, xs, xs.^(1/3), label=L"x^{1/3}")
	xlims!(0, 1)
	ylims!(0, 1)

	axislegend(ax, [[l1, s1], [l2, s2], [l3, s3], [l4, s4], [l5, s5], [l6, s6], [l7, s7]], [L"x", L"x^2", L"x^3", L"x^4", L"\ln(1+x)", L"x^{1/2}", L"x^{1/3}"], nbanks=2, position=:lt, margin=(4,0,0,-5), padding=(0,0,0,0))

	fig[1,1] = ax

	return fig
end

# ╔═╡ 628fac9a-f4ba-4da8-bbe8-7ebb8dc06164
myscatter()

# ╔═╡ 4b6b7880-a391-4661-911a-bb3e144824ef
fig_scatter = with_theme(myscatter, theme_acs())

# ╔═╡ 29faf34e-e1a1-443a-8cda-56ec25a910f8
savefig("test_scatter.pdf", fig_scatter)

# ╔═╡ 372077d2-e2d6-49db-b504-776afea5cef9
fig_scatter2 = with_theme(myscatter, theme_acs(cycle=Cycle([:color], covary=true)))

# ╔═╡ 1a2bf02e-170b-478c-8c33-8529a7832e2e
savefig("test_scatter2.pdf", fig_scatter2)

# ╔═╡ 7e89621b-f8c5-4726-91d7-c32b6ee568b3
fig_web = with_theme(myscatter, theme_web())

# ╔═╡ 51cbde74-93f1-4501-8ee3-6e6e5913782f
save("scatter_web.png", fig_web, px_per_unit=4)

# ╔═╡ Cell order:
# ╠═0da325a8-d6f6-4963-8562-9b7847878c07
# ╠═9e2e4de7-bf62-4bf9-b263-cb4781a5d46c
# ╠═44af131a-0904-4d36-b098-01823d6fb8a8
# ╠═30b0601e-63b9-4986-a311-9685560b55ef
# ╠═5b508e70-8921-407f-8e19-9cd2a3c22e80
# ╠═57ece2db-0068-4d82-8391-ddffefbde28b
# ╠═2d15b715-6441-4b7c-b751-599c0f5894a3
# ╠═b78fde28-81ed-4c31-a161-43ceef0d7c42
# ╠═6c1907e1-5298-438d-b6ef-9adb92d849f1
# ╠═628fac9a-f4ba-4da8-bbe8-7ebb8dc06164
# ╠═4b6b7880-a391-4661-911a-bb3e144824ef
# ╠═29faf34e-e1a1-443a-8cda-56ec25a910f8
# ╠═372077d2-e2d6-49db-b504-776afea5cef9
# ╠═1a2bf02e-170b-478c-8c33-8529a7832e2e
# ╠═7e89621b-f8c5-4726-91d7-c32b6ee568b3
# ╠═51cbde74-93f1-4501-8ee3-6e6e5913782f


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

@doc raw"""CSIDriverSpec is the specification of a CSIDriver.

    IoK8sApiStorageV1beta1CSIDriverSpec(;
        attachRequired=nothing,
        podInfoOnMount=nothing,
        volumeLifecycleModes=nothing,
    )

    - attachRequired::Bool : attachRequired indicates this CSI volume driver requires an attach operation (because it implements the CSI ControllerPublishVolume() method), and that the Kubernetes attach detach controller should call the attach volume interface which checks the volumeattachment status and waits until the volume is attached before proceeding to mounting. The CSI external-attacher coordinates with CSI volume driver and updates the volumeattachment status when the attach operation is complete. If the CSIDriverRegistry feature gate is enabled and the value is specified to false, the attach operation will be skipped. Otherwise the attach operation will be called.
    - podInfoOnMount::Bool : If set to true, podInfoOnMount indicates this CSI volume driver requires additional pod information (like podName, podUID, etc.) during mount operations. If set to false, pod information will not be passed on mount. Default is false. The CSI driver specifies podInfoOnMount as part of driver deployment. If true, Kubelet will pass pod information as VolumeContext in the CSI NodePublishVolume() calls. The CSI driver is responsible for parsing and validating the information passed in as VolumeContext. The following VolumeConext will be passed if podInfoOnMount is set to true. This list might grow, but the prefix will be used. \"csi.storage.k8s.io/pod.name\": pod.Name \"csi.storage.k8s.io/pod.namespace\": pod.Namespace \"csi.storage.k8s.io/pod.uid\": string(pod.UID) \"csi.storage.k8s.io/ephemeral\": \"true\" iff the volume is an ephemeral inline volume                                 defined by a CSIVolumeSource, otherwise \"false\"  \"csi.storage.k8s.io/ephemeral\" is a new feature in Kubernetes 1.16. It is only required for drivers which support both the \"Persistent\" and \"Ephemeral\" VolumeLifecycleMode. Other drivers can leave pod info disabled and/or ignore this field. As Kubernetes 1.15 doesn't support this field, drivers can only support one mode when deployed on such a cluster and the deployment determines which mode that is, for example via a command line parameter of the driver.
    - volumeLifecycleModes::Vector{String} : VolumeLifecycleModes defines what kind of volumes this CSI volume driver supports. The default if the list is empty is \"Persistent\", which is the usage defined by the CSI specification and implemented in Kubernetes via the usual PV/PVC mechanism. The other mode is \"Ephemeral\". In this mode, volumes are defined inline inside the pod spec with CSIVolumeSource and their lifecycle is tied to the lifecycle of that pod. A driver has to be aware of this because it is only going to get a NodePublishVolume call for such a volume. For more information about implementing this mode, see https://kubernetes-csi.github.io/docs/ephemeral-local-volumes.html A driver can support one or more of these modes and more modes may be added in the future.
"""
mutable struct IoK8sApiStorageV1beta1CSIDriverSpec <: SwaggerModel
    attachRequired::Any # spec type: Union{ Nothing, Bool } # spec name: attachRequired
    podInfoOnMount::Any # spec type: Union{ Nothing, Bool } # spec name: podInfoOnMount
    volumeLifecycleModes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: volumeLifecycleModes

    function IoK8sApiStorageV1beta1CSIDriverSpec(;attachRequired=nothing, podInfoOnMount=nothing, volumeLifecycleModes=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1beta1CSIDriverSpec, Symbol("attachRequired"), attachRequired)
        setfield!(o, Symbol("attachRequired"), attachRequired)
        validate_property(IoK8sApiStorageV1beta1CSIDriverSpec, Symbol("podInfoOnMount"), podInfoOnMount)
        setfield!(o, Symbol("podInfoOnMount"), podInfoOnMount)
        validate_property(IoK8sApiStorageV1beta1CSIDriverSpec, Symbol("volumeLifecycleModes"), volumeLifecycleModes)
        setfield!(o, Symbol("volumeLifecycleModes"), volumeLifecycleModes)
        o
    end
end # type IoK8sApiStorageV1beta1CSIDriverSpec

const _property_map_IoK8sApiStorageV1beta1CSIDriverSpec = Dict{Symbol,Symbol}(Symbol("attachRequired")=>Symbol("attachRequired"), Symbol("podInfoOnMount")=>Symbol("podInfoOnMount"), Symbol("volumeLifecycleModes")=>Symbol("volumeLifecycleModes"))
const _property_types_IoK8sApiStorageV1beta1CSIDriverSpec = Dict{Symbol,String}(Symbol("attachRequired")=>"Bool", Symbol("podInfoOnMount")=>"Bool", Symbol("volumeLifecycleModes")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }) = collect(keys(_property_map_IoK8sApiStorageV1beta1CSIDriverSpec))
Swagger.property_type(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1beta1CSIDriverSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1beta1CSIDriverSpec[property_name]

function check_required(o::IoK8sApiStorageV1beta1CSIDriverSpec)
    true
end

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


text:
@testset "Contour" begin
    @polyvar x y
    p1 = 1 + x + y + x * y + y^2 + x^2 * y

    C = PolynomialAmoebas.contour(p1, res=(100, 100))
    @test length(C.ps) != 0
end


text:
include("../src/driver.jl")
using PyPlot
using Printf
function plot_lyapunov_exponents()
	npts = 10
	w = rand(npts)
	les = rand(npts)
	η = [0.001,0.01,0.02,0.025,0.05]
	fig, ax = subplots()
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)

	for (i, ηi) = enumerate(η)
		w .= rand(npts)
		@show i, ηi
		for j = 1:npts
			les[j] = lyap_exp(w[j], ηi)
		end
		ax.plot(w, les, label=latexstring(L"\eta = ",@sprintf("%0.3f", ηi)))
	end
end
function plot_bifurcation()
	npts = 50
	wstar = rand(npts)
	nSteps = 50000
	η = LinRange(0.0001, 0.1, 20)
	fig, ax = subplots()
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)

	for (i, ηi) = enumerate(η)
		@show i, ηi
		wstar .= rand(npts)
		for j = 1:npts
			wstar[j] = run(wstar[j], ηi, nSteps)
		end
		ax.plot(ηi*ones(npts), wstar, "bo", ms=3.0)
	end
	ax.grid(true)
	ax.set_xlabel(L"$ \eta $", fontsize=30)
	ax.set_ylabel(L"$ w^* $", fontsize=30)
end

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET