text:
using TopOpt, LinearAlgebra, StatsFuns
using TopOptMakie
using VTKDataIO
using PyCall
using VTKDataTypes
using DashVtk
using Dash
using DashHtmlComponents

#dash = pyimport("dash")
#dash_html = pyimport("dash_html_components")
dash_vtk = pyimport("dash_vtk")

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 4.0 # filter radius
xmin = 0.0001 # minimum density
problem_size = (60, 20)
V = 0.5 # maximum volume fraction
p = 4.0 # penalty
x0 = fill(V, prod(problem_size)) # initial design

problem = HalfMBB(Val{:Linear}, problem_size, (1.0, 1.0), E, v, f)

solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
cheqfilter = DensityFilter(solver, rmin = rmin)
comp = TopOpt.Compliance(problem, solver)

function obj(x)
    # minimize compliance
    return comp(cheqfilter(x))
end
function constr(x)
    # volume fraction constraint
    return sum(cheqfilter(x)) / length(x) - V
end

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr)

options = MMAOptions(
    maxiter=1000, tol = Tolerance(kkt = 1e-4, x = 1e-4, f = 1e-4),
)
TopOpt.setpenalty!(solver, p)
# Method of Moving Asymptotes
@time r = Nonconvex.optimize(m, MMA87(), x0, options = options);

@show obj(r.minimizer)
@show constr(r.minimizer)
topology = cheqfilter(r.minimizer)

#function GeometryBasics.Mesh(problem, topology)
#    mesh = VTKUnstructuredData(problem)
#    topology = round.(topology)
#    inds = findall(isequal(0), topology)
#    deleteat!(mesh.cell_connectivity, inds)
#    deleteat!(mesh.cell_types, inds)
#    return GeometryBasics.Mesh(mesh)
#end

mesh = VTKUnstructuredData(problem)
topology = round.(topology)
inds = findall(isequal(0), topology)
deleteat!(mesh.cell_connectivity, inds)
deleteat!(mesh.cell_types, inds)

pymesh = PyVTK(mesh)
dash_vtk_utils = pyimport("dash_vtk.utils")
dash_vtk_utils.to_mesh_state
state = dash_vtk_utils.to_mesh_state(pymesh)

content = vtk_view([
        vtk_geometryrepresentation([
        vtk_mesh(state=state)
    ]),
])

app = dash()

app.layout = html_div() do
    html_div(
    style = Dict(
        "width" => "100%",
        "height" => "400px"
    ),
    children=[content],
    )
end

run_server(app, "0.0.0.0", debug = true)


text:
"""
Abstract Experimental Data types
"""
abstract type ExperimentalData end
"""
Base for data in the form of samples
"""
abstract type SampleData <: ExperimentalData end
"""
Base for data in the form of a distribution
"""
abstract type HistogramData <: ExperimentalData end

abstract type AbstractRNAData{hType} <: HistogramData end

"""
Data structures

Do not use underscore "_" in name

arguments:
name: label for the data set
gene: gene name (case sensitive)
nRNA: length of histogram
histRNA: RNA histograms
bins: number of live cell recording time bins
OFF:  OFF time probability density
ON:: ON time probability density
"""
struct RNAData{nType,hType} <: AbstractRNAData{hType}
    name::String
    gene::String
    nRNA::nType
    histRNA::hType
end
struct TransientRNAData{nType,tType,hType} <: AbstractRNAData{hType}
    name::String
    gene::String
    nRNA::nType
    time::tType
    histRNA::hType
end
struct RNAMixedData{hType} <: AbstractRNAData{hType}
    name::String
    gene::String
    nRNA::Array
    fish::Array{Bool,1}
    histRNA::hType
end
struct LiveCellData <: HistogramData
    name::String
    gene::String
    bins::Array
    OFF::Array
    ON::Array
end
struct RNALiveCellData <: HistogramData
    name::String
    gene::String
    nRNA::Int
    histRNA::Array
    bins::Array
    OFF::Array
    ON::Array

end
struct MultiRNALiveCellData <: HistogramData
    name::String
    gene::String
    nRNA::Array
    histRNA::Tuple
    bins::Tuple
    OFF::Tuple
    ON::Tuple
end

"""
Abstract model types
"""
abstract type Model end
abstract type StochasticGRmodel <: Model end
abstract type AbstractGMmodel <: StochasticGRmodel end
abstract type AbstractGRMmodel <: StochasticGRmodel end
abstract type AbstractGMlossmodel <: AbstractGMmodel end

"""
Model structures

arguments:
G: number of G steps
R: number of R steps
nalleles: number of alleles producing RNA
type: type of model (e.g. splice, premature RNA death, etc.)
rates: transition rates
rateprior: prior for rates
proposal: MCMC proposal distribution
fittedparam: indices of rates to be fitted
    randomeffects: indices of rates that are fixed to each other, in the form of a 2tuple of vectors
    with index 1 the tied index vector and 2 the corresponding fitted index vector
method: numerical method for solving Master equation
"""
struct GMmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GMrescaledmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GMmultimodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GMtransientmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GMdelaymodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GMlossmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMlossmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GMfixedeffectslossmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMlossmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    fixedeffects::Tuple
    method::MethodType
end

struct GMmixedmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGMmodel
    G::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GRMmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGRMmodel
    G::Int
    R::Int
    nalleles::Int
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

struct GRSMmodel{RateType,PriorType,ProposalType,ParamType,MethodType} <: AbstractGRMmodel
    G::Int
    R::Int
    nalleles::Int
    type::String
    rates::RateType
    rateprior::PriorType
    proposal::ProposalType
    fittedparam::ParamType
    method::MethodType
end

function write_model(model)

end

"""
Abstract Option types for fitting methods
"""
abstract type Options end
abstract type Results end

# Model and Data dependent functions

"""
datahistogram(data)
Return the RNA histogram data as one vector
"""
# datahistogram(data::RNAData) = data.histRNA
function datahistogram(data::AbstractRNAData{Array{Array,1}})
# function datahistogram(data::TransientRNAData)
    v = data.histRNA[1]
    for i in 2:length(data.histRNA)
        v = vcat(v,data.histRNA[i])
    end
    return v
end
datahistogram(data::AbstractRNAData{Array{Float64,1}}) = data.histRNA
datahistogram(data::RNALiveCellData) = [data.OFF;data.ON;data.histRNA]

datapdf(data::AbstractRNAData{Array{Float64,1}}) = normalize_histogram(data.histRNA)
datapdf(data::RNALiveCellData) = [normalize_histogram(data.OFF);normalize_histogram(data.ON);normalize_histogram(data.histRNA)]
# datapdf(data::AbstractRNAData{Array{Array,1}}) = normalize_histogram(datahistogram(data))
function datapdf(data::AbstractRNAData{Array{Array,1}})
# function datahistogram(data::TransientRNAData)
    v = normalize_histogram(data.histRNA[1])
    for i in 2:length(data.histRNA)
        v = vcat(v,normalize_histogram(data.histRNA[i]))
    end
    return v
end

#### Model likelihoods   ####

"""
likelihoodfn(param,data,model)
model likelihoodfn called by loglikelihood
"""
function likelihoodfn(param,data::RNAData,model::GMmodel)
    r = get_rates(param,model)
    n = model.G-1
    steady_state(r[1:2*n+2],n,data.nRNA,model.nalleles)
end
function likelihoodfn(param,data::RNAData,model::AbstractGMlossmodel)
    r = get_rates(param,model)
    yieldfactor = r[end]
    n = model.G-1
    steady_state(r[1:2*n+2],yieldfactor,n,data.nRNA,model.nalleles)
end
function likelihoodfn(param,data::RNAData{T1,T2},model::AbstractGMlossmodel) where {T1 <: Array, T2 <: Array}
    h = likelihoodarray(get_rates(param,model),data,model)
    hconcat = Array{Float64,1}(undef,0)
    for h in h
        hconcat = vcat(hconcat,h)
    end
    return hconcat
end
function likelihoodfn(param::Vector,data::RNALiveCellData,model::GRSMmodel)
    modelOFF, modelON, histF = likelihoodtuple(get_rates(param,model),data,model)
    return [modelOFF;modelON;histF]
end

"""
likelihoodarray(r,data,model)

Compute likelihoods for multiple distributions

r = full rate vector including yield if it exists
data = data structure
model = model structure

For time dependent model likelihoods
first set of parameters gives the initial histogram
2nd set gives the new parameters at subsequent times
data.histRNA holds array of histograms for time points given by data.time
transient computes the time evolution of the histogram
model.method=1 specifies finite difference solution otherwise use eigendecomposition solution,
"""
function likelihoodarray(r,data::TransientRNAData,model::AbstractGMmodel)
    h=likelihoodarray(r,data,model,maximum(data.nRNA))
    trim(h,data.nRNA)
end
function likelihoodarray(r,data::TransientRNAData,model::AbstractGMlossmodel)
    yieldfactor = r[end]
    nh = nhist_loss(maximum(data.nRNA),yieldfactor)
    h = likelihoodarray(r[1:end-1],data::TransientRNAData,model,yieldfactor,nh)
    # technical_loss!(h,yieldfactor)
    trim(h,data.nRNA)
end
function likelihoodarray(r,data::RNAData{T1,T2},model::AbstractGMlossmodel) where {T1 <: Array, T2 <: Array}
    yieldfactor = r[end]
    n = model.G-1
    h = Array{Array{Float64,1},1}(undef,length(data.nRNA))
    for i in eachindex(data.nRNA)
        h[i] =steady_state(r[(i-1)*2*model.G+1 : i*2*model.G],yieldfactor,n,data.nRNA[i],model.nalleles)
    end
    trim(h,data.nRNA)
end
function likelihoodarray(r,data::TransientRNAData,model::AbstractGMlossmodel,yieldfactor,maxdata)
    G = model.G
    h0 = steady_state_full(r[1:2*G],G-1,maxdata)
    transient(data.time,r[2*G+1:4*G],yieldfactor,G-1,model.nalleles,h0,model.method)
    # transient(t,r,yieldfactor,n,nalleles,P0::Vector,method)
end
function likelihoodarray(r,data::RNAData,model::GMmultimodel)
    G = model.G
    h = Array{Array{Float64,1},1}(undef,length(data.nRNA))
    for i in eachindex(data.nRNA)
        g = steady_state(r[1:2*G],G-1,data.nRNA[i],model.nalleles)
        h[i] = threshold_noise(g,r[2*G+1],r[2*G+1+i],data.nRNA[i])
    end
    return h
end
function likelihoodarray(r,data::RNAMixedData,model::AbstractGMmodel)
    G = model.G
    h = Array{Array{Float64,1},1}(undef,length(data.nRNA))
    j = 1
    for i in eachindex(data.fish)
        g = steady_state(r[1:2*G],G-1,maximum(data.nRNA),model.nalleles)
        if data.fish[i]
            h[i] = threshold_noise(g,r[2*G+j],r[2*G+j+1],data.nRNA[i])
            j += 2
        else
            h[i] = technical_loss(g,r[2*G+j],data.nRNA[i])
            j += 1
        end
    end
    return h
end

function likelihoodtuple(r,data::RNALiveCellData,model::GRSMmodel)
    if model.type == "offdecay"
        modelOFF, modelON = offonPDF(data.bins,r,model.G-1,model.R,model.method)
        histF = steady_state_offpath(r,model.G-1,model.R,data.nRNA,model.nalleles)
    elseif model.type == "offeject"
        modelOFF, modelON = offonPDF_offeject(data.bins,r,model.G-1,model.R,model.method)
        histF = steady_state_offeject(r,model.G-1,model.R,data.nRNA,model.nalleles)
    else
        modelOFF, modelON = offonPDF(data.bins,r,model.G-1,model.R,model.method)
        histF = steady_state(r,model.G-1,model.R,data.nRNA,model.nalleles)
    end
    return modelOFF, modelON, histF
end

"""
get_rates(param,model)
replace fitted rates with new values and return
"""
function get_rates(param,model::AbstractGMmodel)
    r = get_r(model)
    r[model.fittedparam] = param
    return r
end
function get_rates(param,model::GRSMmodel)
    r = get_r(model)
    r[model.fittedparam] = param
    setr(r,model)
end
get_rates(param,model::GMtransientmodel) = param[1:2*model.G]

"""
get_rates(param,model::GMrescaledmodel)

gammas are scaled by nu
"""
function get_rates(param,model::GMrescaledmodel)
    r = get_r(model)
    n = get_n(model)
    nu = n in model.fittedparam ? param[findfirst(model.fittedparam .== n)] : r[n]
    r[1:n-1] /= r[n]
    r[model.fittedparam] = param
    r[1:n-1] *= nu
    if r[2*model.G + 3] > 1
        r[2*model.G + 3] = 1
    end
    return r
