text:
const origin_options = ["sw", "se", "nw", "ne", "center", "centroid"]

"""
    transform_rotate([geojson::T[, angle::Real], pivot::Point=nothing, mutate::Bool=false) where {T <: AbstractGeometry}

Rotates any geojson Feature or Geometry of a specified angle, around its `centroid` or a given `pivot` point;
all rotations follow the right-hand rule.

# Examples
```jldoctest
julia> point = Point([-75.69926351308823,45.43145021122502])
Point([-75.6993, 45.4315])

julia> transform_rotate(geojson=point, angle=80., pivot=Point([-75.6, 45.3]))
Point([-75.433, 45.3915])
```
"""
function transform_rotate(; geojson::T, angle::Real, pivot::Point=nothing, mutate::Bool=false) where {T <: AbstractGeometry}
    if angle === 0
        return geojson
    end

    geo = []

    if mutate
        geo = deepcopy(geojson)
    else
        geo = geojson
    end

    type = geotype(geo)

    pivot === nothing && (pivot = centroid(geo))
    coords = geo.coordinates

    if type === :Point

        initAngle = rhumb_bearing(pivot.coordinates, coords)
        finalAngle = initAngle + angle
        dist = rhumb_distance(pivot.coordinates, coords)
        newCoords = rhumb_destination(pivot.coordinates, dist, finalAngle).coordinates
        coords[1] = newCoords[1]
        coords[2] = newCoords[2]

    elseif type === :Polygon || type === :MultiLineString

        for i in eachindex(coords)
            for j in eachindex(coords[i])
                initAngle = rhumb_bearing(pivot.coordinates, coords[i][j])
                finalAngle = initAngle + angle
                dist = rhumb_distance(pivot.coordinates, coords[i][j])
                newCoords = rhumb_destination(pivot.coordinates, dist, finalAngle).coordinates
                coords[i][j][1] = newCoords[1]
                coords[i][j][2] = newCoords[2]
            end
        end
    elseif type === :LineString

        for i in eachindex(coords)
            initAngle = rhumb_bearing(pivot.coordinates, coords[i])
            finalAngle = initAngle + angle
            dist = rhumb_distance(pivot.coordinates, coords[i])
            newCoords = rhumb_destination(pivot.coordinates, dist, finalAngle).coordinates
            coords[i][1] = newCoords[1]
            coords[i][2] = newCoords[2]
        end
    end

    return geo
end

transform_rotate!(; geojson::T, angle::Real, pivot::Point=nothing) where {T <: AbstractGeometry} = transform_rotate(geojson=geojson, angle=angle, pivot=pivot, mutate=true)

"""
    transform_translate([geojson::T[, distance::R[, direction::R]]], vertical::R=0, mutate::Bool=false, units::String="kilometers") where {T <: Union{AbstractFeature, AbstractGeometry}, R <: Real}

Moves any geojson Feature or Geometry of a specified distance along a Rhumb Line
on the provided direction angle.

# Examples
```jldoctest
julia> poly = Polygon([[[0, 29], [3.5, 29], [2.5, 32], [0, 29]]])
Polygon(Array{Array{Float64,1},1}[[[0.0, 29.0], [3.5, 29.0], [2.5, 32.0], [0.0, 29.0]]])

julia> transform_translate(poly, 300, 70)
Polygon(Array{Array{Float64,1},1}[[[2.91184, 29.9228], [6.41184, 29.9228], [5.50479, 32.9228], [2.91184, 29.9228]]])
```
"""
function transform_translate(geojson::T, distance::R, direction::R, vertical::R=0, mutate::Bool=false, units::String="kilometers") where {T <: Union{AbstractFeature, AbstractGeometry}, R <: Real}
    (distance == 0 && vertical == 0) && return geojson

    if distance < 0
        distance  = -distance
        direction = -direction
    end

    !mutate && (geojson = deepcopy(geojson))

    coords = []
    type = nothing

    if typeof(geojson) <: AbstractGeometry
        coords = geojson.coordinates
        type = geotype(geojson)
    else
        coords = geojson.geometry.coordinates
        type = geotype(geojson.geometry)
    end

    if type === :Point
        newCoords = rhumb_destination(coords, distance, direction, units).coordinates
        coords[1] = newCoords[1]
        coords[2] = newCoords[2]
        (vertical != 0 && length(coords) === 3) && (coords[3] += vertical)

    elseif type === :Polygon || type === :MultiLineString
        for i in eachindex(coords)
            for j in eachindex(coords[i])
                newCoords = rhumb_destination(coords[i][j], distance, direction, units).coordinates
                coords[i][j][1] = newCoords[1]
                coords[i][j][2] = newCoords[2]
                (vertical != 0 && length(coords[i][j]) === 3) && (coords[i][j][3] += vertical)
            end
        end
    elseif type === :LineString
        for i in eachindex(coords)
            newCoords = rhumb_destination(coords[i], distance, direction, units).coordinates
            coords[i][1] = newCoords[1]
            coords[i][2] = newCoords[2]
            (vertical != 0 && length(coords[i]) === 3) && (coords[i][3] += vertical)
        end
    end

    return geojson
end

transform_translate!(geojson::T, distance::R, direction::R, vertical::R=0, units::String="kilometers") where {T <: Union{AbstractFeature, AbstractGeometry}, R <: Real} = transform_translate(geojson, distance, direction, vertical, true, units)

"""
    transform_scale([geojson::T[, factor::Float64]], origin::String="centroid") where {T <: AbstractFeatureCollection}

Scale a GeoJson from a given point by a factor of scaling (ex: factor=2 would make the GeoJson 200% larger).
If a FeatureCollection is provided, the origin point will be calculated based on each individual Feature.

# Examples
```jldoctest
julia> coll = FeatureCollection([Feature(Point([-75.69926351308823,45.43145021122502])), Feature(Polygon([[[0, 29], [3.5, 29], [2.5, 32], [0, 29]]]))])
FeatureCollection{Feature}(Feature[Feature(Point([-75.6993, 45.4315]), Dict{String,Any}()), Feature(Polygon(Array{Array{Float64,1},1}[[[0.0, 29.0], [3.5, 29.0], [2.5, 32.0], [0.0, 29.0]]]), Dict{String,Any}())], nothing, nothing)

julia> transform_scale(coll, 0.1)
FeatureCollection{Feature}(Feature[Feature(Point([-75.6993, 45.4315]), Dict{String,Any}()), Feature(Polygon(Array{Array{Float64,1},1}[[[1.3495, 29.675], [1.70067, 29.675], [1.59896, 29.975], [1.3495, 29.675]]]), Dict{String,Any}())], nothing, nothing)
```
"""
function transform_scale(geojson::T, factor::Real, origin::String="centroid", mutate::Bool=false) where {T <: AbstractFeatureCollection}

    geo = []
    if !mutate
        geo = deepcopy(geojson)
    else
        geo = geojson
    end

    for (i, feat) in enumerate(geo.features)
        geo.features[i] = scale(feat, factor, origin)
    end

    return geo
end

transform_scale!(geojson::T, factor::Real, origin::String="centroid") where {T <: AbstractFeatureCollection} = transform_scale(geojson, factor, origin, true)

"""
    scale([feature::Feature[, factor::Real]], origin::String="centroid")

Scale a Feature.

# Examples
```jldoctest
julia> feature = Feature(Polygon([[[0, 29], [3.5, 29], [2.5, 32], [0, 29]]]))
Feature(Polygon(Array{Array{Float64,1},1}[[[0.0, 29.0], [3.5, 29.0], [2.5, 32.0], [0.0, 29.0]]]), Dict{String,Any}())

julia> scale(feature, 0.1)
Feature(Polygon(Array{Array{Float64,1},1}[[[1.3495, 29.675], [1.70067, 29.675], [1.59896, 29.975], [1.3495, 29.675]]]), Dict{String,Any}())
```
"""
function scale(feature::Feature, factor::Real, origin::String="centroid")

    if !(origin in origin_options)
        throw(error("'$(origin)' is not a valid option. Allowed values are: $(origin_options)"))
    end

    (geotype(feature.geometry) === :Point || factor === 1) && return feature

    type = geotype(feature.geometry)
    coords = feature.geometry.coordinates
    center = getOrigin(feature, origin)

    if type === :LineString

        for i in eachindex(coords)
            start = rhumb_distance(center.coordinates, coords[i])
            bearing = rhumb_bearing(center.coordinates, coords[i])
            distance = start * factor

            newCoords = rhumb_destination(center.coordinates, distance, bearing).coordinates

            coords[i][1] = newCoords[1]
            coords[i][2] = newCoords[2]

            if length(coords) === 3
                coords[i][3] *= factor
            end
        end
    elseif type === :Polygon || type === :MultiLineString

        for i in eachindex(coords)
            for j in eachindex(coords[i])
                start = rhumb_distance(center.coordinates, coords[i][j])
                bearing = rhumb_bearing(center.coordinates, coords[i][j])
                distance = start * factor

                newCoords = rhumb_destination(center.coordinates, distance, bearing).coordinates

                coords[i][j][1] = newCoords[1]
                coords[i][j][2] = newCoords[2]

                if length(coords) === 3
                    coords[i][j][3] *= factor
                end
            end
        end
    end

    return feature

end

function getOrigin(geojson::Feature, origin::String)
    if origin == nothing || origin === ""
        origin = "centroid"
    end

    sw = ["sw", "southwest", "westsouth", "bottomleft"]
    se = ["se", "southeast", "eastsouth", "bottomright"]
    nw = ["nw", "northwest", "westnorth", "topleft"]
    ne = ["ne", "northeast", "eastnorth", "topright"]

    box = bbox(geojson)
    west = box[1]
    south = box[2]
    east = box[3]
    north = box[4]

    if origin in sw
        return Point([west, south])
    elseif origin in se
        return Point([east, south])
    elseif origin in nw
        return Point([west, north])
    elseif origin in ne
        return Point([east, north])
    elseif origin === "center"
        return center(geojson.geometry)
    elseif origin === "centroid"
        return centroid(geojson.geometry)
    else
        throw(error("Invalid origin."))
    end
end

"""
    explode([geojson::T], pointsOnly::Bool=false)::FeatureCollection where {T <: Union{AbstractFeatureCollection, AbstractGeometry}}

Takes a Geometry or a FeatureCollection and returns all positions as Points.

# Examples
```jldoctest
julia> poly = Polygon([[[100, 0], [101, 0], [101, 1], [100, 1], [100, 0]]])
Polygon(Array{Array{Float64,1},1}[[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]]])

julia> explode(poly, true)
5-element Array{Point,1}:
 Point([100.0, 0.0])
 Point([101.0, 0.0])
 Point([101.0, 1.0])
 Point([100.0, 1.0])
 Point([100.0, 0.0])

```
"""
function explode(geojson::T, pointsOnly::Bool=false) where {T <: Union{AbstractFeatureCollection, AbstractGeometry}}
    points::Vector{Point} = []

    if geotype(geojson) === :FeatureCollection
        for feat in geojson.features
            geom = feat.geometry

            if geotype(geom) === :Point
                push!(points, geom)

            elseif geotype(geom) === :Polygon || geotype(geom) === :MultiLineString

                for i in eachindex(geom.coordinates)
                    for j in eachindex(geom.coordinates[i])
                        push!(points, Point(geom.coordinates[i][j]))
                    end
                end
            elseif geotype(geom) === :LineString

                for i in eachindex(geom.coordinates)
                    push!(points, Point(geom.coordinates[i]))
                end

            end
        end
    else
        if geotype(geojson) === :Point
            push!(points, geojson)

        elseif geotype(geojson) === :Polygon || geotype(geojson) === :MultiLineString

            for i in eachindex(geojson.coordinates)
                for j in eachindex(geojson.coordinates[i])
                    push!(points, Point(geojson.coordinates[i][j]))
                end
            end
        elseif geotype(geojson) === :LineString

            for i in eachindex(geojson.coordinates)
                push!(points, Point(geojson.coordinates[i]))
            end

        end

    end

    pointsOnly && return points
    return FeatureCollection([Feature(x) for x in points])
end

"""
    flip([geojson::T], mutate::Bool=false) where {T <: Union{AbstractFeature, AbstractGeometry}}

Take input Features and Geometries and flips all of their coordinates from `[x, y]` to `[y, x]`.

# Examples
```jldoctest
julia> point = Point([77.34374999999999,43.58039085560784,3000])
Point([77.3437, 43.5804, 3000.0])

julia> flip(point)
Point([43.5804, 77.3437, 3000.0])
```
"""
function flip(geojson::T, mutate::Bool=false) where {T <: Union{AbstractFeature, AbstractGeometry}}
    type = geotype(geojson)
    coords = []
    geom = nothing

    geo = []
    if !mutate
        geo = deepcopy(geojson)
    else
        geo = geojson
    end

    if type === :Feature
        geom = geo.geom
        coords = geo.geometry.coordinates
    else
        geom = geo
        coords = geo.coordinates
    end

    if geotype(geom) === :Point
        x = coords[1]
        y = coords[2]
        coords[1] = y
        coords[2] = x
    elseif geotype(geom) === :Polygon || geotype(geom) === :MultiLineString
        for i in eachindex(coords)
            x = coords[i][1]
            y = coords[i][2]
            coords[i][1] = y
            coords[i][2] = x
        end
    elseif geotype(geom) === :LineString
        for i in eachindex(coords[1])
            x = coords[1][i][1]
            y = coords[1][i][2]
            coords[1][i][1] = y
            coords[1][i][2] = x
        end
    end

    return geo