end

function get_rates(param,model::GMfixedeffectslossmodel)
    r = get_r(model)
    n = get_n(model)
    r[model.fittedparam] = param
    r[model.fixedeffects[1]] = r[model.fixedeffects[2]]
    return r
end

get_r(model) = copy(model.rates)
get_n(model) = 2*model.G - 1

"""
get_param(model)

get fitted parameters from model
"""
get_param(model::StochasticGRmodel) = model.rates[model.fittedparam]

function get_param(model::GMrescaledmodel)
    r = copy(model.rates)
    n = 2*model.G - 1
    r[1:n-1] /= r[n]
    r[model.fittedparam]
end

"""
setr(r,model)

"""
function setr(r,model::GRSMmodel)
    n = model.G-1
    nr = model.R
    eta = geteta(r,n,nr)
    r[2*n + 1 + nr + 1:2*n + 1 + nr + nr] = eta
    r
end
"""
logprior(param,model::AbstractGMmodel)

compute log of the prior
"""
function logprior(param,model::AbstractGMmodel)
    d = model.rateprior
    p=0
    for i in eachindex(d)
        p -= logpdf(d[i],param[i])
    end
    return p
end

"""
logprior(x,model::GRSModel)
Compute log prior using distributions in Model.rateprior
called by mhstep() in metropolis_hastings.jl
"""
function logprior(param,model::GRSMmodel)
    r = get_rates(param,model)
    d = model.rateprior
    G = model.G
    R = model.R
    p=0
    j = 1
    #  G rates
    for i in Grange(G)
        p -= logpdf(d[j],r[i])
        j += 1
    end
    # initiation rate
    i = initiation(G)
    p -= logpdf(d[j],r[i])
    j += 1
    # sum of R Steps rates are bounded by length of insertion site to end of gene, i.e sum of reciprocal rates is bounded
    t = sum(1 ./ r[Rrange(G,R)])
    p -= logpdf(d[j],1/t)
    j += 1
    # priors for splice rates
    rs = 0
    for i in Srange(G,R)
        rs += r[i]
        p -= logpdf(d[j],rs)
        j += 1
    end
    return p
end

ext: jl

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
module CTPlots

Base.Experimental.@optlevel 1

export plot_image, plot_sinogram, plot_tomogram
export plot_gray_scale

using RecipesBase: plot, plot!
using RecipesBase
using ..Monads
using ..CTImages: CTImage, CTSinogram, CTTomogram, ctfn, CTImageOrTomog
using ..CTImages: ctimage, ctsinogram, cttomogram
using ..TestImages: AbstractTestImage, ImageParams, gray_scale_indices
using ..CTScan: AbstractCTScanner
using ..Utils: linspace, ORI
using IntervalSets

function __init__()
    ENV["GKS_ENCODING"] = "utf-8"
end

"""
    plot_image(image::AbstractMatrix{T}; <keyword arguments>) where {T}

Plot CT image `image` with predefined options.

Additional keyword arguments are passed to the Plots library functions.

See also: [`plot_sinogram`](@ref), [`plot_tomogram`](@ref)
"""
function plot_image(image::AbstractMatrix{T}; kwargs...) where {T}
    plot(CTImage(image))
end

"""
    plot_tomogram(tomog::AbstractMatrix{T}; <keyword arguments>) where {T}

Plot reconstructed CT image `tomog` with predefined options.

Alias for `plot_image`. See [`plot_image`](@ref) for full reference.
Additional keyword arguments are passed to the Plots library.

See also: [`plot_sinogram`](@ref), [`plot_image`](@ref)
"""
function plot_tomogram(tomog::AbstractMatrix{T}; kwargs...) where {T}
    plot(CTTomogram(tomog); kwargs...)
end

"""
    plot_sinogram(sinog::AbstractMatrix{T}; <keyword arguments>) where {T}

Plot sinogram `sinog` with predefined options.

Additional keyword arguments are passed to the underlying plotting driver.

See also: [`plot_image`](@ref), [`plot_tomogram`](@ref)
"""
function plot_sinogram(sinog::AbstractMatrix{T}; kwargs...) where {T}
    plot(CTSinogram(sinog); kwargs...)
end

# TODO: Check lims for plot again!

@recipe function plot(
    xs::AbstractVector, ys::AbstractVector, image::CTImageOrTomog
)
    rows, cols = size(image)
    proj = get(plotattributes, :projection, nothing)
    aspect_ratio --> :equal
    #seriestype --> :heatmap
    _seriestype = get(plotattributes, :seriestype, nothing)
    if proj === :polar
        yaxis --> false
    elseif _seriestype === :heatmap
        tick_direction --> :out
        # xlims --> ((xs[1], xs[end]) .+ (-0.5, 0.5))
        xlims --> (xs[1], xs[end])
        # ylims --> ((ys[1], ys[end]) .+ (-0.5, 0.5))
        ylims --> (ys[1], ys[end])
    end
    xs, ys, mjoin(image)
end

@recipe function plot(
    xsi::Interval, ysi::Interval, image::CTImageOrTomog
)
    rows, cols = size(image)
    xs = linspace(xsi, cols)
    ys = linspace(ysi, rows)
    xs, ys, image
end

@recipe function plot(
    image::CTImageOrTomog,
    α::Optional{Real} = nothing,
    β::Optional{Real} = nothing,
)
    rows, cols = size(image)
    proj = get(plotattributes, :projection, nothing)
    if proj === :polar
        α = maybe(rows, α)
        β = maybe(2π, β)
        rs = linspace(0..α, rows+1)
        ϕs = linspace(0..β, cols+1)
        ϕs, rs, image
    else
        β = maybe(maybe((rows - 1) / 2, α), β)
        α = maybe((cols - 1) / 2, α)
        xs = linspace(-α..α, cols)
        ys = linspace(-β..β, rows)
        xs, ys, image
    end
end

const _sinog_xticks = [45i for i in 0:8]

@recipe function plot(xs::AbstractVector, ys::AbstractVector, sinog::CTSinogram)
    #seriestype --> :heatmap
    _seriestype = get(plotattributes, :seriestype, nothing)
    if _seriestype === :heatmap
        tick_direction --> :out
        # xlims --> ((xs[begin], xs[end]) .+ (-0.5, 0.5) .* Δϕ)
        xlims --> (xs[begin], xs[end])
        # ylims --> ((ys[begin], ys[end]) .+ (-0.5, 0.5))
        ylims --> (ys[begin], ys[end])
    end
    xs, ys, mjoin(sinog)
end

@recipe function plot(θsi::Interval, tsi::Interval, sinog::CTSinogram)
    nd, nϕ = size(sinog)
    ts = linspace(tsi, nd)
    θs = linspace(θsi, nϕ)
    θs, ts, sinog
end

@recipe function plot(sinog::CTSinogram, α::Optional{Real} = nothing)
    nd, nϕ = size(sinog)
    xs = linspace(ORI(0..360), nϕ) # Now we use midpoints and not edges to support more backends
    α = maybe((nd - 1) / 2, α)
    ys = linspace(-α..α, nd)
    xticks --> _sinog_xticks
    xrotation --> -45
    xs, ys, sinog
end

@recipe function plot(
    xs::Union{AbstractVector,Interval},
    ys::Union{AbstractVector,Interval},
    gs::AbstractTestImage
)
    xs, ys, gs.image
end

@recipe function plot(
    gs::AbstractTestImage,
    α::Optional{Real} = nothing,
    β::Optional{Real} = nothing
)
    gs.image, α, β
end

@recipe function plot(imp::ImageParams, gray_scale_data::AbstractVector...)
    _, xs = gray_scale_indices(imp)
    seriestype --> :scatter
    linewidth --> 1.3
    markersize --> 1.8
    markerstrokewidth --> 0
    legend --> :outertopright
    x --> (xs .- imp.cols ÷ 2)
    gray_scale_data
end

@recipe function plot(gst::AbstractCTScanner, s::Symbol)
    gst, Val(s)
end

for (nm, ctf) ∈ pairs(ctfn)
    v = Val{nm}
    @eval begin
        @recipe function plot(gst::AbstractCTScanner, ::$v)
            mjoin($ctf(gst))
        end
    end
end

"""
    plot_image(gst::AbstractCTScanner)

Plot input image created for the test `gst`.
"""
plot_image(gst::AbstractCTScanner; kwargs...) = plot(gst, :image; kwargs...)

"""
    plot_sinogram(gst::AbstractCTScanner)

Plot the sinogram computed for the test `gst`.
"""
plot_sinogram(gst::AbstractCTScanner; kwargs...) = plot(gst, :sinog; kwargs...)

"""
    plot_tomogram(gst::AbstractCTScanner)

Plot the reconstructed image for the test `gst`.
"""
plot_tomogram(gst::AbstractCTScanner; kwargs...) = plot(gst, :tomog; kwargs...)

function _get_options(options, i)
    if isnothing(options)
        return (;)
    end
    if length(options) == 0
        return (;)
    end
    if i > length(options)
        return (;)
    end
    options[i]
end

"""
    plot_gray_scale(
        imp::ImageParams,
        gray_scale_data...;
        options = nothing,
        <keyword arguments>
    )

Plot gray scale data.

Additional keyword arguments are passed to the Plots library.
"""
function plot_gray_scale(
    imp::ImageParams,
    data::AbstractVector...;
    options = nothing,
    kwargs...)
    p = plot(imp, first(data); st = :line)
    for i ∈ 2:length(data)
        opts = _get_options(options, i)
        plot!(p, imp, data[i]; opts..., kwargs...)
    end
    p
end

end # module


text:
function ontologygraphs(ontol)
    Dict("P" => ontologygraph(ontol, biological_process),
         "F" => ontologygraph(ontol, molecular_function),
         "C" => ontologygraph(ontol, cellular_component))
end

function ontologygraph(ontol, ontology=biological_process)
    gids = collect(Iterators.filter(k -> (ontol[k].namespace == ontology) && (ontol[k].id == k), keys(ontol))) ## select from namespace and confirm its not an alt id
    n = length(gids)
    graphindex = Dict(gids[i] => i for i = 1:n)

    gograph = SimpleDiGraph(n)
    for (id, gt) in ontol
        (gt.namespace != ontology) && continue
        (gt.id != id) && continue

        gi = graphindex[id]

        for j in gt.is_a
            add_edge!(gograph, gi, graphindex[j])
        end
        for j in gt.part_of
            add_edge!(gograph, gi, graphindex[j])
        end
        for j in gt.positively_regulates
            add_edge!(gograph, gi, graphindex[j])
        end
        for j in gt.negatively_regulates
            add_edge!(gograph, gi, graphindex[j])
        end
        for j in gt.regulates
            add_edge!(gograph, gi, graphindex[j])
        end
    end
    gographclosed = transitiveclosure(gograph, true)
    (gids=gids, graphindex=graphindex, gograph=gograph, gographclosed=gographclosed)
end

function extendannotation(annot, gographs)

    exannot = Dict{String, Dict{String, Vector{Int}}}()
    for (ont, ann) in annot
        exannot[ont] = Dict{String, Vector{Int}}()
        for (gid, gis) in ann
            try
                exannot[ont][gid] = goancestors(gographs[ont], gis)
            catch
                @show gid
                @show gis
                @show ont

                error("")
            end
        end
    end

    exannot
end

goancestors(gg, gis::AbstractVector)  = mapreduce(gi -> goancestors(gg, gi), vcat, gis) |> unique
goancestors(gg, gi) = gg.gids[neighbors(gg.gographclosed, gg.graphindex[gi])]

function test_ex(gene, exannot, annot)
    a = annot["P"][gene]
    e = exannot["P"][gene]
    @show setdiff(a, e)
    all([x ∈ e for x in a])

end

#test_ex("gene5966|vegt", exannot, annot)


text:
using NetCDF

# Add array with time in netcdfs with forcing data

path_st = "/data04/jmg/fsm_simulations/netcdf/forcings_st"

files = readdir(path_st)

for file in files

    file = joinpath(path_st, file)

    dim_time = ncread(file, "dim_time")

    time_str = ncread(file, "time_str")

    time = DateTime.(time_str, "yyyy-mm-dd HH:MM:SS")

    time_array = [Dates.year.(time) Dates.month.(time) Dates.day.(time) Dates.hour.(time)]

    nccreate(file, "time_array", "ymdh", 4, Dict("ymdh" => "header"), "dim_time")

    ncclose()

    ncwrite(time_array', file, "time_array")

    ncclose()

    @show file

end


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 2, 16, :x, Quadrilateral(4, 4), false),
  :properties => [10.92 0.3;],
  :x_coords => collect(range(0.0, stop=0.5, length=3)),
  :y_coords => collect(range(0.0, stop=0.5, length=3)),
  :thickness => 1.0,
  :support => [
    (1, [0 0 0 1]),
    (2, [0 0 1 1]),
    (3, [0 0 1 0]),
    (4, [0 1 0 1]),
    (6, [1 0 1 0]),
    (7, [0 1 0 0]),
    (8, [1 1 0 0]),
    (9, [1 0 0 0])
    ],
  :loaded_nodes => [
    (9, [0.25 0.0 0.0 0.0])
    ]
)

@time fm_df, sigma_df = p47(data)

@test round(sigma_df[4,3], digits=4) ≈ -0.0335 atol = 0.001


text:
function lininterp1(x,y,xi)

    x = x[:]

    if isnothing(findfirst(xi.<x))
        placeLow = length(x)-1
    else
        placeLow = findfirst(xi.<x)-1
    end

    if placeLow == 0
        placeLow = 1
    end

    placeHigh = placeLow+1
    xLow = x[placeLow]
    xHigh = x[placeHigh]
    yLow = y[placeLow]
    yHigh = y[placeHigh]

    yi = yLow +(xi-xLow)*(yHigh-yLow)/(xHigh-xLow)

end


text:
# a handler for Proj4
immutable Proj4Handler   <: AbstractPackageHandler; end     # for point handled by the Proj4 package

"""
Points with a full coordinate reference system as defined by an SRID identifier

Its up to the user to determine the what the x / y / z fields actually represent; which is governed by the element order in Proj4

For a quick reference:
    lat long style CRS's,  x -> lon, y -> lat (or get_lat() and get_lon())
    utm style CRS's,  x -> false east, y -> false north, z -> up (or get_east() get_north() get_up())
"""
immutable CRS{T <: AbstractSRID} <: WorldPosition
    x::Float64
    y::Float64
    z::Float64
end

# useful shortcuts
typealias CRS_NULL CRS{UnknownSRID}

default_params{T <: CRS}(::Type{T}) = (UnknownSRID,)

# trait style functions
has_srid{T <: CRS}(::Type{T}) = Val{true}
get_handler{T <: CRS}(::Type{T}) = Proj4Handler

# get the datum
get_datum{T}(::Union{Type{CRS{T}}, CRS{T}}) = get_datum(T)

######################################################################################
### Experimental, compound coordinate reference systems (CCRS)                     ###
######################################################################################

"""
Abstract type for compound coordinate reference system (i.e. height is not ellipsoidal)
"""
abstract AbstractCCRS{T, U} <: WorldPosition

# use the SRID style because we need to Proj4 to handle the Geoid anyway
"""
Compound coordinate reference system where the height is relative to a geoid
"""
immutable CCRS_Geoid{T <: AbstractSRID, G <: AbstractGeoid} <: AbstractCCRS{T, G}
    x::Float64
    y::Float64
    z::Float64
end

# convenient typecasts for compound coordinate reference systems
typealias CCRS_NULL  CCRS_Geoid{UnknownSRID, UnknownGeoid}

default_params{T <: CCRS_Geoid}(::Type{T}) = (UnknownSRID,
                                              UnknownGeoid)

# trait style functions
has_srid{T <: CCRS_Geoid}(::Type{T}) = Val{true}
has_geoid{T <: CCRS_Geoid}(::Type{T}) = Val{true}
get_handler{T <: CCRS_Geoid}(::Type{T}) = Proj4Handler

# get the horizontal datum
get_datum{T,G}(::Union{Type{CCRS_Geoid{T, G}}, CCRS_Geoid{T, G}}) = get_datum(T)

# Use Geodesy to initialize a standard set of methods
Proj4_types = [CRS, CCRS_Geoid]
for (i, t) in enumerate(Proj4_types)
     Geodesy.eval(build_methods(t, fieldnames(t), Geodesy.GeodesyTypes))
end


text:
using PatentsLens
using Test

# [i > length(cites(f)) for (i,f) in zip(indegree(g), fams)] |> any
# [o > length(citedby(f)) for (o,f) in zip(outdegree(g), fams)] |> any


text:
"""Return information in an experiment `e` in column(s) `feature` after filtering with filter `f`.
"""
function diagnostic(e::Experiment,
                    f::AbstractFilter;
                    features = [:Metadata_Row, :Metadata_Column, :Metadata_Field])
    subdata = e.data[filter_entries(e, f),:]
    return(subdata[:,features])
end

"""Get path to images in an experiment `e` stored in feature `s` after filtering with filter `f`.
If `center` is true, the path will be returned together with the center of the filtered entries,
assuming they are stored in variables `AreaShape_Center_X`, `AreaShape_Center_X_1`,
`AreaShape_Center_Y` and `AreaShape_Center_Y_1`.
If `rgx` provides a list of regex substitutions, it will be applied on all file paths in the
output (which is useful if you're in a different file system or environment as the one described
in the `Experiment`'s data).
"""
function diagnostic_path(e::Experiment,
                         f::AbstractFilter,
                         s::Symbol;
                         center = false,
                         rgx = nothing)
    diagCols = [s]
    if center
        centerCols = [:AreaShape_Center_X, :AreaShape_Center_Y,
                      :AreaShape_Center_X_1, :AreaShape_Center_Y_1]
        append!(diagCols, centerCols)
    end
    diagCellCoord = diagnostic(e, f, features = diagCols)

    diagURL = diagCellCoord[:,s]
    unique!(diagURL)
    if rgx !== nothing
        diagURL = map(x -> reduce(replace, rgx, init=x), diagURL)
    end

    if !center
        # Center of regions of interest are not required
        return(diagURL)
    else
        diagURLtoCenters = Dict(x => Array{Int64,2}(undef,0,2) for x in diagURL)
        for (url,cx,cy,nx,ny) in eachrow(diagCellCoord)
            if rgx !== nothing
                url = reduce(replace, rgx, init=url)
            end
            diagURLtoCenters[url] = vcat(diagURLtoCenters[url],
                                         [[round(cx),round(nx)] [round(cy),round(ny)]])
        end
        return((diagURL, map(x -> diagURLtoCenters[x], diagURL)))
    end
end

"""[intended for internal use only]
Return color image (handle merging of 3 channels if needed)
"""
function colimgifrgb(imgPath, rgb)
	if rgb === nothing
		# Images don't have to be merged (e.g. already RGB)
    	return(load(imgPath))
    else
    	# Use the transformations provided in `rgb` to get a
    	# color image
    	imgR = reduce(replace, rgb[1], init=imgPath)
    	imgG = reduce(replace, rgb[2], init=imgPath)
    	imgB = reduce(replace, rgb[3], init=imgPath)
    	return(getColorImage(imgR, imgG, imgB))
    end
end

"""[intended for internal use only]
Assemble output file name for a given path and folder structure
"""
function outimgpath(imgpath, basepath, keepsubfolders)
	name_ending = join(split(imgpath, "/")[end-keepsubfolders:end], "/")
	name = string(basepath, name_ending)
	name_folder = join(split(name, "/")[1:end-1], "/")
	mkpath(name_folder)
	return(string(basepath, name_ending))
end

"""Get images in an experiment `e` whose location is stored in feature `s` after filtering with filter `f`.
Images will be saved at the `path` provided if `saveimages` is set to true,
while keeping up to `keepsubfolders` folders from the original folder structure of the images.
If `center` is true, crosses will indicate the center of selected objects,
assuming they are stored in variables `:AreaShape_Center_X`, `:AreaShape_Center_X_1`,
`:AreaShape_Center_Y` and `:AreaShape_Center_Y_1` (default if you merge nuclear and cytoplasm
measurements from CellProfiler outputs).
Display up to `showlimit` images if `show` is true.
Save up to `savelimit` images if `saveimages` is true.
If `rgx` provides a list of regex substitutions, it will be applied on all image paths
(which is useful if you're in a different file system or environment as the one described
in the `Experiment`'s data).
If `rgb` provides a list of 3 lists of regex substitutions, it will be applied to generate
the path to 3 images.
"""
function diagnostic_images(e::Experiment,
                          f::AbstractFilter,
                          s::Symbol;
                          path = "./",
                          saveimages = true,
                          show = false,
                          center = false,
                          showlimit::Int64 = 20,
                          savelimit::Int64 = 20,
                          rgb = nothing,
                          rgx = nothing,
                          keepsubfolders = 0)

    # Get addresses of images matching criteria
    imagesURL = diagnostic_path(e, f, s; center = center, rgx = rgx)

    # A limit is set to the number of images displayed by default
    # To avoid overflowing notebooks
    showcounter = 0
    savecounter = 0
    if !center
        for imgPath = imagesURL
        	colImg = colimgifrgb(imgPath, rgb)
            if show & (showcounter < showlimit)
                showcounter += 1
                display(colImg)
            end
            if saveimages & (savecounter < savelimit)
                savecounter += 1
            	save(outimgpath(imgPath, path, keepsubfolders), colImg)
            elseif showcounter >= showlimit | !show
            	# Images are not displayed nor saved
            	# so we can safely exit
            	break
            end
        end
    else
        # We highlight the center of region of interests provided
        for (imgPath, imgCenters) = zip(imagesURL[1], imagesURL[2])
            # Combine them into a single image, display if required and save it
            colImg = colimgifrgb(imgPath, rgb)

            # Add white crosses around each center
            for (x,y) = eachrow(imgCenters)
                x += 1
                y += 1
                view(colImg, max(y-10,1):min(y+10,size(colImg)[1]), x) .= RGB(1,1,1)
                view(colImg, y, max(x-10,1):min(x+10,size(colImg)[2])) .= RGB(1,1,1)
            end
            if show & (showcounter < showlimit)
                showcounter += 1
                display(colImg)
            end
            if saveimages & (savecounter < savelimit)
                savecounter += 1
            	save(outimgpath(imgPath, path, keepsubfolders), colImg)
            elseif showcounter >= showlimit | !show
            	# Images are not displayed nor saved
            	# so we can safely exit
            	break
            end
        end
    end
    return true
end

"""[intended for internal use only]
Return an (intensity-normalized) color image, given `R`, `G` and `B` the paths to
3 single-channel images."""
function getColorImage(R::String, G::String, B::String; normalize = true)
    imgR = load(R);
    imgG = load(G);
    imgB = load(B);

    if normalize
        imgR ./= maximum(imgR)
        imgG ./= maximum(imgG)
        imgB ./= maximum(imgB)
    end
    colorview(RGB, imgR, imgG, imgB)
end

# Expand freqtable to support Experiment objects,
# either to know how many entries are selected
# by a filter or the values taken by a feature
# for the subset of the entries selected.
"""
Expand `freqtable` to support `Experiment` objects.
Find the frequency of the values taken by feature `s`
in Experiment `e`.
"""
function FreqTables.freqtable(e::Experiment,
                              s::Symbol;
                              args...)
    return(freqtable(getdata(e)[!,s]; args...))
end

"""
Expand `freqtable` to support `Experiment` objects.
Find the frequency of the values taken by feature `s`
in Experiment `e`.
"""
function FreqTables.freqtable(e::Experiment,
                              f::AbstractFilter;
                              args...)
    entries_kept = [x in filter_entries(e, f) ? "Kept" : "Discarded"
                    for x in e.selected_entries]
    return(freqtable(entries_kept; args...))
end


text:
"""
    importCRS(x::GeoFormatTypes.GeoFormat; [order=:compliant])

Import a coordinate reference system from a `GeoFormat` into GDAL,
returning an [`AbstractSpatialRef`](@ref).

## Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, the default will use axis
  ordering compliant with the relevant CRS authority.
"""
importCRS(x::GFT.GeoFormat; kwargs...) =
    importCRS!(newspatialref(; kwargs...), x)

unsafe_importCRS(x::GFT.GeoFormat; kwargs...) =
    importCRS!(unsafe_newspatialref(; kwargs...), x)

"""
    importCRS!(spref::AbstractSpatialRef, x::GeoFormatTypes.GeoFormat)

Import a coordinate reference system from a `GeoFormat` into the spatial ref.
"""
function importCRS! end

importCRS!(spref::AbstractSpatialRef, x::GFT.EPSG) =
    importEPSG!(spref, GFT.val(x))