end

flip!(geojson::T) where {T <: Union{AbstractFeature, AbstractGeometry}} = flip(geojson, true)

"""
Cohen-Sutherland line clippign algorithm, adapted to efficiently to
handle polylines rather than just segments
"""
function lineclip(points::Vector{P}, bbox::Vector{T}) where {T <: Real, P <: AbstractPoint}
    len = length(points)
    codeA = bitcode(points[1], bbox)
    parts = []
    results = []
    let a, b, codeB, lastCode end

    for i in 2:len
        a = points[i - 1]
        b = points[i]

        codeB = lastCode = bitcode(b, bbox)

        while true
            if !(codeA | codeB)
                push!(parts, a)

                if codeB !== lastCode
                    push!(parts, b)

                    if i < len - 1
                        push!(results, parts)
                        parts = []
                    end
                elseif i === len - 1
                    push!(parts, b)
                end

                break
            elseif (codeA & codeB)
                break

            elseif codeA
                a = intersection(a, b, codeB, bbox)
                codeA = bitcode(b, bbox)

            else
                b = intersection(a, b, codeB, bbox)
                codeB = bitcode(b, bbox)
            end
        end
        codeA = lastCode
    end
    length(part) > 0 && push!(results, parts)

    return results
end

"""
Sutherland-Hodgeman polygon clipping algorithm.
"""
function polygonclip(points::Vector{P}, bbox::Vector{T}) where {P <: AbstractPoint, T <: Real}
    let result, edge, prev, pInside, p, inside end
    ranges = [1, 2, 4, 8]

    for edge in ranges
        result = []

        prev = points[length(points) - 1]
        pInside = !(bitcode(prev, bbox) & edge)

        for i in 1:length(points)
            p = points[i]
            inside = !(bitcode(p, bbox) & edge)

            inside !== pInside && push!(result, intersection(prev, p, edge, bbox))

            inside && push!(result, p)

            prev = p
            pInside = inside
        end

        points = result
        length(points) <= 0 && break
    end

    return result
end

"""
Intersect a segment against one of the 4 lines that make up the bbox
"""
function intersection(a::Point, b::Point, edge, bbox::Vector{T}) where {T <: Real}
    return edge & 8 ? [a[1] + (b[1] - a[1]) * (bbox[4] - a[1]) / (b[1] - a[1]), bbox[4]] : # top
           edge & 4 ? [a[1] + (b[1] - a[1]) * (bbox[2] - a[2]) / (b[2] - a[2]), bbox[2]] : # bottom
           edge & 2 ? [bbox[3], a[2] + (b[2] - a[2]) * (bbox[3] - a[1]) / (b[1] - a[1])] : # right
           edge & 1 ? [bbox[1], a[2] + (b[2] - a[2]) * (bbox[1] - a[1]) / (b[1] - a[1])] : # left
           nothing
end

@inline function bitcode(p::Point, bbox::Vector{T}) where {T <: Real}
    code = 0

    p[1] < bbox[1] && (code |= 1)
    p[1] > bbox[3] && (code |= 2)

    p[2] < bbox[2] && (code |= 4)
    p[2] > bbox[4] && (code |= 8)

    return code
end

"""
    polygon_tangents(pt::Point, poly::Polygon)

Finds the tangents of a Polygon from a Point.

# Examples
```jldoctest
julia> point = Point([92.46093749999999,54.67383096593114])
Point([92.4609, 54.6738])

julia> poly = Polygon([[[100, 0], [101, 0], [101, 1], [100, 1], [100, 0]]])
Polygon(Array{Array{Float64,1},1}[[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]]])

julia> poly = Polygon([[[48.1641, 20.6328], [76.6406, 20.6328], [76.6406, 38.8226], [48.1641, 38.8226], [48.1641, 20.6328]]])
Polygon(Array{Array{Float64,1},1}[[[48.1641, 20.6328], [76.6406, 20.6328], [76.6406, 38.8226], [48.1641, 38.8226], [48.1641, 20.6328]]])

julia> polygon_tangents(point, poly)
FeatureCollection{Feature}(Feature[Feature(Point([48.1641, 38.8226]), Dict{String,Any}()), Feature(Point([76.6406, 20.6328]), Dict{String,Any}())], nothing, nothing)
 ```
"""
function polygon_tangents(pt::Point, poly::Polygon)
    ptCoords = pt.coordinates
    polyCoords = poly.coordinates

    box = bbox(poly)
    nearest = nothing
    nearestIndex = 1

    if ptCoords[1] > box[1] && ptCoords[1] < box[3] && ptCoords[2] > box[2] && ptCoords[2] < box[4]
        nearest = nearestpoint(pt, explode(poly, true))
        nearestIndex = findfirst(x -> x == nearest, explode(poly, true))
    end

    rtan = polyCoords[1][nearestIndex]
    ltan = polyCoords[1][1]

    if nearest != nothing
        nearest.coordinates[2] < ptCoords[2] && (ltan = polyCoords[1][nearestIndex])
    end

    enext = nothing
    eprev = isleft(polyCoords[1][1], polyCoords[1][length(polyCoords[1]) - 1], ptCoords)
    out = process(polyCoords[1], ptCoords, eprev, enext, rtan, ltan)

    rtan = out[1]
    ltan = out[2]

    return FeatureCollection([Feature(Point(rtan)), Feature(Point(ltan))])
end

function isleft(point1, point2, point3)
    return (point2[1] - point1[1]) * (point3[2] - point1[2]) - (point3[1] - point1[1]) * (point2[2] - point1[2])
end

function isabove(point1, point2, point3)
    return isleft(point1, point2, point3) > 0
end

function isbelow(point1, point2, point3)
    return isleft(point1, point2, point3) < 0
end

function process(polyCoords, ptCoords, eprev, enext, rtan, ltan)
    for i in 1:length(polyCoords) - 1
        curr = polyCoords[i]
        next = polyCoords[i + 1]

        (i === length(polyCoords) - 1) && (next = polyCoords[1])

        enext = isleft(curr, next, ptCoords)
        if eprev <= 0 && enext > 0
            !(isbelow(ptCoords, curr, rtan)) && (rtan = curr)
        elseif eprev > 0 && enext <= 0
            !(isabove(ptCoords, curr, ltan)) && (ltan = curr)
        end
        eprev = enext
    end
    return [rtan, ltan]
end

"""
    polygon_to_line([poly::Polygon])

Converts a Polygon to LineString or MultiLineString

# Examples
```jldoctest
julia> poly = Polygon([[[-2.275543, 53.464547],[-2.275543, 53.489271],[-2.215118, 53.489271],[-2.215118, 53.464547],[-2.275543, 53.464547]]])
Polygon(Array{Array{Float64,1},1}[[[-2.27554, 53.4645], [-2.27554, 53.4893], [-2.21512, 53.4893], [-2.21512, 53.4645], [-2.27554, 53.4645]]])

julia> polygon_to_line(poly)
LineString(Array{Float64,1}[[-2.27554, 53.4645], [-2.27554, 53.4893], [-2.21512, 53.4893], [-2.21512, 53.4645], [-2.27554, 53.4645]])
```
"""
function polygon_to_line(poly::Polygon)
    return coordinatesToLine(poly.coordinates)
end

function coordinatesToLine(coords::Vector{Vector{Position}})
    length(coords) > 1 && return MultiLineString(coords)
    return LineString(coords[1])
end

"""
    convert_to([geojson::AbstractGeometry[, projection::String]], mutate::Bool=false)

Convert a GeoJSON geojsonect to the defined `projection`.
"""
function convert_to(geojson::AbstractGeometry, projection::String, mutate::Bool=false)
    allowed_proj = ["mercator", "wgs84"]

    !(projection in allowed_proj) && throw(error("$(projection) is not a valid option."))

    type = geotype(geojson)

    type === :Point && return Point((projection === "mercator") ? to_mercator(geojson) : to_WGS84(geojson))

    geo = []
    if !mutate
        geo = deepcopy(geojson)
    else
        geo = geojson
    end

    coords = geo.coordinates
    if type === :LineString
        for i in eachindex(coords)
            newCoords = (projection === "mercator") ? to_mercator(Point(coords[i])) : to_WGS84(Point(coords[i]))
            coords[i][1] = newCoords[1]
            coords[i][2] = newCoords[2]
        end
    elseif type === :Polygon
        for i in eachindex(coords[1])
            newCoords = (projection === "mercator") ? to_mercator(Point(coords[1][i])) : to_WGS84(Point(coords[1][i]))
            coords[1][i][1] = newCoords[1]
            coords[1][i][2] = newCoords[2]
        end
    else
        throw(error("$(type) is not supported."))
    end

    return geo
end

convert_to(geojson::Feature, projection::String, mutate::Bool=false) = convert_to(geojson.geometry, projection, mutate)
convert_to!(geojson::Feature, projection::String) = convert_to(geojson.geometry, projection, true)
convert_to!(geojson::AbstractGeometry, projection::String) = convert_to(geojson, projection, true)

"""
    combine([ft::FeatureCollection])

Combine a FeatureCollection of Point, LineString, or Polygon features
into MultiPoint, MultiLineString, or MultiPolygon features.

# Examples
```jldoctest
julia> l1 = LineString([[102.0,-10.0],[130.0,4.0]])
LineString(Array{Float64,1}[[102.0, -10.0], [130.0, 4.0]])

julia> l2 = LineString([[40.0,-20.0],[150.0,18.0]])
LineString(Array{Float64,1}[[40.0, -20.0], [150.0, 18.0]])

julia> combine(FeatureCollection([Feature(l1), Feature(l2)]))
FeatureCollection{Feature}(Feature[Feature(MultiLineString(Array{Array{Float64,1},1}[[[102.0, -10.0], [130.0, 4.0], [40.0, -20.0], [150.0, 18.0]]]), Dict{String,Any}())], nothing, nothing)
```
"""
function combine(ft::FeatureCollection)

    groups = Dict{String, Any}([
        "MultiPoint" => [[], Dict{String, Any}()],
        "MultiLineString" => [[], Dict{String, Any}()],
        "MultiPolygon" => [[], Dict{String, Any}()]]
    )

    function add_to_group(ft::Feature, key::String, multi::Bool)
        if multi
            if geotype(ft.geometry) === :MultiPoint
                for i in eachindex(ft.geometry.coordinates)
                    push!(groups[key][1], ft.geometry.coordinates[i])
                end
            elseif geotype(ft.geometry) === :MultiLineString
                for i in eachindex(ft.geometry.coordinates[1])
                    push!(groups[key][1], ft.geometry.coordinates[1][i])
                end
            else
                for i in eachindex(ft.geometry.coordinates[1][1])
                    push!(groups[key][1], ft.geometry.coordinates[1][1][i])
                end
            end
        else
            if geotype(ft.geometry) === :Point
                groups[key][1] = vcat(groups[key][1], [ft.geometry.coordinates])
            elseif geotype(ft.geometry) === :LineString
                groups[key][1] = vcat(groups[key][1], ft.geometry.coordinates)
            elseif geotype(ft.geometry) === :Polygon
                groups[key][1] = vcat(groups[key][1], ft.geometry.coordinates...)
            end
        end # function

        merge!(groups[key][2], ft.properties)
    end

    for feat in ft.features
        type = geotype(feat.geometry)
        try
            add_to_group(feat, string(type), true)
        catch e
            if isa(e, KeyError)
                newType = "Multi" * string(type)
                add_to_group(feat, newType, false)
            else
                throw(e)
            end
        end
    end

    features::Vector{Feature} = []

    groups = filter( x -> length(groups[x.first][1]) > 0, groups)

    for g in keys(groups)
        if g === "MultiPoint"
            point = MultiPoint(values(groups[g][1]))
            push!(features, Feature(point, groups[g][2]))
        elseif g === "MultiLineString"
            point = MultiLineString([values(groups[g][1])])
            push!(features, Feature(point, groups[g][2]))
        elseif g === "MultiPolygon"
            point = MultiPolygon([[values(groups[g][1])]])
            push!(features, Feature(point, groups[g][2]))
        end
    end

    return FeatureCollection(features)
end

"""
    tag(fc1::FeatureCollection, fc2::FeatureCollection, in_field::String, out_field::String)::FeatureCollection

Take a set of Points and a set of Polygons and performs a spatial join.
"""
function tag(fc1::FeatureCollection, fc2::FeatureCollection, in_field::String, out_field::String)::FeatureCollection
    points = deepcopy(fc1)
    polys = deepcopy(fc2)

    for feat in points.features
        geotype(feat.geometry) !== :Point && throw(error("$(feat.geometry) is not a supported geometry."))
        for feat2 in polys.features
            geotype(feat2.geometry) !== :Polygon && throw(error("$(feat2.geometry) is not a supported geometry."))

            if point_in_polygon(feat.geometry, feat2.geometry)
                feat.properties[out_field] = feat2.properties[in_field]
            end
        end
    end
    return points
end