importCRS!(spref::AbstractSpatialRef, x::GFT.AbstractWellKnownText) =
    importWKT!(spref, GFT.val(x))
importCRS!(spref::AbstractSpatialRef, x::GFT.ESRIWellKnownText) =
    importESRI!(spref, GFT.val(x))
importCRS!(spref::AbstractSpatialRef, x::GFT.ProjString) =
    importPROJ4!(spref, GFT.val(x))
importCRS!(spref::AbstractSpatialRef, x::GFT.GML) =
    importXML!(spref, GFT.val(x))
importCRS!(spref::AbstractSpatialRef, x::GFT.KML) =
    importCRS!(spref, GFT.EPSG(4326))

"""
    reproject(points, sourceproj::GeoFormat, destproj::GeoFormat; [order=:compliant])

Reproject points to a different coordinate reference system and/or format.

## Arguments
- `coord`: Vector of Geometry points
- `sourcecrs`: The current coordinate reference system, as a `GeoFormat`
- `targetcrs`: The coordinate reference system to transform to, using any CRS capable `GeoFormat`

## Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, the default will use axis
  ordering compliant with the relevant CRS authority.

## Example
```julia-repl
julia> using ArchGDAL, GeoFormatTypes

julia> ArchGDAL.reproject([[118, 34], [119, 35]], ProjString("+proj=longlat +datum=WGS84 +no_defs"), EPSG(2025))
2-element Array{Array{Float64,1},1}:
 [-2.60813482878655e6, 1.5770429674905164e7]
 [-2.663928675953517e6, 1.56208905951487e7]
```
"""
reproject(coord, sourcecrs::GFT.GeoFormat, targetcrs::Nothing; kwargs...) = coord

# These should be better integrated with geometry packages or follow some standard
const ReprojectCoord = Union{<:NTuple{2,<:Number},<:NTuple{3,<:Number},AbstractVector{<:Number}}

# Vector/Tuple coordinate(s)
reproject(coord::ReprojectCoord, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat; kwargs...) =
    GeoInterface.coordinates(reproject(createpoint(coord...), sourcecrs, targetcrs; kwargs...))
reproject(coords::AbstractArray{<:ReprojectCoord}, sourcecrs::GFT.GeoFormat,
          targetcrs::GFT.GeoFormat; kwargs...) =
    GeoInterface.coordinates.(reproject([createpoint(c...) for c in coords], sourcecrs, targetcrs; kwargs...))

# GeoFormat
reproject(geom::GFT.GeoFormat, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat; kwargs...) =
    convert(typeof(geom), reproject(convert(Geometry, geom), sourcecrs, targetcrs; kwargs...))

# Geometries
function reproject(geom::AbstractGeometry, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat; kwargs...)
    crs2transform(sourcecrs, targetcrs; kwargs...) do transform
        transform!(geom, transform)
    end
end
function reproject(geoms::AbstractArray{<:AbstractGeometry}, sourcecrs::GFT.GeoFormat,
                   targetcrs::GFT.GeoFormat; kwargs...)
    crs2transform(sourcecrs, targetcrs; kwargs...) do transform
        transform!.(geoms, Ref(transform))
    end
end

"""
    crs2transform(f::Function, sourcecrs::GeoFormat, targetcrs::GeoFormat; kwargs...)

Run the function `f` on a coord transform generated from the source and target
crs definitions. These can be any `GeoFormat` (from GeoFormatTypes) that holds
a coordinate reference system.

`kwargs` are passed through to `importCRS`.
"""
function crs2transform(f::Function, sourcecrs::GFT.GeoFormat, targetcrs::GFT.GeoFormat; kwargs...)
    importCRS(sourcecrs; kwargs...) do sourcecrs_ref
        importCRS(targetcrs; kwargs...) do targetcrs_ref
            createcoordtrans(sourcecrs_ref, targetcrs_ref) do transform
                f(transform)
            end
        end
    end
end

"""
    newspatialref(wkt::AbstractString = ""; order=:compliant)

Construct a Spatial Reference System from its WKT.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
newspatialref(wkt::AbstractString = ""; order=:compliant) =
    maybesetaxisorder!(ISpatialRef(GDAL.osrnewspatialreference(wkt)), order)

unsafe_newspatialref(wkt::AbstractString = ""; order=:compliant) =
    maybesetaxisorder!(SpatialRef(GDAL.osrnewspatialreference(wkt)), order)

function maybesetaxisorder!(sr::AbstractSpatialRef, order)
    if order == :trad
        GDAL.osrsetaxismappingstrategy(sr.ptr, GDAL.OAMS_TRADITIONAL_GIS_ORDER)
    elseif order != :compliant
        throw(ArgumentError("order $order is not supported. Use :trad or :compliant"))
    end
    sr
end

function destroy(spref::AbstractSpatialRef)
    GDAL.osrdestroyspatialreference(spref.ptr)
    spref.ptr = C_NULL
end

"""
    clone(spref::AbstractSpatialRef)

Makes a clone of the Spatial Reference System. May return NULL.
"""
function clone(spref::AbstractSpatialRef)
    if spref.ptr == C_NULL
        return ISpatialRef()
    else
        return ISpatialRef(GDAL.osrclone(spref.ptr))
    end
end

function unsafe_clone(spref::AbstractSpatialRef)
    if spref.ptr == C_NULL
        return SpatialRef()
    else
        return SpatialRef(GDAL.osrclone(spref.ptr))
    end
end

"""
    importEPSG!(spref::AbstractSpatialRef, code::Integer)

Initialize SRS based on EPSG GCS or PCS code.

This method will initialize the spatial reference based on the passed in
EPSG GCS or PCS code. It is relatively expensive, and generally involves quite a
bit of text file scanning. Reasonable efforts should be made to avoid calling it
many times for the same coordinate system.

### Additional Remarks
This method is similar to importFromEPSGA() except that EPSG preferred axis
ordering will not be applied for geographic coordinate systems. EPSG normally
defines geographic coordinate systems to use lat/long contrary to typical GIS
use). Since OGR 1.10.0, EPSG preferred axis ordering will also not be applied
for projected coordinate systems that use northing/easting order.

The coordinate system definitions are normally read from the EPSG derived
support files such as pcs.csv, gcs.csv, pcs.override.csv, gcs.override.csv and
falling back to search for a PROJ.4 epsg init file or a definition in epsg.wkt.

These support files are normally searched for in /usr/local/share/gdal or in the
directory identified by the GDAL_DATA configuration option. See CPLFindFile()
for details.
"""
function importEPSG!(spref::AbstractSpatialRef, code::Integer)
    result = GDAL.osrimportfromepsg(spref.ptr, code)
    @ogrerr result "Failed to initialize SRS based on EPSG"
    return spref
end

"""
    importEPSG(code::Integer; [order=:compliant])

Construct a Spatial Reference System from its EPSG GCS or PCS code.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
importEPSG(code::Integer; kwargs...) =
    importEPSG!(newspatialref(; kwargs...), code)

unsafe_importEPSG(code::Integer; kwargs...) =
    importEPSG!(unsafe_newspatialref(; kwargs...), code)

"""
    importEPSGA!(spref::AbstractSpatialRef, code::Integer)

Initialize SRS based on EPSG CRS code.

This method is similar to `importFromEPSG()` except that EPSG preferred axis
ordering will be applied for geographic and projected coordinate systems. EPSG
normally defines geographic coordinate systems to use lat/long, and also there
are also a few projected coordinate systems that use northing/easting order
contrary to typical GIS use). See `importFromEPSG()` for more
details on operation of this method.
"""
function importEPSGA!(spref::AbstractSpatialRef, code::Integer)
    result = GDAL.osrimportfromepsga(spref.ptr, code)
    @ogrerr result "Failed to initializ SRS based on EPSGA"
    return spref
end

"""
    importEPSGA(code::Integer; [order=:compliant])

Construct a Spatial Reference System from its EPSG CRS code.

This method is similar to `importFromEPSG()` except that EPSG preferred axis
ordering will be applied for geographic and projected coordinate systems. EPSG
normally defines geographic coordinate systems to use lat/long, and also there
are also a few projected coordinate systems that use northing/easting order
contrary to typical GIS use). See `importFromEPSG()` for more
details on operation of this method.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
importEPSGA(code::Integer; kwargs...) =
    importEPSGA!(newspatialref(; kwargs...), code)

unsafe_importEPSGA(code::Integer; kwargs...) =
    importEPSGA!(unsafe_newspatialref(; kwargs...), code)

"""
    importWKT!(spref::AbstractSpatialRef, wktstr::AbstractString)

Import from WKT string.

This method will wipe the existing SRS definition, and reassign it based on the
contents of the passed WKT string. Only as much of the input string as needed to
construct this SRS is consumed from the input string, and the input string
pointer is then updated to point to the remaining (unused) input.
"""
function importWKT!(spref::AbstractSpatialRef, wktstr::AbstractString)
    result = GDAL.osrimportfromwkt(spref.ptr, [wktstr])
    @ogrerr result "Failed to initialize SRS based on WKT string"
    return spref
end

"""
    importWKT(wktstr::AbstractString; [order=:compliant])

Create SRS from its WKT string.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
importWKT(wktstr::AbstractString; kwargs...) =
    newspatialref(wktstr; kwargs...)

unsafe_importWKT(wktstr::AbstractString; kwargs...) =
    unsafe_newspatialref(wktstr; kwargs...)

"""
    importPROJ4!(spref::AbstractSpatialRef, projstr::AbstractString)

Import PROJ.4 coordinate string.

The OGRSpatialReference is initialized from the passed PROJ.4 style coordinate
system string. In addition to many `+proj` formulations which have OGC
equivalents, it is also possible to import `"+init=epsg:n"` style definitions.
These are passed to `importFromEPSG()`. Other init strings (such as the state
plane zones) are not currently supported.

Example: `pszProj4 = \"+proj=utm +zone=11 +datum=WGS84\"`

Some parameters, such as grids, recognized by PROJ.4 may not be well understood
and translated into the OGRSpatialReference model. It is possible to add the
`+wktext` parameter which is a special keyword that OGR recognized as meaning
\"embed the entire PROJ.4 string in the WKT and use it literally when converting
back to PROJ.4 format\".

For example: `\"+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150
+ellps=intl +units=m +nadgrids=nzgd2kgrid0005.gsb +wktext\"`
"""
function importPROJ4!(spref::AbstractSpatialRef, projstr::AbstractString)
    result = GDAL.osrimportfromproj4(spref.ptr, projstr)
    @ogrerr result "Failed to initialize SRS based on PROJ4 string"
    return spref
end

"""
    importPROJ4(projstr::AbstractString; [order=:compliant])

Create SRS from its PROJ.4 string.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
importPROJ4(projstr::AbstractString; kwargs...) =
    importPROJ4!(newspatialref(; kwargs...), projstr)

unsafe_importPROJ4(projstr::AbstractString; kwargs...) =
    importPROJ4!(unsafe_newspatialref(; kwargs...), projstr)

"""
    importESRI!(spref::AbstractSpatialRef, esristr::AbstractString)

Import coordinate system from ESRI .prj format(s).

This function will read the text loaded from an ESRI .prj file, and translate it
into an OGRSpatialReference definition. This should support many (but by no
means all) old style (Arc/Info 7.x) .prj files, as well as the newer pseudo-OGC
WKT .prj files. Note that new style .prj files are in OGC WKT format, but
require some manipulation to correct datum names, and units on some projection
parameters. This is addressed within `importFromESRI()` by an automatic call to
`morphFromESRI()`.

Currently only `GEOGRAPHIC`, `UTM`, `STATEPLANE`, `GREATBRITIAN_GRID`, `ALBERS`,
`EQUIDISTANT_CONIC`, `TRANSVERSE (mercator)`, `POLAR`, `MERCATOR` and
`POLYCONIC` projections are supported from old style files.

At this time there is no equivalent `exportToESRI()` method. Writing old style
.prj files is not supported by OGRSpatialReference. However the `morphToESRI()`
and `exportToWkt()` methods can be used to generate output suitable to write to
new style (Arc 8) .prj files.
"""
function importESRI!(spref::AbstractSpatialRef, esristr::AbstractString)
    result = GDAL.osrimportfromesri(spref.ptr, [esristr])
    @ogrerr result "Failed to initialize SRS based on ESRI string"
    return spref
end

"""
    importESRI(esristr::AbstractString; kwargs...)

Create SRS from its ESRI .prj format(s).