"""
    simplify(geojson::FeatureCollection, tolerance::Real=1., hq::Bool=false, mutate::Bool=false)

Take a FeatureCollection and return a simplified version.
Internally uses an adaptation of [simplify-js](http://mourner.github.io/simplify-js/) to perform simplification using the Ramer-Douglas-Peucker algorithm.

# Examples
```jldoctest
julia> poly = Polygon([[
    [-70.603637, -33.399918],
    [-70.614624, -33.395332],
    [-70.639343, -33.392466],
    [-70.659942, -33.394759],
    [-70.683975, -33.404504],
    [-70.697021, -33.419406],
    [-70.701141, -33.434306],
    [-70.700454, -33.446339],
    [-70.694274, -33.458369],
    [-70.682601, -33.465816],
    [-70.668869, -33.472117],
    [-70.646209, -33.473835],
    [-70.624923, -33.472117],
    [-70.609817, -33.468107],
    [-70.595397, -33.458369],
    [-70.587158, -33.442901],
    [-70.587158, -33.426283],
    [-70.590591, -33.414248],
    [-70.594711, -33.406224],
    [-70.603637, -33.399918]]])

julia> simplify(poly, 0.01)
Polygon(Array{Array{Float64,1},1}[[[-70.6036, -33.3999], [-70.684, -33.4045], [-70.7011, -33.4343], [-70.6943, -33.4584], [-70.6689, -33.4721], [-70.6098, -33.4681], [-70.5872, -33.4429], [-70.6036, -33.3999]]])
```
"""
function simplify(geojson::FeatureCollection, tolerance::Real=1., hq::Bool=false, mutate::Bool=false)
    tolerance < 0 && throw(error("Invalid tolerance. It must be a positive number!"))

    !mutate && (geojson = deepcopy(geojson))

    for feat in geojson.features
        geometry = feat.geometry

        simplify_geometry(geometry, tolerance, hq)
    end

    return geojson
end

"""
    simplify(geojson::FeatureCollection, tolerance::Real=1., hq::Bool=false)

Take a FeatureCollection and return a simplified version, modifying the original FeatureCollection.
Internally uses an adaptation of [simplify-js](http://mourner.github.io/simplify-js/) to perform simplification using the Ramer-Douglas-Peucker algorithm.
"""
simplify!(geojson::FeatureCollection, tolerance::Real=1., hq::Bool=false) = simplify(geojson, tolerance, hq, true)

"""
    simplify(geojson::AbstractGeometry, tolerance::Real=1., hq::Bool=false, mutate::Bool=false)

Take a GeoJSON Geometry and return a simplified version.
Internally uses an adaptation of [simplify-js](http://mourner.github.io/simplify-js/) to perform simplification using the Ramer-Douglas-Peucker algorithm.
"""
function simplify(geojson::AbstractGeometry, tolerance::Real=1., hq::Bool=false, mutate::Bool=false)
    !mutate && (geojson = deepcopy(geojson))
    return simplify_geometry(geojson, tolerance, hq)
end

"""
    simplify(geojson::AbstractGeometry, tolerance::Real=1., hq::Bool=false)

Take a GeoJSON Geometry and return a simplified version, modifying the original Geometry.
Internally uses an adaptation of [simplify-js](http://mourner.github.io/simplify-js/) to perform simplification using the Ramer-Douglas-Peucker algorithm.
"""
simplify!(geojson::AbstractGeometry, tolerance::Real=1., hq::Bool=false) = simplify(geojson, tolerance, hq, true)

"""
    simplify(geojson::AbstractFeature, tolerance::Real=1., hq::Bool=false, mutate::Bool=false)

Take a GeoJSON Feature and return a simplified version.
Internally uses an adaptation of [simplify-js](http://mourner.github.io/simplify-js/) to perform simplification using the Ramer-Douglas-Peucker algorithm.
"""
simplify(geojson::AbstractFeature, tolerance::Real=1., hq::Bool=false, mutate::Bool=false) = simplify(geojson.geometry, tolerance, hq, mutate)

"""
    simplify(geojson::AbstractFeature, tolerance::Real=1., hq::Bool=false)

Take a GeoJSON Feature and return a simplified version, modifying the original Feature.
Internally uses an adaptation of [simplify-js](http://mourner.github.io/simplify-js/) to perform simplification using the Ramer-Douglas-Peucker algorithm.
"""
simplify!(geojson::AbstractFeature, tolerance::Real=1., hq::Bool=false) = simplify(geojson, tolerance, hq, true)

function simplify_geometry(geom::AbstractGeometry, tol::Real, hq::Bool)
    type = geotype(geom)

    isequal(type, :Point) || isequal(type, :MultiPoint) && return geom

    clean!(geom)

    coords = geom.coordinates

    if isequal(type, :LineString)
        geom.coordinates = simplify_line(coords, tol, hq)
    elseif isequal(type, :MultiLineString)
        geom.coordinates = map(lines -> simplify_line(lines, tol, hq), coords[1])

    elseif isequal(type, :Polygon)
        geom.coordinates = simplify_polygon(coords, tol, hq)
    else
        geom.coordinates = map(rings -> simplify_polygon(rings, tol, hq), coords[1])
    end

    return geom
end

function simplify_line(coords, tol::Real, hq::Bool)
    return map(coords -> (coords[3]) ? [coords[1], coords[2], coords[3]] : [coords[1], coords[2]],
        DP_simplify(map(coord -> [coord[1], coord[2], coord[3]], coords), tol, hq))
end

function simplify_polygon(coords, tol::Real, hq::Bool)
    return map(ring -> begin
        pts = map(coord -> [coord[1], coord[2]], ring)

        length(pts) < 4 && throw(error("Invalid Polygon."))

        simple = map(p -> [p.coordinates[1], p.coordinates[2]], DP_simplify([Point(p) for p in pts], tol, hq))

        while !valid(simple)
            tol -= tol * 0.01
            simple = map(p -> [p.coordinates[1], p.coordinates[2]], DP_simplify([Point(p) for p in pts], tol, hq))
        end

        ((simple[length(simple)][1] != simple[1][1]) ||
            (simple[length(simple)][2] != simple[1][2])) && push!(simple, simple[1])

        return simple

        end, coords)
end

function valid(ring)
    length(ring) < 3 && return false

    return !(isequal(length(ring), 3) && (isequal(ring[3][1], ring[1][1]) && isequal(ring[3][2], ring[1][2])))
end

function dist²(p1::Point, p2::Point)
    dx = p1.coordinates[1] - p2.coordinates[1]
    dy = p1.coordinates[2] - p2.coordinates[2]

    return dx * dx + dy * dy
end

function segdist²(p::Point, l1::Point, l2::Point)
    x, y = l1.coordinates
    dx = l2.coordinates[1] - x
    dy = l2.coordinates[2] - y

    if !iszero(dx) || !iszero(dy)
        t = ((p.coordinates[1] - x) * dx + (p.coordinates[2] - y) * dy) / (dx * dx + dy * dy)

        if t > 1
            x = l2.coordinates[1]
            y = l2.coordinates[2]
        elseif t > 0
            x += dx * t
            y += dy * t
        end
    end

    dx = p.coordinates[1] - x
    dy = p.coordinates[2] - y

    return dx * dx + dy * dy
end

function simplify_radial_distance(points::Vector, tol::Real)
    previous = points[1]
    new_points = [previous]
    point = nothing

    for i in eachindex(points)
        point = points[i]
        if dist²(point, previous) > tol
            push!(new_points, point)
            previous = point
        end
    end

    !isequal(previous, point) && push!(new_points, point)

    return new_points
end

function simplify_DP_step(points::Vector, first::Integer, last::Integer, tol::Real, simplified)
    max_dist = tol
    index = 0

    for i = first+1:last
        dist = segdist²(points[i], points[first], points[last])

        if dist > max_dist
            index = i
            max_dist = dist
        end
    end

    if max_dist > tol
        (index - first > 1) && simplify_DP_step(points, first, index, tol, simplified)
        push!(simplified, points[index])
        (last - index > 1) && simplify_DP_step(points, index, last, tol, simplified)
    end
end

function douglas_peucker(points::Vector, tol::Real)
    last = length(points)

    simplified = [points[1]]
    simplify_DP_step(points, 1, last, tol, simplified)
    push!(simplified, points[last])

    return simplified
end

function DP_simplify(points::Vector, tol::Real, hq::Bool)
    length(points) <= 2 && return points
    sq_tol = tol * tol

    points = hq ? points : simplify_radial_distance(points, sq_tol)
    points = douglas_peucker(points, sq_tol)

    return points
end


text:
using Outliers
using Base.Test

a = randn(10)
o, idc = Outliers.outliers(a)
@test isempty(o)

b = [randn(10), 100]
o, idc = Outliers.outliers(b)
@test !isempty(o)

c = randn(50,2)

ext: jl


text:
#=
IO functions for Str types

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#
@inline _fastwrite(io, str) =
    @preserve str unsafe_write(io, pointer(str), reinterpret(UInt, sizeof(str)))

## outputting Str strings and Chr characters ##

# Todo: make this more common with print code

@inline _write(::Type{C}, io, str::MaybeSub{<:Str{C}}) where {C<:CSE} =
    _fastwrite(io, str)

@inline _write(::Type{LatinCSE}, io,
               str::MaybeSub{<:Str{<:Union{ASCIICSE,Binary_CSEs,_LatinCSE}}}) =
    _fastwrite(io, str)
@inline _write(::Type{_LatinCSE}, io,
               str::MaybeSub{<:Str{<:Union{ASCIICSE,Binary_CSEs,LatinCSE}}}) =
    _fastwrite(io, str)

@inline _write(::Type{UCS2CSE}, io, str::MaybeSub{<:Str{<:Union{_UCS2CSE,Text2CSE}}}) =
    _fastwrite(io, str)
@inline _write(::Type{_UCS2CSE}, io, str::MaybeSub{<:Str{<:Union{UCS2CSE,Text2CSE}}}) =
    _fastwrite(io, str)

@inline _write(::Type{UTF32CSE}, io, str::MaybeSub{<:Str{<:Union{_UTF32CSE,Text4CSE}}}) =
    _fastwrite(io, str)
@inline _write(::Type{_UTF32CSE}, io, str::MaybeSub{<:Str{<:Union{UTF32CSE,Text4CSE}}}) =
    _fastwrite(io, str)

@inline _write(::Type{UTF8CSE}, io, str::MaybeSub{<:Str{ASCIICSE}}) =
    _fastwrite(io, str)
@inline _write(::Type{UTF16CSE}, io, str::MaybeSub{<:Str{<:UCS2_CSEs}}) =
    _fastwrite(io, str)
@inline _write(::Type{UTF8CSE}, io, str::MaybeSub{<:Str{RawUTF8CSE}}) =
    _fastwrite(io, str)

function _write(::Type{<:UTF8_CSEs}, io, str::MaybeSub{<:Str{<:Union{Latin_CSEs}}})
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + sizeof(str)
        # Skip and write out ASCII sequences together
        cnt = 0
        while pnt < fin
            # Skip to first non-ASCII sequence
            # Todo: Optimize this to look at chunks at a time to find first non-ASCII
            beg = pnt
            ch = 0x00
            while (ch = get_codeunit(pnt)) < 0x80 && (pnt += 1) < fin ; end
            # Now we have from beg to < pnt that are ASCII
            unsafe_write(io, beg, pnt - beg)
            cnt += (pnt - beg)
            pnt < fin || break
            # Todo: Optimize sequences of more than one character > 0x7f
            # Write out two bytes of Latin1 character encoded as UTF-8
            _write_utf8_2(io, ch)
            cnt += 2
            pnt += 1
        end
        cnt
    end
end

function _write(::Type{<:UTF8_CSEs}, io, str::MaybeSub{<:Str{<:UCS2_CSEs}})
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + sizeof(str)
        cnt = 0
        while pnt < fin
            ch = get_codeunit(pnt)
            pnt += 2
            if ch <= 0x7f
                write(io, ch%UInt8)
                cnt += 1
            elseif ch <= 0x7ff
                _write_utf8_2(io, ch)
                cnt += 2
            else
                _write_utf8_3(io, ch)
                cnt += 3
            end
        end
        cnt
    end
end

function _write(::Type{<:UTF8_CSEs}, io, str::MaybeSub{<:Str{UTF16CSE}})
    @preserve str begin
        pnt = pointer(str)
        # Skip and write out ASCII sequences together
        fin = pnt + sizeof(str)
        cnt = 0
        while pnt < fin
            ch = get_codeunit(pnt)
            # Handle 0x80-0x7ff
            if ch <= 0x7f
                write(io, ch%UInt8)
                cnt += 1
            elseif ch <= 0x7ff
                _write_utf8_2(io, ch)
                cnt += 2
            elseif is_surrogate_lead(ch)
                _write_utf8_4(io, get_supplementary(ch, get_codeunit(pnt += 2)))
                cnt += 4
            else
                _write_utf8_3(io, ch)
                cnt += 3
            end
            pnt += 2
        end
        cnt
    end
end

function _write(::Type{<:UTF8_CSEs}, io, str::MaybeSub{<:Str{<:UTF32_CSEs}})
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + sizeof(str)
        cnt = 0
        while pnt < fin
            ch = get_codeunit(pnt)
            # Handle 0x80-0x7ff
            if ch <= 0x7f
                write(io, ch%UInt8)
                cnt += 1
            elseif ch <= 0x7ff
                _write_utf8_2(io, ch)
                cnt += 2
            elseif ch <= 0xffff
                _write_utf8_3(io, ch)
                cnt += 3
            else
                _write_utf8_4(io, ch)
                cnt += 4
            end
            pnt += 4
        end
        cnt
    end