Passing the keyword argument `order=:compliant` or `order=:trad` will set the
mapping strategy to return compliant axis order or traditional lon/lat order.
"""
importESRI(esristr::AbstractString; kwargs...) =
    importESRI!(newspatialref(; kwargs...), esristr)

unsafe_importESRI(esristr::AbstractString; kwargs...) =
    importESRI!(unsafe_newspatialref(; kwargs...), esristr)

"""
    importXML!(spref::AbstractSpatialRef, xmlstr::AbstractString)

Import SRS from XML format (GML only currently).
"""
function importXML!(spref::AbstractSpatialRef, xmlstr::AbstractString)
    result = GDAL.osrimportfromxml(spref.ptr, xmlstr)
    @ogrerr result "Failed to initialize SRS based on XML string"
    return spref
end

"""
    importXML(xmlstr::AbstractString; [order=:compliant])

Construct SRS from XML format (GML only currently).

Passing the keyword argument `order=:compliant` or `order=:trad` will set the
mapping strategy to return compliant axis order or traditional lon/lat order.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
importXML(xmlstr::AbstractString; kwargs...) =
    importXML!(newspatialref(; kwargs...), xmlstr)

unsafe_importXML(xmlstr::AbstractString; kwargs...) =
    importXML!(unsafe_newspatialref(; kwargs...), xmlstr)

"""
    importURL!(spref::AbstractSpatialRef, url::AbstractString)

Set spatial reference from a URL.

This method will download the spatial reference at a given URL and feed it into
SetFromUserInput for you.
"""
function importURL!(spref::AbstractSpatialRef, url::AbstractString)
    result = GDAL.osrimportfromurl(spref.ptr, url)
    @ogrerr result "Failed to initialize SRS from URL"
    return spref
end

"""
    importURL(url::AbstractString; [order=:compliant])

Construct SRS from a URL.

This method will download the spatial reference at a given URL and feed it into
SetFromUserInput for you.

# Keyword Arguments
- `order`: Sets the axis mapping strategy. `:trad` will use traditional lon/lat axis
  ordering in any actions done with the crs. `:compliant`, will use axis
  ordering compliant with the relevant CRS authority.
"""
importURL(url::AbstractString; kwargs...) =
    importURL!(newspatialref(; kwargs...), url)

unsafe_importURL(url::AbstractString; kwargs...) =
    importURL!(unsafe_newspatialref(; kwargs...), url)

"""
    toWKT(spref::AbstractSpatialRef)

Convert this SRS into WKT format.
"""
function toWKT(spref::AbstractSpatialRef)
    wktptr = Ref{Cstring}()
    result = GDAL.osrexporttowkt(spref.ptr, wktptr)
    @ogrerr result "Failed to convert this SRS into WKT format"
    return unsafe_string(wktptr[])
end

"""
    toWKT(spref::AbstractSpatialRef, simplify::Bool)

Convert this SRS into a nicely formatted WKT string for display to a person.

### Parameters
* `spref`:      the SRS to be converted
* `simplify`:   `true` if the `AXIS`, `AUTHORITY` and `EXTENSION` nodes should be
                stripped off.
"""
function toWKT(spref::AbstractSpatialRef, simplify::Bool)
    wktptr = Ref{Cstring}()
    result = GDAL.osrexporttoprettywkt(spref.ptr, wktptr, simplify)
    @ogrerr result "Failed to convert this SRS into pretty WKT"
    return unsafe_string(wktptr[])
end

"""
    toPROJ4(spref::AbstractSpatialRef)

Export coordinate system in PROJ.4 format.
"""
function toPROJ4(spref::AbstractSpatialRef)
    projptr = Ref{Cstring}()
    result = GDAL.osrexporttoproj4(spref.ptr, projptr)
    @ogrerr result "Failed to export this SRS to PROJ.4 format"
    return unsafe_string(projptr[])
end

"""
    toXML(spref::AbstractSpatialRef)

Export coordinate system in XML format.

Converts the loaded coordinate reference system into XML format to the extent
possible. LOCAL_CS coordinate systems are not translatable. An empty string will
be returned along with OGRERR_NONE.
"""
function toXML(spref::AbstractSpatialRef)
    xmlptr = Ref{Cstring}()
    result = GDAL.osrexporttoxml(spref.ptr, xmlptr, C_NULL)
    @ogrerr result "Failed to convert this SRS into XML"
    return unsafe_string(xmlptr[])
end

"""
    toMICoordSys(spref::AbstractSpatialRef)

Export coordinate system in Mapinfo style CoordSys format.
"""
function toMICoordSys(spref::AbstractSpatialRef)
    ptr = Ref{Cstring}()
    result = GDAL.osrexporttomicoordsys(spref.ptr, ptr)
    @ogrerr result "Failed to convert this SRS into XML"
    return unsafe_string(ptr[])
end

"""
    morphtoESRI!(spref::AbstractSpatialRef)

Convert in place to ESRI WKT format.

The value nodes of this coordinate system are modified in various manners more
closely map onto the ESRI concept of WKT format. This includes renaming a
variety of projections and arguments, and stripping out nodes note recognised by
ESRI (like AUTHORITY and AXIS).
"""
function morphtoESRI!(spref::AbstractSpatialRef)
    result = GDAL.osrmorphtoesri(spref.ptr)
    @ogrerr result "Failed to convert in place to ESRI WKT format"
    return spref
end

"""
    morphfromESRI!(spref::AbstractSpatialRef)

Convert in place from ESRI WKT format.

The value notes of this coordinate system are modified in various manners to
adhere more closely to the WKT standard. This mostly involves translating a
variety of ESRI names for projections, arguments and datums to "standard" names,
as defined by Adam Gawne-Cain's reference translation of EPSG to WKT for the CT
specification.

Missing parameters in `TOWGS84`, `DATUM` or `GEOGCS` nodes can be added to the
`WKT`, comparing existing `WKT` parameters to GDAL's databases. Note that this
optional procedure is very conservative and should not introduce false
information into the WKT definition (although caution should be advised when
activating it). Needs the Configuration Option `GDAL_FIX_ESRI_WKT` be set to one
of the following (`TOWGS84` recommended for proper datum shift calculations)

`GDAL_FIX_ESRI_WKT` values:

* `TOWGS84` Adds missing TOWGS84 parameters (necessary for datum
        transformations), based on named datum and spheroid values.
* `DATUM`   Adds EPSG AUTHORITY nodes and sets SPHEROID name to OGR spec.
* `GEOGCS`  Adds EPSG AUTHORITY nodes and sets `GEOGCS`, `DATUM` and `SPHEROID`
        names to OGR spec. Effectively replaces `GEOGCS` node with the result of
        `importFromEPSG(n)`, using `EPSG` code `n` corresponding to the existing
        `GEOGCS`. Does not impact `PROJCS` values.
"""
function morphfromESRI!(spref::AbstractSpatialRef)
    result = GDAL.osrmorphfromesri(spref.ptr)
    @ogrerr result "Failed to convert in place from ESRI WKT format"
    return spref
end

"""
    setattrvalue!(spref::AbstractSpatialRef, path::AbstractString, value::AbstractString)

Set attribute value in spatial reference.

Missing intermediate nodes in the path will be created if not already in
existence. If the attribute has no children one will be created and assigned
the value otherwise the zeroth child will be assigned the value.

### Parameters
* `path`: full path to attribute to be set. For instance "PROJCS|GEOGCS|UNIT".
* `value`: (optional) to be assigned to node, such as "meter". This may be left
            out if you just want to force creation of the intermediate path.
"""
function setattrvalue!(
        spref::AbstractSpatialRef,
        path::AbstractString,
        value::AbstractString
    )
    result = GDAL.osrsetattrvalue(spref.ptr, path, value)
    @ogrerr result "Failed to set attribute path to value"
    return spref
end

function setattrvalue!(spref::AbstractSpatialRef, path::AbstractString)
    result = GDAL.osrsetattrvalue(spref.ptr, path, C_NULL)
    @ogrerr result "Failed to set attribute path"
    return spref
end

"""
    getattrvalue(spref::AbstractSpatialRef, name::AbstractString, i::Integer)

Fetch indicated attribute of named node.

This method uses GetAttrNode() to find the named node, and then extracts the
value of the indicated child. Thus a call to `getattrvalue(spref,"UNIT",1)`
would return the second child of the UNIT node, which is normally the length of
the linear unit in meters.

Parameters
`name` the tree node to look for (case insensitive).
`i`    the child of the node to fetch (zero based).

Returns
the requested value, or NULL if it fails for any reason.
"""
getattrvalue(spref::AbstractSpatialRef, name::AbstractString, i::Integer) =
    GDAL.osrgetattrvalue(spref.ptr, name, i)

"""
    unsafe_createcoordtrans(source::AbstractSpatialRef, target::AbstractSpatialRef)

Create transformation object.

### Parameters
* `source`: source spatial reference system.
* `target`: target spatial reference system.

### Returns
NULL on failure or a ready to use transformation object.
"""
function unsafe_createcoordtrans(
        source::AbstractSpatialRef,
        target::AbstractSpatialRef
    )
    return CoordTransform(GDAL.octnewcoordinatetransformation(source.ptr,
        target.ptr))
end

"OGRCoordinateTransformation destructor."
function destroy(obj::CoordTransform)
    GDAL.octdestroycoordinatetransformation(obj.ptr)
    obj.ptr = C_NULL
end

"""
    transform!(xvertices, yvertices, zvertices, obj::CoordTransform)

Transform points from source to destination space.

### Parameters
* `xvertices`   array of nCount X vertices, modified in place.
* `yvertices`   array of nCount Y vertices, modified in place.
* `zvertices`   array of nCount Z vertices, modified in place.

### Returns
`true` on success, or `false` if some or all points fail to transform.
"""
function transform!(
        xvertices::Vector{Cdouble},
        yvertices::Vector{Cdouble},
        zvertices::Vector{Cdouble},
        obj::CoordTransform
    )
    # The method TransformEx() allows extended success information to be captured
    # indicating which points failed to transform.
    n = length(xvertices)
    @assert length(yvertices) == n
    @assert length(zvertices) == n
    return Bool(GDAL.octtransform(obj.ptr, n, pointer(xvertices),
        pointer(yvertices), pointer(zvertices)))
end


text:
eval_at(args...) = f -> f(args...)
findallin(xs, ys) = map(x -> findfirst(y -> x == y, ys), xs)

"""
    catalan_number(n)

calculate the n-th catalan number
"""
catalan_number(n) = div(binomial(2n, n), n+1)

"""
    number_binary_trees(n)

calculate the number of binary trees with n leafs
"""
number_binary_trees(n) = catalan_number(n - 1)

"""
    categorical_sample(tokens, weights)

sample one token according to the given weights
"""
function categorical_sample(tokens, weights)
    T = eltype(weights)
    x = convert(T, rand()) * sum(weights)
    cum_weights = zero(T)
    for (t, w) in zip(tokens, weights)
        cum_weights += w
        if cum_weights >= x
            return t
        end
    end
    error("categorical sample failed with length(tokens)=$(length(tokens)), weights=$weights")
end

categorical_sample(d::Dict)   = categorical_sample(keys(d), values(d))
categorical_sample(v::Vector) = categorical_sample(1:length(v), v)


text:
module WMATA
include("authentication.jl")
include("rail.jl")

import DataFrames.DataFrame, JSON.parse, HTTP.request, DataFrames.rename!, DataFrames.names, DataFrames.nrow
export WMATA_auth
export get_rail_predictions, get_station_list, get_station_timings, get_path_between, get_station_to_station
export get_rail_incidents, get_train_positions
export wmata_request

end


text:
using InteractiveUtils

function print_tree(t, indent)
  println("$(repeat("    ", indent))$t")
  for s in subtypes(t)
    if s != Any && s != Function
      print_tree(s, indent + 1)
    end
  end
end

print_tree(Any, 0)

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using DataStructures
using Base.Test

# construction

@test typeof(OrderedDict()) == OrderedDict{Any,Any}
@test typeof(OrderedDict('a',1)) == OrderedDict{Char,Int}
@test typeof(OrderedDict([("a",1),("b",2)])) == OrderedDict{ASCIIString,Int}

# empty dictionary
d = OrderedDict(Char, Int)
@test length(d) == 0
@test isempty(d)
@test_throws KeyError d['c'] == 1
d['c'] = 1
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification

for c in 'a':'z'
    d[c] = c-'a'+1
end

@test (d['a'] += 1) == 2
@test 'a' in keys(d)
@test haskey(d, 'a')
@test get(d, 'B', 0) == 0
@test !('B' in keys(d))
@test !haskey(d, 'B')
@test pop!(d, 'a') == 2