end

function _write(::Type{UTF16CSE}, io, str::MaybeSub{<:Str{<:CSE}})
    cnt = 0
    @inbounds for ch in str
        cnt += write_utf16(io, ch)
    end
    cnt
end

function _write(::Type{T}, io, str::MaybeSub{<:Str{C}}
                ) where {T<:Union{UTF16CSE,UCS2_CSEs,UTF32_CSEs,Text2CSE,Text4CSE},
                         C<:Union{ASCIICSE,Latin_CSEs,Binary_CSEs}}
    cnt = 0
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + sizeof(str)
        while pnt < fin
            cnt += write(io, get_codeunit(pnt)%codeunit(T))
            pnt += sizeof(codeunit(C))
        end
    end
    cnt
end

function _write(::Type{UTF32CSE}, io, str::MaybeSub{<:Str{<:CSE}})
    cnt = 0
    @inbounds for ch in str
        write(io, ch%UInt32)
        cnt += 4
    end
    cnt
end

function _write(::Type{<:UTF32_CSEs}, io, str::MaybeSub{<:Str{C}}
                ) where {C<:Union{UCS2_CSEs,Text2CSE}}
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + sizeof(str)
        while pnt < fin
            write(io, get_codeunit(pnt)%UInt32)
            pnt += sizeof(codeunit(C))
        end
    end
    sizeof(str)*2
end

function _write(::Type{<:UTF32_CSEs}, io,
                str::MaybeSub{<:Str{C}}) where {C<:Union{UTF8_CSEs,UTF16CSE}}
    cnt = 0
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + sizeof(str)
        while pnt < fin
            ch, pnt = _nextcp(C, pnt)
            write(io, ch)
            cnt += 4
        end
    end
    cnt
end

@inline write(io::IO, str::MaybeSub{<:Str{C}}) where {C<:CSE} = _write(C, io, str)

# Printing bytes

print(io::IO, str::MaybeSub{<:Str{<:CSE}}) = (_write(UTF8CSE, io, str) ; nothing)

# optimized methods to avoid iterating over chars
print(io::IO, str::MaybeSub{T}) where {T<:Str{<:Union{Binary_CSEs,ASCIICSE,UTF8_CSEs}}} =
    (_fastwrite(io, str); nothing)

function Base.sprint(f::Function, ::Type{T}, args...;
                context=nothing, sizehint::Integer=0) where {T<:Union{String, Str}}
    s = IOBuffer(sizehint=sizehint)
    if context !== nothing
        f(IOContext(s, context), args...)
        S = get(context, :type, T)
        S(resize!(s.data, s.size))
    else
        f(s, args...)
        T(resize!(s.data, s.size))
    end
end

read(io::IO, ::Type{T}) where {C<:CSE,T<:Str{C}} = Str(C, read(io, String))

IOBuffer(str::T) where {T<:Str} =
    IOContext(IOBuffer(unsafe_wrap(Vector{UInt8}, str.data)), (:type => T))
IOBuffer(s::SubString{T}) where {T<:Str} =
    IOContext(IOBuffer(view(unsafe_wrap(Vector{UInt8}, s.string.data),
                            s.offset + 1 : s.offset + sizeof(s))), (:type => T))

function _joinio(io, ::Type{C}, strings) where {C}
    @inbounds for str in strings
        _write(C, io, str)
    end
    io
end

function _joinio(io, ::Type{C}, strings, delim::T) where {C,T}
    # Could speed this up in the common case where delim ===
    delbuf = (T <: AbstractString && cse(delim) === C) ? delim : convert(Str{C}, delim)
    @inbounds for str in strings
        _write(C, io, str)
        (num -= 1) == 0 || _write(C, io, delbuf)
    end
    io
end

function _joinio(io, ::Type{C}, strings, delim::T, last::S) where {C,T,S}
    # Could speed this up in the common case where delim ===
    delbuf = (T <: AbstractString && cse(delim) === C) ? delim : convert(Str{C}, delim)
    lastbuff = (S <: AbstractString && cse(last) === C) ? last : convert(Str{C}, last)
    @inbounds for str in strings
        _write(C, io, str)
        (num -= 1) > 1 ? _write(C, io, delbuf) : (num == 0 || _write(C, io, lastbuff))
    end
    io
end

function _joincvt(::Type{C}, strings) where {C}
    (num = length(strings)) < 2 && return num == 0 ? empty_str(C) : convert(Str{C}, strings[1])
    Str(C, String(take!(_joinio(IOBuffer(), C, strings))))
end

function _joincvt(::Type{C}, strings, delim) where {C}
    (num = length(strings)) < 2 && return num == 0 ? empty_str(C) : convert(Str{C}, strings[1])
    Str(C, String(take!(_joinio(IOBuffer(), C, strings, delim))))
end

function _joincvt(::Type{C}, strings, delim, last) where {C}
    (num = length(strings)) < 2 && return num == 0 ? empty_str(C) : convert(Str{C}, strings[1])
    Str(C, String(take!(_joinio(IOBuffer(), C, strings, delim, last))))
end

function _join(::Type{C}, strings) where {C}
    (num = length(strings)) < 2 && return num == 0 ? empty_str(C) : convert(Str{C}, strings[1])
    len = 0
    @inbounds for str in strings
        len += ncodeunits(str)
    end
    buf, pnt = _allocate(codeunit(C), len)
    @inbounds for str in strings
        len = ncodeunits(str)
        _memcpy(pnt, pointer(str), len)
        pnt += len
    end
    Str(C, buf)
end

function _join(::Type{C}, strings, delim::T) where {C,T}
    (num = length(strings)) < 2 && return num == 0 ? empty_str(C) : convert(Str{C}, strings[1])

    @preserve delim begin
        # Could speed this up in the common case where delim ===
        delbuf = (T <: AbstractString && cse(delim) === C) ? delim : convert(Str{C}, delim)
        delpnt = pointer(delbuf)
        dellen = ncodeunits(delbuf)
        len = 0
        @inbounds for str in strings
            len += ncodeunits(str) + dellen
        end
        len -= dellen

        buf, pnt = _allocate(codeunit(C), len)
        @inbounds for str in strings
            len = ncodeunits(str)
            _memcpy(pnt, pointer(str), len)
            pnt += len
            if (num -= 1) != 0
                _memcpy(pnt, delpnt, dellen)
                pnt += dellen
            end
        end
        Str(C, buf)
    end
end

function _join(::Type{C}, strings, delim::T, last::S) where {C,T,S}
    (num = length(strings)) < 2 && return num == 0 ? empty_str(C) : convert(Str{C}, strings[1])

    @preserve delim begin
        # Could speed this up in the common case where delim ===
        delbuf = (T <: AbstractString && cse(delim) === C) ? delim : convert(Str{C}, delim)
        delpnt = pointer(delbuf)
        dellen = ncodeunits(delbuf)
        lastbuff = (S <: AbstractString && cse(last) === C) ? last : convert(Str{C}, last)
        lastpnt = pointer(lastbuff)
        lastlen = ncodeunits(lastbuff)
        len = 0
        @inbounds for str in strings
            len += ncodeunits(str) + dellen
        end
        len = len - 2*dellen + lastlen

        buf, pnt = _allocate(codeunit(C), len)
        @inbounds for str in strings
            len = ncodeunits(str)
            _memcpy(pnt, pointer(str), len)
            pnt += len
            if (num -= 1) > 1
                _memcpy(pnt, delpnt, dellen)
                pnt += dellen
            elseif num != 0
                _memcpy(pnt, lastpnt, lastlen)
                pnt += lastlen
            end
        end
        Str(C, buf)
    end
end

@inline function calc_type(strings)
    C = Union{}
    for str in strings
        C = promote_type(C, cse(str))
    end
    C
end

join(io::IO, strings::AbstractArray{<:MaybeSub{<:Str}}) =
    (_joinio(io, calc_type(strings), strings) ; nothing)
join(io::IO, strings::AbstractArray{<:MaybeSub{<:Str}}, delim) =
    (_joinio(io, calc_type(strings), strings, delim) ; nothing)
join(io::IO, strings::AbstractArray{<:MaybeSub{<:Str}}, delim, last) =
    (_joinio(io, calc_type(strings), strings, delim, last) ; nothing)

join(strings::AbstractArray{<:MaybeSub{<:Str}}) =
    _joincvt(calc_type(strings), strings)
join(strings::AbstractArray{<:MaybeSub{T}}) where {C<:Union{ASCIICSE, Latin_CSEs},T<:Str{C}} =
    _join(C, strings)
join(strings::AbstractArray{<:MaybeSub{T}}) where {C<:Word_CSEs,T<:Str{C}} =
    _join(C, strings)
join(strings::AbstractArray{<:MaybeSub{T}}) where {C<:Quad_CSEs,T<:Str{C}} =
    _join(C, strings)

join(strings::AbstractArray{<:MaybeSub{<:Str}}, delim) =
    _joincvt(_calc_type(strings), strings, delim)
join(strings::AbstractArray{<:MaybeSub{T}},
     delim) where {C<:Union{Text1CSE, BinaryCSE, ASCIICSE, Latin_CSEs},T<:Str{C}} =
    _join(C, strings, delim)
join(strings::AbstractArray{<:MaybeSub{T}}, delim) where {C<:Word_CSEs,T<:Str{C}} =
    _join(C, strings, delim)
join(strings::AbstractArray{<:MaybeSub{T}}, delim) where {C<:Quad_CSEs,T<:Str{C}} =
    _join(C, strings, delim)

join(strings::AbstractArray{<:MaybeSub{<:Str}}, delim, last) =
    _joincvt(calc_type(strings), strings, delim, last)
join(strings::AbstractArray{<:MaybeSub{T}},
     delim, last) where {C<:Union{Text1CSE, BinaryCSE, ASCIICSE, Latin_CSEs},T<:Str{C}} =
         _join(C, strings, delim, last)
join(strings::AbstractArray{<:MaybeSub{T}},
     delim, last) where {C<:Word_CSEs,T<:Str{C}} =
         _join(C, strings, delim, last)
join(strings::AbstractArray{<:MaybeSub{T}},
     delim, last) where {C<:Quad_CSEs,T<:Str{C}} =
         _join(C, strings, delim, last)


text:
using ADCME
using LinearAlgebra
using AdFem
using PyPlot
using SparseArrays

function Q_exact(x, y)
    # (2*x^2*y*(x - 1)*(2*y - 1) + 1.0*(y^2 + 1)^2*(x*y*(x - 1) + x*y*(y - 1) + x*(x - 1)*(y - 1) + y*(x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
    # (2*x^2*y*(x - 1)*(2*y - 1) + x*(y^2 + 1)^2*(-y^2*(2*x - 1)*(y - 1)^2 + 1.0*y*(x - 1) + 1.0*(x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
    (2*x^2*y*(x - 1)*(2*y - 1)*(x + y + 1) - x*y^2*(2*x - 1)*(y - 1)^2*(y^2 + 1)^2*(x + y + 1) + x*(x - 1)*(2*y - 1)*(y^2 + 1)^2 - (y^2 + 1)*(x + y + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/((y^2 + 1)^2*(x + y + 1))
end

function T_exact(x, y)
    x * (1-x) * y * (1-y)
end

function u_exact(x, y)
    x * y * (1-y)
end

function v_exact(x, y)
    1 / (1 + x + y)
end

function k_exact(x, y)
    1 + x^2 + x / (1 + y^2)
end

# function k_nn(xy) # xy shape N (=m*n) x 2
#     out = fc(xy, [20,20,20,1])^2 + 0.1 # N x 1
#     squeeze(out)
# end

#---------------------------------------------
# grid setup

######## mesh 1 ########
# m = 20
# n = 20
# h = 1/n
# mesh = Mesh(m, n, h)

######## mesh 2 ########
filename = "CHT_2D.stl"
file_format = "stl"
mesh = Mesh(filename, file_format = file_format)
mesh = Mesh(mesh.nodes ./ 0.0305, mesh.elems)

nnode = size(mesh.nodes, 1)
nelem = size(mesh.elems, 1)

#---------------------------------------------
# discretized governing equation
#  J * u^{n+1} - J * u^n = - v ⋅ grad u + K * u^{n+1} + F^{n+1}
#---------------------------------------------
xy = mesh.nodes
x, y = xy[:,1], xy[:,2]
T0 = @. T_exact(x, y)
u = @. u_exact(x,y)
v = @. v_exact(x,y)
Q = @. Q_exact(x,y)
k = @. k_exact(x, y)
# k0 = @. k_nn(xy)

# ---------------------------------------------

bd = Array{Int64, 1}([])
eps = 1e-6
for j = 1:nnode
    nodex, nodey = mesh.nodes[j, 1], mesh.nodes[j, 2]
    if abs(nodex-0.0) <= eps || abs(nodex-1.0) <= eps || abs(nodey-0.0) <= eps || abs(nodey-1.0) <= eps
        global bd = [bd; j]
    end
end

ugauss = fem_to_gauss_points(u, mesh)
vgauss = fem_to_gauss_points(v, mesh)
kgauss = fem_to_gauss_points(k, mesh)
Qgauss = fem_to_gauss_points(Q, mesh)

Advection = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh))
Laplace = compute_fem_laplace_matrix1(constant(kgauss), mesh)
A = Advection + Laplace
A, _ = fem_impose_Dirichlet_boundary_condition1(A, bd, mesh)
b = constant(compute_fem_source_term1(constant(Qgauss), mesh))
b = scatter_update(b, bd, zeros(length(bd)))
sol = A\b