@test collect(keys(d)) == ['b':'z']
@test collect(values(d)) == [2:26]
@test collect(d) == [(a,i) for (a,i) in zip('b':'z', 2:26)]

# Test for #60

od60 = OrderedDict{Int,Int}()
od60[1] = 2

ranges = [2:5,6:9,10:13]
for range in ranges
    for i = range
        od60[i] = i+1
    end
    for i = range
        delete!( od60, i )
    end
end
od60[14]=15

@test od60[14] == 15


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash(suppress_callback_exceptions = true)

app.layout = html_div() do
    html_button("Add Filter", id="add-filter", n_clicks=0),
    html_div([], id="container")
end

callback!(app,
        Output("container", "children"),
        Input("add-filter", "n_clicks"),
        State("container", "children")
    ) do n_clicks, children
    new_element = html_div() do
        dcc_dropdown(
            id = (type = "dropdown", index = n_clicks),
            options = [(label = i, value = i) for i in ["NYC", "MTL", "LA", "TOKYO"]]
        ),
        html_div(
            id = (type = "output", index = n_clicks)
        )
    end
    return vcat(children, new_element)
end

callback!(app,
    Output((type = "output", index = MATCH), "children"),
    Input((type = "dropdown", index = MATCH), "value"),
    State((type = "dropdown", index = MATCH), "id"),
) do value, id
    return html_div("Dropdown $(id.index) = $(value)")
end

run_server(app)


text:
using BinDeps
using Compat
import Compat: String

@BinDeps.setup

function validate_mbed(name, handle)
    try
        get_version = Libdl.dlsym(handle, :mbedtls_version_get_string)
        version_ptr = Vector{UInt8}(9)
        ccall(get_version, Void, (Ptr{Void},), version_ptr)
        version = VersionNumber(unsafe_string(pointer(version_ptr)))
        version >= v"2.1.1"
    catch err
        warn("Could not check MbedTLS version: $err")
        false
    end
end

mbed = library_dependency("libmbedtls", aliases=["libmbedtls", "libmbedtls.2.1.1"])
mbed_crypto = library_dependency("libmbedcrypto", aliases=["libmbedcrypto", "libmbedcrypto.2.1.1"], validate=validate_mbed)
mbed_x509 = library_dependency("libmbedx509", aliases=["libmbedx509", "libmbedx509.2.1.1"])

mbed_all = [mbed, mbed_crypto, mbed_x509]

if haskey(ENV, "USE_GPL_MBEDTLS")  # The source code is identical except for the license text
    source_uri = URI("https://cache.julialang.org/https://tls.mbed.org/download/mbedtls-2.1.1-gpl.tgz")
    srcsha = "22c76e9d8036a76e01906423b3e8a02ab0ef84027f791bd719fff8edee9c61a9"
else
    source_uri = URI("https://cache.julialang.org/https://tls.mbed.org/download/mbedtls-2.1.1-apache.tgz")
    srcsha = "8f25b6f156ae5081e91bcc58b02455926d9324035fe5f7028a6bb5bc0139a757"
end

provides(Sources,
        source_uri,
        mbed_all, unpacked_dir="mbedtls-2.1.1",
        SHA = srcsha)

if is_unix()
    mbed_dir = joinpath(BinDeps.depsdir(mbed), "src", "mbedtls-2.1.1")
    provides(BuildProcess,
        (@build_steps begin
            `./cmake_check.sh`
            GetSources(mbed)
            @build_steps begin
                ChangeDirectory(mbed_dir)
                 @build_steps begin
                    `cmake -DUSE_SHARED_MBEDTLS_LIBRARY=On .`
                    `make lib`
                end
            end
        end), mbed_all, installed_libpath=joinpath(mbed_dir, "library"))
end

if is_apple()
    if Pkg.installed("Homebrew") === nothing
		error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
	end
	using Homebrew
	provides(Homebrew.HB, "mbedtls", mbed_all)
end

if is_windows()
    unpacked_dir = Int==Int32 ? "usr/bin32" : "usr/bin64"
    provides(
        Binaries,
        URI("https://malmaud.github.io/files/mbedtls-2.1.1-r1.zip"),
        mbed_all,
        unpacked_dir=unpacked_dir,
        SHA = "ab5a86d6c35d478082722e08747742fe04bf761a8e3ac4f3c960159244bbd8d8")
end

@BinDeps.install Dict("libmbedtls"=>"MBED_TLS",
                      "libmbedcrypto"=>"MBED_CRYPTO",
                      "libmbedx509"=>"MBED_X509")


text:
import Celeste.ArgumentParse

@testset "argument parse" begin

@testset "positional argument" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "my_arg")

    @test_throws ArgumentParse.ArgumentParsingError ArgumentParse.parse_args(parser, String[])

    parsed_args = ArgumentParse.parse_args(parser, ["hello"])
    @test parsed_args["my_arg"] == "hello"

    @test_throws ArgumentParse.ArgumentParsingError ArgumentParse.parse_args(parser, ["foo", "bar"])
end

@testset "positional argument type conversion" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "my_arg", arg_type=Int)
    parsed_args = ArgumentParse.parse_args(parser, ["123"])
    @test parsed_args["my_arg"] == 123
end

@testset "keyword argument" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "--my-arg")

    parsed_args = ArgumentParse.parse_args(parser, String[])
    @test !haskey(parsed_args, "my-arg")

    parsed_args = ArgumentParse.parse_args(parser, ["--my-arg", "hello"])
    @test parsed_args["my-arg"] == "hello"
end

@testset "keyword argument type conversion" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "--my-arg", arg_type=Int)

    parsed_args = ArgumentParse.parse_args(parser, ["--my-arg", "123"])
    @test parsed_args["my-arg"] == 123
end

@testset "required keyword argument fails when missing" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "--my-arg", required=true)
    @test_throws ArgumentParse.ArgumentParsingError ArgumentParse.parse_args(parser, String[])
end

@testset "keyword argument default value when missing" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "--my-arg", default="hi there")
    parsed_args = ArgumentParse.parse_args(parser, String[])
    @test parsed_args["my-arg"] == "hi there"
end

@testset "keyword argument :store_true" begin
    parser = ArgumentParse.ArgumentParser(propagate_errors=true)
    ArgumentParse.add_argument(parser, "--my-arg", action=:store_true)

    parsed_args = ArgumentParse.parse_args(parser, String[])
    @test parsed_args["my-arg"] == false

    parsed_args = ArgumentParse.parse_args(parser, String["--my-arg"])
    @test parsed_args["my-arg"] == true
end

@testset "help message" begin
    parser = ArgumentParse.ArgumentParser(program_name="test.jl", propagate_errors=true)
    ArgumentParse.add_argument(parser, "--my-keyword", default="hello")
    ArgumentParse.add_argument(parser, "my_arg")
    help_message = nothing
    try
        ArgumentParse.parse_args(parser, ["--help"])
        @test false
    catch exc
        if !isa(exc, ArgumentParse.ShowHelp)
            rethrow()
        end
        help_message = exc.message
    end
    @test help_message == "Usage: test.jl [--my-keyword my-keyword] my_arg"
end

end


text:
export label2vec,vec2label

"label=2,class=3 => label2vec(outputBase,class) = [0,1,0]"
function label2vec(label,class)
  v=zeros(class,length(label))
  for idx=1:length(label)
    v[round(Int,label[idx]),idx]=1
  end
  return v
end

"vec2label([1. 0 1;0 1 0])"
function vec2label(vector)
~,index=findmax(vector,1)
sub=ind2sub(size(vector),vec(index))[1]
end


text:
mutable struct RelaxedRef{T}
    @atomic value::T
    RelaxedRef{T}() where {T} = new{T}()
    RelaxedRef{T}(value) where {T} = new{T}(value)
end

Base.getindex(ref::RelaxedRef) = @atomic :monotonic ref.value

function Base.setindex!(ref::RelaxedRef{T}, value) where {T}
    value = convert(T, value)
    @atomic :monotonic ref.value = value
end

# Use a better algorithm?
# https://en.wikipedia.org/wiki/String-searching_algorithm
function has_subseq(seq::AbstractVector{T}, sub::AbstractVector{T}) where {T}
    for i in firstindex(seq):lastindex(seq)-length(sub)
        d = 0
        for j in eachindex(sub)
            if !(@inbounds isequal(seq[i+d], sub[j]))
                @goto failed
            end
            d += 1
        end
        return true
        @label failed
    end
    return false
end


text:
"""
    find_leaves(lnr::OptimalTreeLearner)

Finds all leaves of OptimalTreeLearner.
"""
find_leaves(lnr::OptimalTreeLearner) = [i for i=1:IAI.get_num_nodes(lnr) if IAI.is_leaf(lnr, i)]

"""
    pwl_constraint_data(lnr::IAI.OptimalTreeLearner, vks)

Returns the regression weights from an OptimalTreeLearner by leaf.
Arguments:
- lnr: OptimalTreeLearner
- vks: headers of DataFrame X, i.e. varkeys
Returns Dict[leaf_number] containing [B0, B].
"""
function pwl_constraint_data(lnr::IAI.OptimalTreeLearner, vks)
    all_leaves = find_leaves(lnr)
    pwlConstraintDict = Dict()
    for i = 1:size(all_leaves, 1)
        β0 = IAI.get_regression_constant(lnr, all_leaves[i])
        weights = IAI.get_regression_weights(lnr, all_leaves[i])[1]
        β = []
        for i = 1:size(vks, 1)
            if vks[i] in keys(weights)
                append!(β, weights[vks[i]])
            else
                append!(β, 0.0)
            end
        end
        pwlConstraintDict[all_leaves[i]] = [β0, β]
    end
    return pwlConstraintDict
end

"""
    trust_region_data(lnr:: IAI.OptimalTreeLearner, vks)

Returns the hyperplane data from an OptimalTreeLearner.
Arguments:
- lnr: OptimalTreeLearner
- vks: headers of DataFrame X, i.e. varkeys
Returns Dict[leaf_number] containing [B0, B].
"""
function trust_region_data(lnr:: IAI.OptimalTreeLearner, vks)
    all_leaves = find_leaves(lnr)
    upperDict = Dict()
    lowerDict = Dict()
    for i = 1:size(all_leaves, 1)
        # Find all parents
        parents = [all_leaves[i]]
        while IAI.get_depth(lnr, parents[end]) > 0
            append!(parents, IAI.get_parent(lnr, parents[end]))
        end
        upperDict[all_leaves[i]] = []
        lowerDict[all_leaves[i]] = []
        for j in parents[2:end]
            # For each parent, define trust region with binary variables
            threshold = IAI.get_split_threshold(lnr, j)
            if IAI.is_hyperplane_split(lnr, j)
                weights = IAI.get_split_weights(lnr, j)
                weights = weights[1]
            else
                feature = IAI.get_split_feature(lnr, j)
                weights = Dict(feature => 1)
            end
            upper = IAI.get_upper_child(lnr, j) in parents # Checking upper vs. lower split
            α = []
            for i = 1:size(vks, 1)
                if vks[i] in keys(weights)
                    append!(α, weights[vks[i]])
                else
                    append!(α, 0.0)
                end
            end
            if upper
                append!(upperDict[all_leaves[i]], [[threshold, α]])
            else
                append!(lowerDict[all_leaves[i]], [[threshold, α]])
            end
        end
    end
    return upperDict, lowerDict
end

""" Checks if a learner is trained. """
function check_if_trained(lnr::IAI.OptimalTreeLearner)
    try
        n_nodes = IAI.get_num_nodes(lnr);
    catch err
        if isa(err, UndefRefError)
            throw(OCTHaGOnException("Trees require training before being used in constraints!"))
        else
            rethrow(err)
        end
    end
end


text:
module PyCallTweaks

import Base: print

using PyCall

"""
    print(io, pyobj)

Print the `str()` form of the Python object.
"""
Base.print(io::IO, obj::PyCall.PyObject) = print(io, py"str"(obj))

end


text:
##############################################################################
#
# In this test problem we consider a quadratic penalty method.
# This example features an algorithm with the 3 steps:
# the penalization - the unconstrained min - the 1d min
#
# Note that there is no optimization of the evaluations here.
# The penalization gives an approximation of the gradients, multipliers...
#
# Note the use of a structure for the algorithmic parameters which is
# forwarded to all the 3 steps. If a parameter is not mentioned, then the default
# entry in the algorithm will be taken.
#
#############################################################################

#include("uncons.jl")