# ---------------------------------------------------
# create a session and run
sess = Session(); init(sess)
T_computed = run(sess, sol)

#---------------------------------------------
# visualize numerical solution and exact solution
figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fem_points(T_computed, mesh, with_mesh=true)
subplot(122)
visualize_scalar_on_fem_points(T0, mesh, with_mesh=true)
savefig("forward_solution_unstructured.png")
close("all")


text:
module SigMF

using JSON3
using Tar
using SHA
using StructTypes
using Dates
using Mmap

const _sigmf_default_version = "0.0.2"

include("errors.jl")
include("datatype.jl")
include("extension_signal.jl")
include("core.jl")
include("io.jl")

export SigMFFile, SigMFMeta, SigMFGlobal, SigMFCaptures, SigMFAnnotations

end


text:
export Env, os
"""
Custom type to store *os-environment sensitive parameters* and commands
"""
type Env
    newline::AbstractString
    xppcall::AbstractString
    options::Array
end

#set os-dependent variables
os = Dict()
os[:Linux] = Env("\n", "xppaut", ["-silent",  "-logfile",  "log.txt"])
os[:Windows] = Env("\r\n", "/xppall/xppaut", ["-silent",  "-logfile",  "log.txt"])
os[:Darwin] = Env("\n", "xppaut", ["-silent",  "-logfile",  "log.txt"])

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
module Database

import Revise
import YAML, DataFrames, Logging
using SearchLight, SearchLight.Configuration

const ADAPTERS_FOLDER = joinpath(@__DIR__, "database_adapters")

function setup_adapter(adapter = SearchLight.config.db_config_settings["adapter"] * "DatabaseAdapter") :: Union{Bool,Nothing}
  Core.eval(@__MODULE__, Meta.parse("""include(joinpath(ADAPTERS_FOLDER, "$adapter.jl"))"""))
  Core.eval(@__MODULE__, Meta.parse("using .$adapter"))

  Core.eval(@__MODULE__, :(db_adapter = Symbol($adapter)))
  ! Core.isdefined(@__MODULE__, :DatabaseAdapter) && Core.eval(@__MODULE__, :(const DatabaseAdapter = $db_adapter))
  ! Core.isdefined(@__MODULE__, :DatabaseHandle) && Core.eval(@__MODULE__, :(const DatabaseHandle = $db_adapter.DatabaseHandle))
  ! Core.isdefined(@__MODULE__, :ResultHandle) && Core.eval(@__MODULE__, :(const ResultHandle = $db_adapter.ResultHandle))

  Core.eval(@__MODULE__, :(export DatabaseAdapter))

  true
end

"""
    connect()::DatabaseHandle
    connect(conn_settings::Dict{String,Any})::DatabaseHandle

Connects to the DB and returns a database handler. If used without arguments, it defaults to using `config.db_config_settings`

# Examples
```julia
julia> Database.connect()
PostgreSQL.PostgresDatabaseHandle(Ptr{Nothing} @0x00007fbf3839f360,0x00000000,false)

julia> dict = config.db_config_settings
Dict{String,Any} with 6 entries:
  "host"     => "localhost"
  "password" => "adrian"
  "username" => "adrian"
  "port"     => 5432
  "database" => "blogjl_dev"
  "adapter"  => "PostgreSQL"

julia> Database.connect(dict)
PostgreSQL.PostgresDatabaseHandle(Ptr{Nothing} @0x00007fbf3839f360,0x00000000,false)
```
"""
function connect() #::DatabaseHandle
  connect(SearchLight.config.db_config_settings)
end

function connect(conn_settings::Dict)
  isdefined(@__MODULE__, :DatabaseAdapter) || connect!(conn_settings)

  c = Base.invokelatest(DatabaseAdapter.connect, conn_settings) #::DatabaseHandle
  Core.eval(@__MODULE__, :(CONNECTION = $c))

  c
end

function connect!(conn_settings::Dict)
  SearchLight.config.db_config_settings["adapter"] = conn_settings["adapter"]
  setup_adapter() && Database.connect(conn_settings) #::DatabaseHandle
end

"""

"""
function disconnect(conn)
  DatabaseAdapter.disconnect(conn)
end

function disconnect()
  DatabaseAdapter.disconnect(CONNECTION)
end

"""
    create_database()::Bool
    create_database(db_name::String)::Bool

Invokes the database adapter's create database method. If invoked without param, it defaults to the
database name defined in `config.db_config_settings`
"""
function create_database() :: Bool
  create_database(SearchLight.config.db_config_settings["database"])
end

function create_database(db_name::String)::Bool
  DatabaseAdapter.create_database(db_name)
end

"""
    create_migrations_table(table_name::String)::Bool

Invokes the database adapter's create migrations table method. If invoked without param, it defaults to the
database name defined in `config.db_migrations_table_name`
"""
function create_migrations_table(table_name::String) :: Bool
  DatabaseAdapter.create_migrations_table(table_name)
end

"""
    init()::Bool

Sets up the DB tables used by SearchLight.
"""
function init() :: Bool
  DatabaseAdapter.create_migrations_table(SearchLight.config.db_migrations_table_name)
end

function escape_column_name(c::String) :: String
  DatabaseAdapter.escape_column_name(c, CONNECTION)
end

function escape_value(v::T)::T where {T}
  result =  try
              DatabaseAdapter.escape_value(v, CONNECTION)
            catch ex
              @error ex
            end

  result
end

function table_columns(table_name::String) :: DataFrames.DataFrame
  query(DatabaseAdapter.table_columns_sql(table_name), suppress_output = true)
end

"""
    query(sql::String; suppress_output::Bool = false)::DataFrames.DataFrame

Executes the `sql` query against the database adapter and returns a DataFrame result.
Optionally logs the result DataFrame.

# Examples:
```julia
julia> SearchLight.to_fetch_sql(Article, SQLQuery(limit = 5)) |> Database.query;

2017-01-16T21:33:40.079 - info: SQL QUERY: SELECT \"articles\".\"id\" AS \"articles_id\", \"articles\".\"title\" AS \"articles_title\", \"articles\".\"summary\" AS \"articles_summary\", \"articles\".\"content\" AS \"articles_content\", \"articles\".\"updated_at\" AS \"articles_updated_at\", \"articles\".\"published_at\" AS \"articles_published_at\", \"articles\".\"slug\" AS \"articles_slug\" FROM \"articles\" LIMIT 5

  0.001172 seconds (16 allocations: 576 bytes)

2017-01-16T21:33:40.089 - info: 5×7 DataFrames.DataFrame
...
```
"""
function query(sql::String; suppress_output::Bool = false, system_query::Bool = false) :: DataFrames.DataFrame
  df::DataFrames.DataFrame =  DatabaseAdapter.query(sql, (suppress_output || system_query), CONNECTION)
  (! suppress_output && ! system_query && SearchLight.config.log_db) && @info(df)

  df
end

"""

"""
function to_find_sql(m::Type{T}, q::SearchLight.SQLQuery, joins::Union{Nothing,Vector{SearchLight.SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  DatabaseAdapter.to_find_sql(m, q, joins)
end

const to_fetch_sql = to_find_sql

"""

"""
function to_store_sql(m::T; conflict_strategy = :error)::String where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.to_store_sql(m, conflict_strategy = conflict_strategy)
end

"""

"""
function delete_all(m::Type{T}; truncate::Bool = true, reset_sequence::Bool = true, cascade::Bool = false)::Nothing where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.delete_all(m, truncate = truncate, reset_sequence = reset_sequence, cascade = cascade)
end

"""

"""
function delete(m::T)::T where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.delete(m)
end

"""

"""
function count(m::Type{T}, q::SearchLight.SQLQuery = SearchLight.SQLQuery())::Int where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.count(m, q)
end

"""

"""
function update_query_part(m::T)::String where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.update_query_part(m)
end

"""

"""
function to_select_part(m::Type{T}, cols::Vector{SearchLight.SQLColumn}, joins::Union{Nothing,Vector{SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  DatabaseAdapter.to_select_part(m, cols, joins)
end
"""

"""
function _to_select_part(m::Type{T}, cols::Vector{SearchLight.SQLColumn}, joins::Union{Nothing,Vector{SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  _m::T = m()

  if ! isempty(cols)
    table_columns = []
    cols = vcat(cols, columns_from_joins(joins))

    for column in cols
      push!(table_columns, prepare_column_name(column, _m))
    end

    return join(table_columns, ", ")
  else
    tbl_cols = join(SearchLight.to_fully_qualified_sql_column_names(_m, SearchLight.persistable_fields(_m), escape_columns = true), ", ")
    table_columns = isempty(tbl_cols) ? String[] : vcat(tbl_cols, map(x -> prepare_column_name(x, _m), columns_from_joins(joins)))

    return join(table_columns, ", ")
  end
end

"""

"""
function to_from_part(m::Type{T})::String where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.to_from_part(m)
end

"""

"""
function to_where_part(w::Vector{SearchLight.SQLWhereEntity})::String
  DatabaseAdapter.to_where_part(w)
end

"""

"""
function to_order_part(m::Type{T}, o::Vector{SearchLight.SQLOrder})::String where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.to_order_part(m, o)
end

"""

"""
function to_group_part(group::Vector{SearchLight.SQLColumn}) :: String
  DatabaseAdapter.to_group_part(group)
end

"""

"""
function to_limit_part(limit::SearchLight.SQLLimit) :: String
  DatabaseAdapter.to_limit_part(limit)
end

"""

"""
function to_offset_part(offset::Int) :: String
  DatabaseAdapter.to_offset_part(offset)
end

"""

"""
function to_having_part(having::Vector{SearchLight.SQLHaving}) :: String
  DatabaseAdapter.to_having_part(having)
end

"""

"""
function to_join_part(m::Type{T}, joins::Union{Nothing,Vector{SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  DatabaseAdapter.to_join_part(m, joins)
end

"""
  columns_from_joins(joins::Vector{SQLJoin})::Vector{SQLColumn}

Extracts columns from joins param and adds to be used for the SELECT part
"""
function columns_from_joins(joins::Union{Nothing,Vector{SQLJoin{N}}} = nothing)::Vector{SearchLight.SQLColumn} where {N<:Union{Nothing,SearchLight.AbstractModel}}
  jcols = SearchLight.SQLColumn[]

  joins === nothing && return jcols

  for j in joins
    jcols = vcat(jcols, j.columns)
  end

  jcols
end

"""

"""
function prepare_column_name(column::SearchLight.SQLColumn, _m::T)::String where {T<:SearchLight.AbstractModel}
  if column.raw
    column.value |> string
  else
    column_data::Dict{Symbol,Any} = SearchLight.from_literal_column_name(column.value)
    if ! haskey(column_data, :table_name)
      column_data[:table_name] = SearchLight.table_name(_m)
    end
    if ! haskey(column_data, :alias)
      column_data[:alias] = ""
    end

    DatabaseAdapter.column_data_to_column_name(column, column_data)
  end
end

"""

"""
function rand(m::Type{T}; limit = 1)::Vector{T} where {T<:SearchLight.AbstractModel}
  DatabaseAdapter.rand(m, limit = limit)
end

"""

"""
function index_name(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}) :: String
  string(table_name) * "__" * "idx_" * string(column_name)
end

end


text:
using Jig

tally = 0
include("foo.jl")
tally += fails
include("bar.jl")
tally += fails

tally > 0 ? error("There are $tally failing tests in this suite") : nothing


text:
module FillArrays

using LinearAlgebra, SparseArrays
import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
                +, -, *, /, \, sum, cumsum, maximum, minimum
import LinearAlgebra: rank

import Base.Broadcast: broadcasted, DefaultArrayStyle

export Zeros, Ones, Fill, Eye

abstract type AbstractFill{T, N} <: AbstractArray{T, N} end

@inline function getindex(F::AbstractFill, k::Integer)
    @boundscheck checkbounds(F, k)
    getindex_value(F)
end

@inline function getindex(F::AbstractFill{T, N}, kj::Vararg{<:Integer, N}) where {T, N}
    @boundscheck checkbounds(F, kj...)
    getindex_value(F)
end

rank(F::AbstractFill) = iszero(getindex_value(F)) ? 0 : 1

IndexStyle(::Type{<:AbstractFill}) = IndexLinear()

struct Fill{T, N, SZ} <: AbstractFill{T, N}
    value::T
    size::SZ

    @inline function Fill{T,N,SZ}(x::T, sz::SZ) where SZ<:Tuple{Vararg{Integer,N}} where {T, N}
        @boundscheck any(k -> k < 0, sz) && throw(BoundsError())
        new{T,N,SZ}(x,sz)
    end
    @inline Fill{T, N}(x::T, sz::SZ) where SZ<:Tuple{Vararg{Integer,N}} where {T, N} = Fill{T,N,SZ}(x, sz)
    @inline Fill{T, N}(x, sz::SZ) where SZ<:Tuple{Vararg{Integer,N}} where {T, N} = Fill{T,N}(convert(T, x)::T, sz)
    @inline Fill{T, N}(x, sz::Vararg{Integer, N}) where {T, N} = Fill{T,N}(convert(T, x)::T, sz)
end

@inline Fill{T}(x, sz::Vararg{<:Integer,N}) where {T, N} = Fill{T, N}(x, sz)
@inline Fill{T}(x, sz::Tuple{Vararg{Integer,N}}) where {T, N} = Fill{T, N}(x, sz)
@inline Fill(x::T, sz::Vararg{<:Integer,N}) where {T, N}  = Fill{T, N}(x, sz)
@inline Fill(x::T, sz::Tuple{Vararg{Integer,N}}) where {T, N}  = Fill{T, N}(x, sz)

@inline size(F::Fill) = F.size
@inline getindex_value(F::Fill) = F.value

AbstractArray{T}(F::Fill{T}) where T = F
AbstractArray{T,N}(F::Fill{T,N}) where {T,N} = F
AbstractArray{T}(F::Fill{V,N}) where {T,V,N} = Fill{T}(convert(T, F.value)::T, F.size)
AbstractArray{T,N}(F::Fill{V,N}) where {T,V,N} = Fill{T}(convert(T, F.value)::T, F.size)

convert(::Type{AbstractArray{T}}, F::Fill{T}) where T = F
convert(::Type{AbstractArray{T,N}}, F::Fill{T,N}) where {T,N} = F
convert(::Type{AbstractArray{T}}, F::Fill) where {T} = AbstractArray{T}(F)
convert(::Type{AbstractArray{T,N}}, F::Fill) where {T,N} = AbstractArray{T,N}(F)
convert(::Type{AbstractFill}, F::AbstractFill) = F
convert(::Type{AbstractFill{T}}, F::AbstractFill) where T = convert(AbstractArray{T}, F)
convert(::Type{AbstractFill{T,N}}, F::AbstractFill) where {T,N} = convert(AbstractArray{T,N}, F)

function getindex(F::Fill, kj::Vararg{AbstractVector{II},N}) where {II<:Integer,N}
    checkbounds(F, kj...)
    Fill(getindex_value(F), length.(kj))
end

function getindex(A::Fill, kr::AbstractVector{Bool})
   length(A) == length(kr) || throw(DimensionMismatch())
   Fill(getindex_value(A), count(kr))
end
function getindex(A::Fill, kr::AbstractArray{Bool})
   size(A) == size(kr) || throw(DimensionMismatch())
   Fill(getindex_value(A), count(kr))
end

+(a::AbstractFill) = a
-(a::AbstractFill) = Fill(-getindex_value(a), size(a))

# Fill +/- Fill
function +(a::AbstractFill{T, N}, b::AbstractFill{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return Fill(getindex_value(a) + getindex_value(b), size(a))
end
-(a::AbstractFill, b::AbstractFill) = a + (-b)

function +(a::Fill{T, 1}, b::AbstractRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return (a.value + first(b)):convert(Tout, step(b)):(a.value + last(b))
end
function +(a::Fill{T, 1}, b::UnitRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return (a.value + first(b)):(a.value + last(b))
end
+(a::AbstractRange, b::AbstractFill) = b + a

-(a::AbstractFill, b::AbstractRange) = a + (-b)
-(a::AbstractRange, b::AbstractFill) = a + (-b)

for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
    @eval begin
        struct $Typ{T, N, SZ} <: AbstractFill{T, N}
            size::SZ
            @inline function $Typ{T, N}(sz::SZ) where SZ<:Tuple{Vararg{Integer,N}} where {T, N}
                @boundscheck any(k -> k < 0, sz) && throw(BoundsError())
                new{T,N,SZ}(sz)
            end
            @inline $Typ{T, N}(sz::Vararg{<:Integer, N}) where {T, N} = $Typ(sz)
        end

        @inline $Typ{T}(sz::Vararg{Integer,N}) where {T, N} = $Typ{T, N}(sz)
        @inline $Typ{T}(sz::SZ) where SZ<:Tuple{Vararg{Integer,N}} where {T, N} = $Typ{T, N}(sz)
        @inline $Typ(sz::Vararg{Integer,N}) where N = $Typ{Float64,N}(sz)
        @inline $Typ(sz::SZ) where SZ<:Tuple{Vararg{Integer,N}} where N = $Typ{Float64,N}(sz)

        @inline $Typ{T,N}(A::AbstractArray{V,N}) where{T,V,N} = $Typ{T,N}(size(A))
        @inline $Typ{T}(A::AbstractArray) where{T} = $Typ{T}(size(A))
        @inline $Typ(A::AbstractArray) = $Typ(size(A))

        @inline size(Z::$Typ) = Z.size
        @inline getindex_value(Z::$Typ{T}) where T = $func(T)

        AbstractArray{T}(F::$Typ{T}) where T = F
        AbstractArray{T,N}(F::$Typ{T,N}) where {T,N} = F
        AbstractArray{T}(F::$Typ) where T = $Typ{T}(F.size)
        AbstractArray{T,N}(F::$Typ{V,N}) where {T,V,N} = $Typ{T}(F.size)
        convert(::Type{AbstractArray{T}}, F::$Typ{T}) where T = AbstractArray{T}(F)
        convert(::Type{AbstractArray{T,N}}, F::$Typ{T,N}) where {T,N} = AbstractArray{T,N}(F)
        convert(::Type{AbstractArray{T}}, F::$Typ) where T = AbstractArray{T}(F)
        convert(::Type{AbstractArray{T,N}}, F::$Typ) where {T,N} = AbstractArray{T,N}(F)

        function getindex(F::$Typ{T}, kj::Vararg{AbstractVector{II},N}) where {T,II<:Integer,N}
            checkbounds(F, kj...)
            $Typ{T}(length.(kj))
        end
        function getindex(A::$Typ{T}, kr::AbstractVector{Bool}) where T
            length(A) == length(kr) || throw(DimensionMismatch())
            $Typ{T}(count(kr))
        end
        function getindex(A::$Typ{T}, kr::AbstractArray{Bool}) where T
            size(A) == size(kr) || throw(DimensionMismatch())
            $Typ{T}(count(kr))
        end
    end
end

rank(F::Zeros) = 0
rank(F::Ones) = 1

struct Eye{T, SZ} <: AbstractMatrix{T}
    size::SZ
    @inline function Eye{T}(sz::SZ) where {T,SZ<:Tuple{Vararg{Integer,2}}}
        @boundscheck any(k -> k < 0, sz) && throw(BoundsError())
        new{T,SZ}(sz)
    end

    Eye{T}(sz::Vararg{Integer,2}) where {T} = Eye{T}(sz)
end

Eye{T}(n::Integer) where T = Eye{T}(n, n)
Eye(n::Integer, m::Integer) = Eye{Float64}(n, m)
Eye(sz::Tuple{Vararg{Integer,2}}) = Eye{Float64}(sz)
Eye(n::Integer) = Eye(n, n)

@inline Eye{T}(A::AbstractMatrix) where T = Eye{T}(size(A))
@inline Eye(A::AbstractMatrix) = Eye{eltype(A)}(size(A))

size(E::Eye) = E.size
rank(E::Eye) = minimum(size(E))

@inline function getindex(E::Eye{T}, k::Integer, j::Integer) where T
    @boundscheck checkbounds(E, k, j)
    ifelse(k == j, one(T), zero(T))
end

IndexStyle(::Type{<:Eye}) = IndexCartesian()

AbstractArray{T}(E::Eye{T}) where T = E
AbstractMatrix{T}(E::Eye{T}) where T = E
AbstractArray{T}(E::Eye) where T = Eye{T}(E.size)
AbstractMatrix{T}(E::Eye) where T = Eye{T}(E.size)
convert(::Type{AbstractArray{T}}, E::Eye{T}) where T = E
convert(::Type{AbstractMatrix{T}}, E::Eye{T}) where T = E
convert(::Type{AbstractArray{T}}, E::Eye) where T = Eye{T}(E.size)
convert(::Type{AbstractMatrix{T}}, E::Eye) where T = Eye{T}(E.size)

#########
#  Special matrix types
#########

## Array
convert(::Type{Array}, F::AbstractFill) = fill(getindex_value(F), size(F))
convert(::Type{Array{T}}, F::AbstractFill) where T = fill(convert(T, getindex_value(F)), size(F))
convert(::Type{Array{T,N}}, F::AbstractFill{V,N}) where {T,V,N} = fill(convert(T, getindex_value(F)), size(F))

# These are in case `zeros` or `ones` are ever faster than `fill`
for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
    @eval begin
        convert(::Type{Array}, F::$Typ{T}) where T = $funcs(T, size(F))
        convert(::Type{Array{T}}, F::$Typ{T}) where T = $funcs(T, size(F))
        convert(::Type{Array{T,N}}, F::$Typ{V,N}) where {T,V,N} = $funcs(T,size(F))
    end
end

convert(::Type{Array},     E::Eye{T}) where T = Matrix{T}(I, E.size[1], E.size[2])
convert(::Type{Array{T}},  E::Eye)    where T = Matrix{T}(I, E.size[1], E.size[2])
convert(::Type{Matrix{T}}, E::Eye)    where T = Matrix{T}(I, E.size[1], E.size[2])

function convert(::Type{Diagonal}, Z::Zeros{T,2}) where T
    n,m = size(Z)
    n ≠ m && throw(BoundsError(Z))
    Diagonal(zeros(T, n))
end

function convert(::Type{Diagonal{T}}, Z::Zeros{V,2}) where {T,V}
    n,m = size(Z)
    n ≠ m && throw(BoundsError(Z))
    Diagonal(zeros(T, n))
end

function convert(::Type{Diagonal}, E::Eye{T}) where T
    n,m = size(E)
    n ≠ m && throw(BoundsError(E))
    Diagonal(ones(T, n))
end

function convert(::Type{Diagonal{T}}, E::Eye{V}) where {T,V}
    n,m = size(E)
    n ≠ m && throw(BoundsError(E))
    Diagonal(ones(T, n))
end

## Sparse arrays

convert(::Type{SparseVector}, Z::Zeros{T,1}) where T = spzeros(T, length(Z))
convert(::Type{SparseVector{Tv}}, Z::Zeros{T,1}) where {T,Tv} = spzeros(Tv, length(Z))
convert(::Type{SparseVector{Tv,Ti}}, Z::Zeros{T,1}) where {T,Tv,Ti} = spzeros(Tv, Ti, length(Z))

convert(::Type{AbstractSparseVector}, Z::Zeros{T,1}) where T = spzeros(T, length(Z))
convert(::Type{AbstractSparseVector{Tv}}, Z::Zeros{T,1}) where {Tv,T}= spzeros(Tv, length(Z))

convert(::Type{SparseMatrixCSC}, Z::Zeros{T,2}) where T = spzeros(T, size(Z)...)
convert(::Type{SparseMatrixCSC{Tv}}, Z::Zeros{T,2}) where {T,Tv} = spzeros(Tv, size(Z)...)
convert(::Type{SparseMatrixCSC{Tv,Ti}}, Z::Zeros{T,2}) where {T,Tv,Ti} = spzeros(Tv, Ti, size(Z)...)

convert(::Type{AbstractSparseMatrix}, Z::Zeros{T,2}) where {T} = spzeros(T, size(Z)...)
convert(::Type{AbstractSparseMatrix{Tv}}, Z::Zeros{T,2}) where {T,Tv} = spzeros(Tv, size(Z)...)

convert(::Type{AbstractSparseArray}, Z::Zeros{T}) where T = spzeros(T, size(Z)...)
convert(::Type{AbstractSparseArray{Tv}}, Z::Zeros{T}) where {T,Tv} = spzeros(Tv, size(Z)...)
convert(::Type{AbstractSparseArray{Tv,Ti}}, Z::Zeros{T}) where {T,Tv,Ti} = spzeros(Tv, Ti, size(Z)...)
convert(::Type{AbstractSparseArray{Tv,Ti,N}}, Z::Zeros{T,N}) where {T,Tv,Ti,N} = spzeros(Tv, Ti, size(Z)...)

convert(::Type{SparseMatrixCSC}, Z::Eye{T}) where T = SparseMatrixCSC{T}(I, size(Z)...)
convert(::Type{SparseMatrixCSC{Tv}}, Z::Eye{T}) where {T,Tv} = SparseMatrixCSC{Tv}(I, size(Z)...)
# works around missing `speye`:
convert(::Type{SparseMatrixCSC{Tv,Ti}}, Z::Eye{T}) where {T,Tv,Ti<:Integer} =
    convert(SparseMatrixCSC{Tv,Ti}, SparseMatrixCSC{Tv}(I, size(Z)...))

convert(::Type{AbstractSparseMatrix}, Z::Eye{T}) where {T} = SparseMatrixCSC{T}(I, size(Z)...)
convert(::Type{AbstractSparseMatrix{Tv}}, Z::Eye{T}) where {T,Tv} = SparseMatrixCSC{Tv}(I, size(Z)...)

convert(::Type{AbstractSparseArray}, Z::Eye{T}) where T = SparseMatrixCSC{T}(I, size(Z)...)
convert(::Type{AbstractSparseArray{Tv}}, Z::Eye{T}) where {T,Tv} = SparseMatrixCSC{Tv}(I, size(Z)...)

convert(::Type{AbstractSparseArray{Tv,Ti}}, Z::Eye{T}) where {T,Tv,Ti} =
    convert(SparseMatrixCSC{Tv,Ti}, Z)
convert(::Type{AbstractSparseArray{Tv,Ti,2}}, Z::Eye{T}) where {T,Tv,Ti} =
    convert(SparseMatrixCSC{Tv,Ti}, Z)

## Algebraic identities

function mult_zeros(a, b::AbstractMatrix)
    size(a, 2) ≠ size(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Zeros{promote_type(eltype(a), eltype(b))}(size(a, 1), size(b, 2))
end
function mult_zeros(a, b::AbstractVector)
    size(a, 2) ≠ size(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Zeros{promote_type(eltype(a), eltype(b))}(size(a, 1))
end

const ZerosVecOrMat{T} = Union{Zeros{T,1}, Zeros{T,2}}
*(a::ZerosVecOrMat, b::AbstractMatrix) = mult_zeros(a, b)
*(a::AbstractMatrix, b::ZerosVecOrMat) = mult_zeros(a, b)
*(a::ZerosVecOrMat, b::AbstractVector) = mult_zeros(a, b)
*(a::AbstractVector, b::ZerosVecOrMat) = mult_zeros(a, b)
*(a::ZerosVecOrMat, b::ZerosVecOrMat) = mult_zeros(a, b)

function *(a::Adjoint{T, <:AbstractVector{T}}, b::Zeros{S, 1}) where {T, S}
    la, lb = length(a), length(b)
    if la ≠ lb
        throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
    end
    return zero(promote_type(T, S))
end
*(a::Adjoint{T, <:AbstractMatrix{T}} where T, b::Zeros{<:Any, 1}) = mult_zeros(a, b)

function *(a::Transpose{T, <:AbstractVector{T}}, b::Zeros{T, 1}) where T<:Real
    la, lb = length(a), length(b)
    if la ≠ lb
        throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
    end
    return zero(T)
end
*(a::Transpose{T, <:AbstractMatrix{T}}, b::Zeros{T, 1}) where T<:Real = mult_zeros(a, b)

+(a::Zeros) = a
-(a::Zeros) = a

# Zeros +/- Zeros
function +(a::Zeros{T}, b::Zeros{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return Zeros{promote_type(T,V)}(size(a)...)
end
-(a::Zeros, b::Zeros) = -(a + b)

# Zeros +/- Fill and Fill +/- Zeros
function +(a::AbstractFill{T}, b::Zeros{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return convert(AbstractFill{promote_type(T, V)}, a)
end
+(a::Zeros, b::AbstractFill) = b + a
-(a::AbstractFill, b::Zeros) = a + b
-(a::Zeros, b::AbstractFill) = a + (-b)

# Zeros +/- Array and Array +/- Zeros
function +(a::Zeros{T, N}, b::Array{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractArray{promote_type(T,V),N}(b)
end
function +(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractArray{promote_type(T,V),N}(a)
end

function -(a::Zeros{T, N}, b::Array{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N} = a + b

+(a::AbstractRange, b::Zeros) = b + a

function +(a::Zeros{T, 1}, b::AbstractRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return convert(Tout, first(b)):convert(Tout, step(b)):convert(Tout, last(b))
end
function +(a::Zeros{T, 1}, b::UnitRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return convert(Tout, first(b)):convert(Tout, last(b))
end

function -(a::Zeros{T, 1}, b::AbstractRange{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::AbstractRange{T}, b::Zeros{V, 1}) where {T, V} = a + b

#########
# broadcasting
#########

for op in (:+, :-)
    @eval broadcasted(::DefaultArrayStyle{N}, ::typeof($op), r1::AbstractFill{T,N}, r2::AbstractFill{V,N}) where {T,V,N} =
        $op(r1, r2)
end

broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}) where {T,N} = Fill(op(getindex_value(r)), size(r))
broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}, x::Number) where {T,N} = Fill(op(getindex_value(r),x), size(r))
broadcasted(::DefaultArrayStyle{N}, op, x::Number, r::AbstractFill{T,N}) where {T,N} = Fill(op(x, getindex_value(r)), size(r))
function broadcasted(::DefaultArrayStyle{N}, op, r1::AbstractFill{T,N}, r2::AbstractFill{V,N}) where {T,V,N}
    size(r1) ≠ size(r2) && throw(DimensionMismatch("dimensions must match."))
    Fill(op(getindex_value(r1),getindex_value(r2)), size(r1))
end

for op in (:*, :/, :\)
    @eval function broadcasted(::DefaultArrayStyle{N}, ::typeof($op), r1::Ones{T,N}, r2::Ones{V,N}) where {T,V,N}
        size(r1) ≠ size(r2) && throw(DimensionMismatch("dimensions must match."))
        Ones{promote_type(T,V)}(size(r1))
    end
end

#########
# maximum/minimum
#########

for op in (:maximum, :minimum)
    @eval $op(x::AbstractFill) = getindex_value(x)
end

#########
# Cumsum
#########

sum(x::AbstractFill) = getindex_value(x)*length(x)
sum(x::Zeros) = getindex_value(x)

cumsum(x::AbstractFill) = range(getindex_value(x); step=getindex_value(x),
                                                    length=length(x))

cumsum(x::Zeros) = x
cumsum(x::Zeros{Bool}) = x
cumsum(x::Ones{II}) where II<:Integer = Base.OneTo{II}(length(x))
cumsum(x::Ones{Bool}) = Base.OneTo{Int}(length(x))
cumsum(x::AbstractFill{Bool}) = cumsum(convert(AbstractFill{Int}, x))

end # module


text:
# Command: Switch to Frame
"""
    frame!(frame::Element)::Nothing

The Switch To Frame command is used to select the current top-level browsing context or a child browsing context of the current browsing context to use as the current browsing context for subsequent commands.
"""
function frame! end
# function frame!(session::Session)::Nothing
#     @unpack addr, id = session
#     response = HTTP.post(
#         "$addr/session/$id/frame",
#         [("Content-Type" => "application/json")],
#         JSON3.write(Dict("id" => nothing))),
#     )
#     @assert response.status == 200
# end
function frame!(frame::Element)::Nothing
    @unpack addr, id = frame.session
    response = HTTP.post(
        "$addr/session/$id/frame",
        [("Content-Type" => "application/json")],
        JSON3.write(Dict("id" => ("ELEMENT" => frame.id))),
    )
    @assert response.status == 200
end


text:
"""
Encoding and decoding Modia models as JSON.

* Developer: Hilding Elmqvist, Mogram AB
* First version: August 2021
* License: MIT (expat)

"""
module JSONModel

export modelToJSON, JSONToModel, cloneModel, writeModel, readModel

import JSON
using Modia
using Unitful
using OrderedCollections
using ModiaBase.Symbolic: removeBlock
using Measurements
import MonteCarloMeasurements

"""
    jModel = encodeModel(model, expressionsAsStrings=true)

Encodes a model suitable to convert to a JSON string.

* `model`: model (declarations and equations)
* `expressionsAsStrings`: Julia expressions stored as strings or ast (Expr struct, etc) stored as dictionary
* `return transformed model`
"""
function encodeModel(model, expressionsAsStrings=true)
    Base.remove_linenums!(model)
    model = removeBlock(model)
    if typeof(model) <: Quantity
        return Dict(:_value => ustrip.(model), :_unit=>Modia.unitAsString(unit(model)))
    elseif expressionsAsStrings && typeof(model) == Expr
        return Dict(:_Expr=>string(model))
    elseif typeof(model) == Expr
        return Dict("head"=>model.head, "args"=>[encodeModel(a, expressionsAsStrings) for a in model.args])
    elseif typeof(model) == Symbol
        return Dict(:_Symbol=>string(model))
    elseif typeof(model) <: AbstractDict
        jModel = OrderedDict()
        for (n,v) in model
            jModel[n] = encodeModel(v, expressionsAsStrings)
        end
        return jModel
    elseif typeof(model) <: QuoteNode || typeof(model) <: Measurement || typeof(model) <: MonteCarloMeasurements.StaticParticles || typeof(model) <: MonteCarloMeasurements.Particles
        return Dict(:_type=>"$(typeof(model))", :_value => string(model))
    elseif typeof(model) <: AbstractArray
        return [encodeModel(m, expressionsAsStrings) for m in model]
    else
        return model
    end
end

"""
    json = modelToJSON(model; expressionsAsStrings=true)

Encodes a model as JSON. Expressions can be stored as strings or as ASTs.

* `model`: model (declarations and equations)
* `expressionsAsStrings`: Julia expressions stored as strings or Expr struct stored as dictionary
* `return JSON string for the model`
"""
function modelToJSON(model; expressionsAsStrings=true)
    jModel = encodeModel(model, expressionsAsStrings)
    return JSON.json(jModel, 2)
end

"""
    model = decodeModel(jModel)

Encodes a model suitable to convert to a JSON string.

* `jmodel`: encoded model
* `return model`
"""
function decodeModel(jModel)
    if typeof(jModel) <: AbstractDict && length(keys(jModel)) == 2 && "_value" in keys(jModel) && "_unit" in keys(jModel)
        unit = jModel["_unit"]
        #unit = replace(unit, ' ' => '*')  # Workaround for https://github.com/PainterQubits/Unitful.jl/issues/391  # Remove, since already correctly done in encodeModel(...)
        return jModel["_value"]*uparse(unit)
    elseif typeof(jModel) <: AbstractDict && length(keys(jModel)) == 1 && "_Expr" in keys(jModel)
        return Base.remove_linenums!( Meta.parse(jModel["_Expr"]) )  # without remove_linenums!, there are errors later with if-expressions.
    elseif typeof(jModel) <: AbstractDict && length(keys(jModel)) == 1 && "_Symbol" in keys(jModel)
        return Symbol(jModel["_Symbol"])
    elseif typeof(jModel) <: AbstractDict && length(keys(jModel)) == 2 && "head" in keys(jModel) && "args" in keys(jModel)
        ast = Expr(Symbol(jModel["head"]), decodeModel.(jModel["args"])...)
        return ast
    elseif typeof(jModel) <: AbstractDict && length(keys(jModel)) == 2 && "_type" in keys(jModel) && "_value" in keys(jModel)
        if jModel["_type"] == "QuoteNode"
            return QuoteNode(eval(Meta.parse(jModel["_value"])))
        else
            return eval(Meta.parse(jModel["_value"]))
        end
    elseif typeof(jModel) <: AbstractDict
        model = OrderedDict()
        for (n,v) in jModel
            model[Symbol(n)] = decodeModel(v)
        end
        return model
    elseif typeof(jModel) <: AbstractArray
        return [decodeModel(m) for m in jModel]
    else
        return jModel
    end
end

"""
    jModel = JSONToModel(json)

Decodes JSON string as model.

* `json`: JSON string for the model
* `return model`
"""
function JSONToModel(json)
    jModel = JSON.parse(json, dicttype=OrderedDict)
    model = decodeModel(jModel)
end

"""
    writeModel(filename::String, model; log=true)

Write model in JSON format on file `filename`.
"""
function writeModel(filename::String, model::AbstractDict; log=true)
    if log
        println("  Write model in JSON format on file \"$filename\"")
    end
    write(filename, modelToJSON(model))
end

"""
    model = readModel(filename::AbstractString; log=true)

Read a model that is stored in JSON format on file `filename` and return it.
"""
function readModel(filename::String; log=true)
    if log
        println("  Read model from JSON file \"$filename\"")
    end
    JSONToModel( read(filename,String) )
end

"""
    modelClone = cloneModel(model; expressionsAsStrings=true)

Clones a model by encoding/decoding as JSON.

* `model`: model
* `expressionsAsStrings`: Julia expressions stored as strings or Expr struct stored as dictionary
* `return modelClone`
"""
function cloneModel(model; expressionsAsStrings=true)
#    @showModel(model)
    jsonModel = modelToJSON(model, expressionsAsStrings=expressionsAsStrings)
#    @show jsonModel
    modelClone = JSONToModel(jsonModel)
#    @showModel(modelClone)
    return modelClone
end

# --------------------------------------------------------------------

#=
function test()

    model = Model(
        p = 3u"m/s",
        q = 4.5,
        r = :(p+q),
        v = Var(min=0),
        s = "asdf qwerty",
        t = 1 ± 0.1,
#        u = 1 ∓ 0.1,
        equations = :[
            2*der(x) = -p*x
            y = 2*x
        ]
    )

    model2 = Model(
#        t = 1 ± 0.1,
#        u = 1 ∓ 0.1,
    )

    println("Original model:")
    @showModel(model)
    println()
    @instantiateModel(model, log=true)

    json = modelToJSON(model)

    println("JSON for model:")
    println(json)

    println()

    modelClone = JSONToModel(json)
    @show modelClone
    println("Cloned model:")
    @showModel(modelClone)
    @instantiateModel(modelClone, log=true)

    println()

    json = modelToJSON(modelClone)
    println("JSON for cloned model:")
    println(json)

    # ------------------------

    Pin = Model( v = potential, i = flow )
    Ground = Model( p = Pin, equations = :[ p.v = 0.0u"V" ] )
    Gro = cloneModel(Ground)
    @showModel(Gro)
    @instantiateModel(Gro)
end

test()
=#

end


text:
using IcanProj
using NetCDF
using DataFrames
using PyPlot
using JFSM2
using Statistics

function compute_statistics(data)

    res = Dict("mean" => mean(data, dims = 2),
               "std" => std(data, dims = 2),
               "cv" => std(data, dims = 2) ./ mean(data, dims = 2))

    return res

end

function plot_hydrol_exchng(time, data, df_cfg, fig_title)

    hydrol_off = convert(Array{Bool}, df_cfg[:hydrol] .== 0)
    hydrol_on  = convert(Array{Bool}, df_cfg[:hydrol] .== 1)

    exchng_off = convert(Array{Bool}, df_cfg[:exchng] .== 0)
    exchng_on  = convert(Array{Bool}, df_cfg[:exchng] .== 1)

    figure(figsize = (7, 5))

    #figure()

    fill_between(time,
                 maximum(data[hydrol_off .& exchng_off,:], dims = 1)[:],
                 minimum(data[hydrol_off .& exchng_off,:], dims = 1)[:],
                 facecolor = "gray", edgecolor = "gray", alpha = 0.5,
                 label = "Hydrol=0 & Exchng=0")

    fill_between(time,
                 maximum(data[hydrol_on .& exchng_off,:], dims = 1)[:],
                 minimum(data[hydrol_on .& exchng_off,:], dims = 1)[:],
                 facecolor = "green", edgecolor = "green", alpha = 0.5,
                 label = "Hydrol=1 & Exchng=0")

    fill_between(time,
                 maximum(data[hydrol_off .& exchng_on,:], dims = 1)[:],
                 minimum(data[hydrol_off .& exchng_on,:], dims = 1)[:],
                 facecolor = "yellow", edgecolor = "yellow", alpha = 0.5,
                 label = "Hydrol=0 & Exchng=1")

    fill_between(time,
                 maximum(data[hydrol_on .& exchng_on,:], dims = 1)[:],
                 minimum(data[hydrol_on .& exchng_on,:], dims = 1)[:],
                 facecolor = "magenta", edgecolor = "magenta", alpha = 0.5,
                 label = "Hydrol=1 & Exchng=1")

    legend(loc=2)

    ylabel("SWE (mm)")

    title(fig_title)

end

println("USING OPEN DATA!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")

respath = "/data04/jmg/fsm_simulations/netcdf/fsmres_open"

figpath = joinpath(dirname(pathof(IcanProj)), "..", "plots", "config_ts")

variable = "swe"

iexp = 1:32

for spatial_res in ["1km", "50km"]

    time = load_time(respath, variable, spatial_res)

    data = load_time_slice(respath, variable, spatial_res, iexp)

    stats = compute_statistics(data)

    imin = argmin(dropdims(stats["cv"], dims = 2))

    imax = argmax(dropdims(stats["cv"], dims = 2))

    datamin = load_space_slice(respath, variable, spatial_res, iexp, imin)

    datamax = load_space_slice(respath, variable, spatial_res, iexp, imax)

    df_cfg = cfg_table()

    fig_title = "Spatial resolution $(spatial_res)"

    plot_hydrol_exchng(time, datamin, df_cfg, fig_title)

    savefig(joinpath(figpath, "datamin_$(spatial_res).png"), bbox_inches="tight", dpi=200)

    plot_hydrol_exchng(time, datamax, df_cfg, fig_title)

    savefig(joinpath(figpath, "datamax_$(spatial_res).png"), bbox_inches="tight", dpi=200)

end

#=

ipos = 62
figure()
for icfg in 1:32
    data = datamax[icfg, :]
    plot(collect(1:length(data)), data)
    annotate(string([df_cfg[icfg, i] for i in 1:6]), xy=(ipos, data[ipos]))
    @show data[ipos]
end

=#

#=
cfg_vec = [:albedo, :condct, :densty, :exchng, :hydrol]

for cfg in cfg_vec

    ioff = convert(Array{Bool}, df_cfg[cfg] .== 0)
    ion = convert(Array{Bool}, df_cfg[cfg] .== 1)

    figure()
    plot(time, datamin[ioff,:]', label = "off", color = "gray")
    plot(time, datamin[ion,:]', label = "on", color = "blue")
    ylabel("SWE (mm)")
    title(string(cfg))

    figure()
    plot(time, datamax[ioff,:]', label = "off", color = "gray")
    plot(time, datamax[ion,:]', label = "on", color = "blue")
    ylabel("SWE (mm)")
    title(string(cfg))

end

=#

#=

plot(time, datamin[hydrol_off .& exchng_on,:]', label = "on", color = "blue")
plot(time, datamin[hydrol_on .& exchng_off,:]', label = "off", color = "pink")
plot(time, datamin[hydrol_on .& exchng_on,:]', label = "on", color = "yellow")
ylabel("SWE (mm)")

figure()
plot(time, datamax[hydrol_off .& exchng_off,:]', label = "off", color = "gray")
plot(time, datamax[hydrol_off .& exchng_on,:]', label = "on", color = "blue")
plot(time, datamax[hydrol_on .& exchng_off,:]', label = "off", color = "pink")
plot(time, datamax[hydrol_on .& exchng_on,:]', label = "on", color = "yellow")
ylabel("SWE (mm)")

=#

#=
hydrol_off = convert(Array{Bool}, df_cfg[:hydrol] .== 0)
hydrol_on  = convert(Array{Bool}, df_cfg[:hydrol] .== 1)

exchng_off = convert(Array{Bool}, df_cfg[:exchng] .== 0)
exchng_on  = convert(Array{Bool}, df_cfg[:exchng] .== 1)

figure()
plot(time, datamin[hydrol_off .& exchng_off,:]', label = "off", color = "gray")
plot(time, datamin[hydrol_off .& exchng_on,:]', label = "on", color = "blue")
plot(time, datamin[hydrol_on .& exchng_off,:]', label = "off", color = "pink")
plot(time, datamin[hydrol_on .& exchng_on,:]', label = "on", color = "yellow")
ylabel("SWE (mm)")

figure()
plot(time, datamax[hydrol_off .& exchng_off,:]', label = "off", color = "gray")
plot(time, datamax[hydrol_off .& exchng_on,:]', label = "on", color = "blue")
plot(time, datamax[hydrol_on .& exchng_off,:]', label = "off", color = "pink")
plot(time, datamax[hydrol_on .& exchng_on,:]', label = "on", color = "yellow")
ylabel("SWE (mm)")
=#

#=

# Test for plotting maps

path = "/data02/Ican/vic_sim/fsm_past_1km/netcdf/tmp"

variable = "swe"

spatial_res = "50km"

iexp = 1:32

itime = 200

data = load_time_slice(path, variable, spatial_res, iexp)

stats = compute_statistics(data)

df = link_table(path, variable, spatial_res)

map_values = project_map(stats["mean"], df)

figure()
imshow(map_values)
cb = colorbar()
=#

#=
# Test for plotting time series

map_values = fill(NaN, 1550, 1195)

map_values[df[:ind_julia]] = df[:nc]

figure()
imshow(map_values)
cb = colorbar()

ispace = 169

spatial_res = "10km"

data = load_space_slice(path, variable, spatial_res, iexp, ispace)

figure()
plot(data')
=#


text:
using Test, TumorGrowth

@testset "Simulation" begin

    simoutput = birth_death_pushing(100; b=1., d=0.0, μ=0.1, ρ=1., dim=2, seed=1010, showprogress=false)
    @test simoutput isa Dict
    @test map((:index,:mutation,:tumor,:time)) do key
        haskey(simoutput,key)
    end |> all
    @test try
        birth_death_pushing(100; b=1., d=1., μ=0.5, ρ=1., dim=2, seed=1010, showprogress=false)
    catch e
        e == ErrorException("Tumor died")
    end

    birth_death_pushing!(simoutput[:tumor], simoutput[:mutations], length(simoutput[:tumor])+1; b=0.69, d=0.0, μ=0.3, dim=2,
     t=simoutput[:time], cur_id=simoutput[:index], cur_mutation=simoutput[:mutation])

     tumor = simoutput[:tumor] |> DataFrame
     @test map(eachcol(tumor)) do col
         eltype(eltype(col)) <: Number
     end |> all

     mutations = simoutput[:mutations] |> DataFrame
     function lambda()
         while true
             m = rand(1:length(mutations.origin))
             ori = mutations.origin[m]
             if ori in tumor.index
                 cell = findfirst(isequal(ori), tumor.index)
                 loc = findfirst(isequal(m), tumor.mutations[cell])

                 return if mutations.ancestor[m] == 0
                     loc == 1
                 else
                     mutations.ancestor[m] == tumor.mutations[loc-1]
                 end
             end
         end
     end
     @test lambda()

     TumorGrowth.b_linear()
     @test iszero(TumorGrowth.b_curve(1.1; bup=1., ρc=1.))
     TumorGrowth.b_hill(6)
     @test !iszero(TumorGrowth.b_curve(1.1; bup=1., ρc=1.))
     reduced_μ!(tumor, 1/2)

     tumor = birth_death_pushing(2000; b=1., d=0.0, μ=0.3, ρ=6., dim=3, seed=1234)[:tumor]
     TumorGrowth.CSV.write("temp.csv", DataFrame(tumor), delim=',')
     data_import("temp.csv", delim=",")
end

@testset "Sampling" begin
    tumor = data_import("temp.csv", delim=",")
    @test isa.( [
        cross_section(tumor; x=10., width=3.),
        radial_sample(tumor; r=30., width=3.),
        bulk(tumor; pos=(10.,0.,10.), box=(20,20,20)),
        punch(tumor; pos=(-10,-10,0), r=5)
        ], DataFrame) |> all
    @test isa.([
        tumor |> mutation_freqs,
        tumor |> stochastic_sequencing
        ], DataFrame) |> all

    plane = cross_section(tumor; x=0., width=3., reduce_dim=true)
    lattice, samples, sample_r = multi_region_sampling(plane; n = 20, cells_per_sample = 20)
    @test lattice isa Vector{Vector{Float64}} &&
        samples isa Vector{DataFrame}
    samples, sampletumor = multi_region_sequencing(plane; n = 20, cells_per_sample = 20, stochastic = true)
    @test samples isa Vector{DataFrame} && sampletumor isa DataFrame
end

@testset "Plotting" begin
    tumor = data_import("temp.csv", delim=",")

    sc = plotting(tumor, shading=false, inline=false, size=(500,500))
    plotting_colored_mutations!(sc, tumor, colorpalette = TumorGrowth.palette(:tab20), show_warning = true, shading=false, path = "temp.png")

    plane = cross_section(tumor; x=0., width=3., reduce_dim=true)
    samples, sampletumor = multi_region_sequencing(plane; n = 20)
    plotting_colored_mutations(sampletumor; markersize=sampletumor.sample_r, show_axis=true)
    plotting_sampletumor_pies(sampletumor, first=10, colorpalette = TumorGrowth.palette(:tab10))
    plotting_sampletumor_pies(sampletumor, mutations = [1,10,55])

    @test sc isa TumorGrowth.Scene
    display(sc)
end

@testset "Timeseries" begin
    N, simparams... = (N=1000, b=1., d=0.0, μ=0.3, ρ=6., dim=3, seed=1234)
    tumor = birth_death_pushing(N; simparams...)[:tumor] |> DataFrame

    time_series = tumor_stepper(range(0., last( tumor.t_birth), length=50); simparams...)

    @test last(time_series) == tumor

    record_growth(time_series; path="temp.gif",
            frames=1, shading=true,
            points_colors = t-> colors_by_mutations(t; colorpalette = TumorGrowth.palette(:tab20), show_warning=false)
            )
end
rm("temp.png")
rm("temp.csv")
rm("temp.gif")

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
# deactivate plot GUI, which is not available in Docker
ENV["GKSwstype"] = "100"

# instantiate project
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

const TARGET_FOLDER = "result"
const RESULTS_FILE = "results.csv"

function main()
    if !isdir(TARGET_FOLDER)
        mkdir(TARGET_FOLDER)
    end
    global io = open(joinpath(TARGET_FOLDER, RESULTS_FILE), "w")

    println("Running NLN benchmarks...")

    # Robertson benchmark (ROBE21)
    println("###\nRunning Robertson benchmark\n###")
    include("models/Robertson/robertson_benchmark.jl")

    # Coupled Van der Pol benchmark (CVP20)
    println("###\nRunning Van der Pol benchmark\n###")
    include("models/VanDerPol/vanderpol_benchmark.jl")

    # Laub-Loomis benchmark (LALO20)
    println("###\nRunning Laub-Loomis benchmark\n###")
    include("models/LaubLoomis/laubloomis_benchmark.jl")

    # Lotka-Volterra tangential crossing benchmark (LOVO21)
    println("###\nRunning Lotka-Volterra tangential crossing benchmark\n###")
    include("models/LotkaVolterra/lotka_volterra_benchmark.jl")

    # Spacecraft benchmark (SPRE21)
    println("###\nSpacecraft benchmark\n###")
    include("models/Spacecraft/spacecraft_benchmark.jl")

    print(io, "\n")
    println("Finished running benchmarks.")
    close(io)
    nothing
end

main()