##############################################################################
#
# Quadratic penalty algorithm
# fill_in! used instead of update! (works but usually more costly in evaluations)
# subproblems are solved via Newton method
#
#############################################################################
function penalty(stp :: NLPStopping; rho0 = 1.0, rho_min = 1e-10,
                                     rho_update = 0.5, prms = nothing)

 #algorithm's parameters
 rho = rho0

 #First call to the stopping
 #Becareful here, some entries must be filled in first.
 fill_in!(stp, stp.current_state.x)
 OK = start!(stp)

 #prepare the subproblem stopping:
 sub_nlp_at_x = NLPAtX(stp.current_state.x)
 sub_pb  = ADNLPModel(x -> obj(stp.pb, x)
                      + 1/rho * norm(max.(cons(stp.pb, x) - stp.pb.meta.ucon, 0.0))^2
                      + 1/rho * norm(max.(- cons(stp.pb, x) + stp.pb.meta.lcon, 0.0))^2,  x0)
 sub_stp = NLPStopping(sub_pb, unconstrained_check,
                               sub_nlp_at_x, main_stp = stp)

 #main loop
 while !OK

  #solve the subproblem
  reinit!(sub_stp)
  sub_stp.meta.atol = min(rho, sub_stp.meta.atol)
  global_newton(sub_stp, prms)

  #Update all the entries of the State
  fill_in!(stp, sub_stp.current_state.x)

  #Either stop! is true OR the penalty parameter is too small
  if rho < rho_min stp.meta.suboptimal = true end
  OK = stop!(stp)

  @show stp.meta.nb_of_stop, OK, rho

  #update the penalty parameter if necessary
  if !OK
   rho = rho * rho_update
   sub_stp.pb  = ADNLPModel(x -> obj(stp.pb, x)
                            + 1/rho * norm(max.(cons(stp.pb, x) - stp.pb.meta.ucon, 0.0))^2
                            + 1/rho * norm(max.(- cons(stp.pb, x) + stp.pb.meta.lcon, 0.0))^2,  x0)
  end
 end

 return stp
end

##############################################################################
#
# Quadratic penalty algorithm: buffer function
#
#############################################################################
function penalty(stp :: NLPStopping, prms)

 #extract required values in the prms file
 r0 = :rho0       ∈ fieldnames(typeof(prms)) ? prms.rho0       : 1.0
 rm = :rho_min    ∈ fieldnames(typeof(prms)) ? prms.rho_min    : 1e-10
 ru = :rho_update ∈ fieldnames(typeof(prms)) ? prms.rho_update : 0.5

 return penalty(stp, rho0 = r0, rho_min = rm, ru = 0.5, prms = prms)
end

##############################################################################
#
# Algorithmic parameters structure
#
#############################################################################

mutable struct Param

    #parameters for the penalty
    rho0       :: Float64 #initial value of the penalty parameter
    rho_min    :: Float64 #smallest possible parameter
    rho_update :: Float64 #update of the penalty parameter

    #parameters of the unconstrained minimization
    armijo_prm  :: Float64 #Armijo parameter
    wolfe_prm   :: Float64 #Wolfe parameter
    onedsolve   :: Function #1D solver
    ls_func     :: Function

    #parameters of the 1d minimization
    back_update :: Float64 #backtracking update

    function Param(;rho0        :: Float64 = 1.0,
                    rho_min     :: Float64 = sqrt(eps(Float64)),
                    rho_update  :: Float64 = 0.5,
                    armijo_prm  :: Float64 = 0.01,
                    wolfe_prm   :: Float64 = 0.99,
                    onedsolve   :: Function = backtracking_ls,
                    ls_func     :: Function = (x,y)-> armijo(x,y, τ₀ = armijo_prm),
                    back_update :: Float64 = 0.5)
        return new(rho0, rho_min, rho_update,
                   armijo_prm, wolfe_prm, onedsolve, ls_func,
                   back_update)
    end
end


text:
using MonteCarlo, MonteCarloObservable
using Test
using Random

@testset "All Tests" begin
    @testset "Models" begin
        include("modeltests_IsingModel.jl")
        include("modeltests_HubbardModelAttractive.jl")
    end

    @testset "Flavors" begin
        # include("flavortests_MC.jl")
        include("flavortests_DQMC.jl")
    end

    @testset "Examples" begin
        @testset "Ising Model + MC" begin
            Random.seed!(123)
            m = IsingModel(dims=2, L=8);
            mc = MC(m, beta=0.35);
            run!(mc, sweeps=1000, thermalization=10, verbose=false);
            m = mc.obs["m"] # magnetization

            @test isapprox(0.398, round(mean(m), digits=3))
            @test isapprox(0.013, round(std_error(m), digits=3))
            @test typeof(observables(mc)) == Dict{String, String}
        end
    end

    @testset "Exact Diagonalization" begin
        include("ED/ED_tests.jl")
    end
end

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
### CONCRETE TYPE: ITERATIVE PROX EVALUATION
# prox! is computed using CG on a system with matrix lambda*A'A + I/gamma
# or lambda*AA' + I/gamma, according to which matrix is smaller.

using LinearAlgebra

struct LeastSquaresIterative{N, R <: Real, RC <: RealOrComplex{R}, M, V <: AbstractArray{RC, N}, O} <: LeastSquares
    A::M # m-by-n operator
    b::V # m (by-p)
    lambda::R
    lambdaAtb::V
    shape::Symbol
    S::O
    res::Array{RC, N} # m (by-p)
    res2::Array{RC, N} # m (by-p)
    q::Array{RC, N} # n (by-p)
end

is_prox_accurate(f::LeastSquaresIterative) = false
is_convex(f::LeastSquaresIterative) = f.lambda >= 0
is_concave(f::LeastSquaresIterative) = f.lambda <= 0

function LeastSquaresIterative(A::M, b::V, lambda::R) where {N, R <: Real, RC <: RealOrComplex{R}, M, V <: AbstractArray{RC, N}}
    if size(A, 1) != size(b, 1)
        error("A and b have incompatible dimensions")
    end
    m, n = size(A)
    x_shape = infer_shape_of_x(A, b)
    if m >= n
        shape = :Tall
        S = AcA(A)
        LeastSquaresIterative{N, R, RC, M, V, AcA}(A, b, lambda, lambda*(A'*b), shape, S, zero(b), [], zeros(RC, x_shape))
    else
        shape = :Fat
        S = AAc(A)
        LeastSquaresIterative{N, R, RC, M, V, AAc}(A, b, lambda, lambda*(A'*b), shape, S, zero(b), zero(b), zeros(RC, x_shape))
    end
end

function (f::LeastSquaresIterative{N, R, RC, M, V})(x::AbstractArray{RC, N}) where {N, R, RC, M, V}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function prox!(y::AbstractArray{D, N}, f::LeastSquaresIterative{N, R, RC, M, V}, x::AbstractArray{D, N}, gamma::R=R(1)) where {N, R, RC, M, V, D <: RealOrComplex{R}}
    f.q .= f.lambdaAtb .+ x./gamma
    # two cases: (1) tall A, (2) fat A
    if f.shape == :Tall
        y .= x
        op = ScaleShift(RC(f.lambda), f.S, RC(1)/gamma)
        IterativeSolvers.cg!(y, op, f.q)
    else # f.shape == :Fat
        # y .= gamma*(f.q - lambda*(f.A'*(f.fact\(f.A*f.q))))
        mul!(f.res, f.A, f.q)
        op = ScaleShift(RC(f.lambda), f.S, RC(1)/gamma)
        IterativeSolvers.cg!(f.res2, op, f.res)
        mul!(y, adjoint(f.A), f.res2)
        y .*= -f.lambda
        y .+= f.q
        y .*= gamma
    end
    mul!(f.res, f.A, y)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function gradient!(y::AbstractArray{D, N}, f::LeastSquaresIterative{N, R, RC, M, V}, x::AbstractArray{D, N}) where {N, R, RC, M, V, D <: Union{R, Complex{R}}}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    mul!(y, adjoint(f.A), f.res)
    y .*= f.lambda
    fy = (f.lambda/2)*dot(f.res, f.res)
end

function prox_naive(f::LeastSquaresIterative{N}, x::AbstractArray{D, N}, gamma::R=R(1)) where {N, R, D <: RealOrComplex{R}}
    y = IterativeSolvers.cg(f.lambda*f.A'*f.A + I/gamma, f.lambda*f.A'*f.b + x/gamma)
    fy = (f.lambda/2)*norm(f.A*y-f.b)^2
    return y, fy
end


text:
# these functions return vectors of states, actions and observations, ordered according to stateindex, actionindex, etc.

"""
    ordered_actions(mdp)

Return an `AbstractVector` of actions ordered according to `actionindex(mdp, a)`.

`ordered_actions(mdp)` will always return an `AbstractVector{A}` `v` containing all of the actions in `actions(mdp)` in the order such that `actionindex(mdp, v[i]) == i`. You may wish to override this for your problem for efficiency.
"""
ordered_actions(mdp::Union{MDP,POMDP}) = ordered_vector(actiontype(typeof(mdp)), a->actionindex(mdp,a), actions(mdp), n_actions(mdp), "action")

"""
    ordered_states(mdp)

Return an `AbstractVector` of states ordered according to `stateindex(mdp, a)`.

`ordered_states(mdp)` will always return a `AbstractVector{A}` `v` containing all of the states in `states(mdp)` in the order such that `stateindex(mdp, v[i]) == i`. You may wish to override this for your problem for efficiency.
"""
ordered_states(mdp::Union{MDP,POMDP}) = ordered_vector(statetype(typeof(mdp)), s->stateindex(mdp,s), states(mdp), n_states(mdp), "state")

"""
    ordered_observations(pomdp)

Return an `AbstractVector` of observations ordered according to `obsindex(pomdp, a)`.

`ordered_observations(mdp)` will always return a `AbstractVector{A}` `v` containing all of the observations in `observations(pomdp)` in the order such that `obsindex(pomdp, v[i]) == i`. You may wish to override this for your problem for efficiency.
"""
ordered_observations(pomdp::POMDP) = ordered_vector(obstype(typeof(pomdp)), o->obsindex(pomdp,o), observations(pomdp), n_observations(pomdp), "observation")

function ordered_vector(T::Type, index::Function, space, len, singular, plural=singular*"s")
    a = Array{T}(undef, len)
    gotten = falses(len)
    for x in space
        id = index(x)
        if id > len || id < 1
            error("""
                  $(singular)index(...) returned an index that was out of bounds for $singular $x.

                  index was $id.

                  n_$plural(...) was $len.
                  """)
        end
        a[id] = x
        gotten[id] = true
    end
    if !all(gotten)
        missing = findall(.!gotten)
        @warn """
             Problem creating an ordered vector of $plural in ordered_$plural(...). There is likely a mistake in $(singular)index(...) or n_$plural(...).

             n_$plural(...) was $len.

             $plural corresponding to the following indices were missing from $plural(...): $missing
             """
    end
    return a
end

@POMDP_require ordered_actions(mdp::Union{MDP,POMDP}) begin
    P = typeof(mdp)
    @req actionindex(::P, ::actiontype(P))
    @req n_actions(::P)
    @req actions(::P)
    as = actions(mdp)
end

@POMDP_require ordered_states(mdp::Union{MDP,POMDP}) begin
    P = typeof(mdp)
    @req stateindex(::P, ::statetype(P))
    @req n_states(::P)
    @req states(::P)
    as = states(mdp)
end

@POMDP_require ordered_observations(mdp::Union{MDP,POMDP}) begin
    P = typeof(mdp)
    @req obsindex(::P, ::obstype(P))
    @req n_observations(::P)
    @req observations(::P)
    as = observations(mdp)
end


text:
###############################################################################
#
# Install cdsapi package for Conda Python
#
###############################################################################
"""
    install_cdsapi!()

Install cdsapi package if it does not exist
"""
function install_cdsapi!()
    @info "Install cdsapi from conda-forge channel";
    Conda.add("cdsapi"; channel="conda-forge");

    return nothing
end


text:
using JDF
using Test
using DataFrames
using Random:randstring
using WeakRefStrings

include("test-categorical-ararys.jl")
include("test-stringarray.jl")
include("test-symbol.jl")
include("test-jdf-str.jl")
include("test-eachcol.jl")
include("test-ZonedDateTime.jl")
include("test-dataframe-syntax.jl")

@testset "JDF.jl parallel" begin
    df = DataFrame([collect(1:100) for i =1:3000])
    df[!, :int_missing] =
        rand(
            [rand(rand([UInt, Int, Float64, Float32, Bool])), missing],
            nrow(df)
        )

    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i in 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i in 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i in 1:nrow(df)]
    df[!, :stringarray_missing] = StringVector([rand([missing, randstring(8)]) for i in 1:nrow(df)])
    df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i in 1:nrow(df)]
    df[!, :char] = getindex.(df[!, :strs], 1)
    df[!, :char_missing] = allowmissing(df[!, :char])
    df[rand(1:nrow(df), 10), :char_missing] .= missing

    savejdf("a.jdf", df)

    df2 = loadjdf("a.jdf")
    @test ncol(df2) == 3009
    @test nrow(df2) == 100

    df2[!, :stringarray_missing]

    @test all(all(isequal(df[!,n], df2[!, n])) for n in names(df))

    # clean up
    rm("a.jdf", force=true, recursive=true)
end

@testset "JDF.jl serial" begin
    df = DataFrame([collect(1:100) for i =1:3000])
    df[!, :int_missing] =
        rand(
            [rand(rand([UInt, Int, Float64, Float32, Bool])), missing],
            nrow(df)
        )

    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i in 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i in 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i in 1:nrow(df)]
    df[!, :stringarray_missing] = StringVector([rand([missing, randstring(8)]) for i in 1:nrow(df)])
    df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i in 1:nrow(df)]
    df[!, :char] = getindex.(df[!, :strs], 1)
    df[!, :char_missing] = allowmissing(df[!, :char])
    df[rand(1:nrow(df), 10), :char_missing] .= missing

    ssavejdf("a.jdf", df)

    df2 = sloadjdf("a.jdf")
    @test ncol(df2) == 3009
    @test nrow(df2) == 100

    @test all(all(isequal(df[!,n], df2[!, n])) for n in names(df))

    # clean up
    rm("a.jdf", force=true, recursive=true)
end


text:
module ProfileCanvas

using Profile, JSON, REPL, Pkg.Artifacts

export @profview, @profview_allocs

struct ProfileData
    data
    typ
end

mutable struct ProfileFrame
    func::String
    file::String # human readable file name
    path::String # absolute path
    line::Int # 1-based line number
    count::Int # number of samples in this frame
    countLabel::Union{Missing,String} # defaults to `$count samples`
    flags::UInt8 # any or all of ProfileFrameFlag
    taskId::Union{Missing,UInt}
    children::Vector{ProfileFrame}
end

struct ProfileDisplay <: Base.Multimedia.AbstractDisplay end

function __init__()
    pushdisplay(ProfileDisplay())

    atreplinit(
        i -> begin
            while ProfileDisplay() in Base.Multimedia.displays
                popdisplay(ProfileDisplay())
            end
            pushdisplay(ProfileDisplay())
        end
    )
end

function Base.show(io::IO, ::MIME"text/html", canvas::ProfileData)
    id = "profiler-container-$(round(Int, rand()*100000))"

    rootpath = artifact"jlprofilecanvas"
    path = joinpath(rootpath, "jl-profile.js-0.5.2", "dist", "profile-viewer.js")

    println(
        io,
        """

        """
    )
end

function Base.display(_::ProfileDisplay, canvas::ProfileData)
    rootpath = artifact"jlprofilecanvas"
    path = joinpath(rootpath, "jl-profile.js-0.5.2", "dist", "profile-viewer.js")

    file = string(tempname(), ".html")
    open(file, "w") do io
        id = "profiler-container-$(round(Int, rand()*100000))"

        println(
            io,
            """

            """
        )
    end
    url = "file://$file"

    if Sys.iswindows()
        run(`cmd /c "start $url"`)
    elseif Sys.isapple()
        run(`open $url`)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $url`)
    end
end
using Profile

# https://github.com/timholy/FlameGraphs.jl/blob/master/src/graph.jl
const ProfileFrameFlag = (
    RuntimeDispatch=UInt8(2^0),
    GCEvent=UInt8(2^1),
    REPL=UInt8(2^2),
    Compilation=UInt8(2^3),
    TaskEvent=UInt8(2^4)
)

function view(data=Profile.fetch(); C=false, kwargs...)
    d = Dict{String,ProfileFrame}()

    if VERSION >= v"1.8.0-DEV.460"
        threads = ["all", 1:Threads.nthreads()...]
    else
        threads = ["all"]
    end

    if isempty(data)
        Profile.warning_empty()
        return
    end

    lidict = Profile.getdict(unique(data))
    data_u64 = convert(Vector{UInt64}, data)
    for thread in threads
        graph = stackframetree(data_u64, lidict; thread=thread, kwargs...)
        d[string(thread)] = make_tree(
            ProfileFrame(
                "root", "", "", 0, graph.count, missing, 0x0, missing, ProfileFrame[]
            ), graph; C=C, kwargs...)
    end

    return ProfileData(d, "Thread")
end

function stackframetree(data_u64, lidict; thread=nothing, combine=true, recur=:off)
    root = combine ? Profile.StackFrameTree{StackTraces.StackFrame}() : Profile.StackFrameTree{UInt64}()
    if VERSION >= v"1.8.0-DEV.460"
        thread = thread == "all" ? (1:Threads.nthreads()) : thread
        root, _ = Profile.tree!(root, data_u64, lidict, true, recur, thread)
    else
        root = Profile.tree!(root, data_u64, lidict, true, recur)
    end
    if !isempty(root.down)
        root.count = sum(pr -> pr.second.count, root.down)
    end

    return root
end

function status(sf::StackTraces.StackFrame)
    st = UInt8(0)
    if sf.from_c && (sf.func === :jl_invoke || sf.func === :jl_apply_generic || sf.func === :ijl_apply_generic)
        st |= ProfileFrameFlag.RuntimeDispatch
    end
    if sf.from_c && startswith(String(sf.func), "jl_gc_")
        st |= ProfileFrameFlag.GCEvent
    end
    if !sf.from_c && sf.func === :eval_user_input && endswith(String(sf.file), "REPL.jl")
        st |= ProfileFrameFlag.REPL
    end
    if !sf.from_c && occursin("./compiler/", String(sf.file))
        st |= ProfileFrameFlag.Compilation
    end
    if !sf.from_c && occursin("task.jl", String(sf.file))
        st |= ProfileFrameFlag.TaskEvent
    end
    return st
end

function status(node::Profile.StackFrameTree, C::Bool)
    st = status(node.frame)
    C && return st
    # If we're suppressing C frames, check all C-frame children
    for child in values(node.down)
        child.frame.from_c || continue
        st |= status(child, C)
    end
    return st
end

function add_child(graph::ProfileFrame, node, C::Bool)
    name = string(node.frame.file)
    func = String(node.frame.func)

    if func == ""
        func = "unknown"
    end

    frame = ProfileFrame(
        func,
        basename(name),
        name,
        node.frame.line,
        node.count,
        missing,
        status(node, C),
        missing,
        ProfileFrame[]
    )

    push!(graph.children, frame)

    return frame
end

function make_tree(graph, node::Profile.StackFrameTree; C=false)
    for child_node in sort!(collect(values(node.down)); rev=true, by=node -> node.count)
        # child not a hidden frame
        if C || !child_node.frame.from_c
            child = add_child(graph, child_node, C)
            make_tree(child, child_node; C=C)
        else
            make_tree(graph, child_node)
        end
    end

    return graph
end

"""
    @profview f(args...) [C = false]

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.

The default of `C = false` will only show Julia frames in the profile graph.
"""
macro profview(ex, args...)
    return quote
        Profile.clear()
        Profile.@profile $(esc(ex))
        view(; $(esc.(args)...))
    end
end

## Allocs

"""
    @profview_allocs f(args...) [sample_rate=0.0001] [C=false]

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.
"""
macro profview_allocs(ex, args...)
    sample_rate_expr = :(sample_rate = 0.0001)
    for arg in args
        if Meta.isexpr(arg, :(=)) && length(arg.args) > 0 && arg.args[1] === :sample_rate
            sample_rate_expr = arg
        end
    end
    if isdefined(Profile, :Allocs)
        return quote
            Profile.Allocs.clear()
            Profile.Allocs.@profile $(esc(sample_rate_expr)) $(esc(ex))
            view_allocs()
        end
    else
        return :(@error "This version of Julia does not support the allocation profiler.")
    end
end

function view_allocs(_results=Profile.Allocs.fetch(); C=false)
    results = _results::Profile.Allocs.AllocResults
    allocs = results.allocs

    allocs_root = ProfileFrame("root", "", "", 0, 0, missing, 0x0, missing, ProfileFrame[])
    counts_root = ProfileFrame("root", "", "", 0, 0, missing, 0x0, missing, ProfileFrame[])
    for alloc in allocs
        this_allocs = allocs_root
        this_counts = counts_root

        for sf in Iterators.reverse(alloc.stacktrace)
            if !C && sf.from_c
                continue
            end
            file = string(sf.file)
            this_counts′ = ProfileFrame(
                string(sf.func), basename(file), file,
                sf.line, 0, missing, 0x0, missing, ProfileFrame[]
            )
            ind = findfirst(c -> (
                    c.func == this_counts′.func &&
                    c.path == this_counts′.path &&
                    c.line == this_counts′.line
                ), this_allocs.children)

            this_counts, this_allocs = if ind === nothing
                push!(this_counts.children, this_counts′)
                this_allocs′ = deepcopy(this_counts′)
                push!(this_allocs.children, this_allocs′)

                (this_counts′, this_allocs′)
            else
                (this_counts.children[ind], this_allocs.children[ind])
            end
            this_allocs.count += alloc.size
            this_allocs.countLabel = memory_size(this_allocs.count)
            this_counts.count += 1
        end

        alloc_type = replace(string(alloc.type), "Profile.Allocs." => "")
        ind = findfirst(c -> (c.func == alloc_type), this_allocs.children)
        if ind === nothing
            push!(this_allocs.children, ProfileFrame(
                alloc_type, "", "",
                0, this_allocs.count, memory_size(this_allocs.count), ProfileFrameFlag.GCEvent, missing, ProfileFrame[]
            ))
            push!(this_counts.children, ProfileFrame(
                alloc_type, "", "",
                0, 1, missing, ProfileFrameFlag.GCEvent, missing, ProfileFrame[]
            ))
        else
            this_counts.children[ind].count += 1
            this_allocs.children[ind].count += alloc.size
            this_allocs.children[ind].countLabel = memory_size(this_allocs.count)
        end

        counts_root.count += 1
        allocs_root.count += alloc.size
        allocs_root.countLabel = memory_size(allocs_root.count)
    end

    d = Dict{String,ProfileFrame}(
        "size" => allocs_root,
        "count" => counts_root
    )

    return ProfileData(d, "Allocation")
end

const prefixes = ["bytes", "KB", "MB", "GB", "TB", "PB", "EB", "ZB", "YB"]
function memory_size(size)
    i = 1
    while size > 1000 && i + 1 < length(prefixes)
        size /= 1000
        i += 1
    end
    return string(round(Int, size), " ", prefixes[i])
end

end


text:
@testset "Test Costs" begin
    M = Sphere(2)
    N = Power(M,(3,3))
    f = PowPoint(repeat([SnPoint([1.,0.,0.])],3,3))
    x = PowPoint(repeat([SnPoint([1.,0.,0.])],3,3))
    x[2,1] = SnPoint([0., 1., 0.])

    @test costIntrICTV12(N,f,f,f,0.,0.) == 0.
    @test costIntrICTV12(N,f,x,x,0.,0.) == 1/2*distance(N, midPoint(N,x,x,f),f)^2
    @test costIntrICTV12(N, x, x, x, 2., .5) ≈ costTV2(N,x) + costTV(N,x)
    @test costIntrICTV12(N, x, x, x, 1., 0.) ≈ costTV2(N,x)
    @test costIntrICTV12(N, x, x, x, 1., 1.) ≈ costTV(N,x)
    @test costL2TV2(N,f,1.,x) == 1/2*distance(N,f,x)^2 + 1. *costTV2(N,x)
    #
    @test costL2TV(N,f,1.,f) ≈ 0.
    @test costL2TV(N,x,1.,x) ≈ 3*π/2
    @test costL2TVTV2(N, f, 0., 1., x) ≈ 1/2*distance(N,x,f)^2 + costTV2(N,x)
    @test costL2TVTV2(N, f, 1., 1., x) ≈ 1/2*distance(N,x,f)^2 + costTV(N,x) + costTV2(N,x)

    @test costTV2(M, Tuple(getValue(x[1:3,1]))) ≈ π/2
    @test costTV(N,x,1,2) ≈ sqrt(5/4)*π
    @test sum( costTV2(N,x,1,false) ) == costTV2(N,x)
    @test costTV2(N,f,2) == 0
end