text:
export VC_model
"""
    VC_model(df, y, grouping, d, link)

Form the variance component model (VCM) for intercept only regression with a random intercept covariance matrix and the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance components vector.
    One can put true (e.g. `penalized = true`) to add this penalty for numerical stability (default `penalized = false`).
"""
function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        V = [ones(di, di)]
        gcs[i] = GaussianCopulaVCObs(Y, X, V)
    end
    gcm = GaussianCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        V = [ones(di, di)]
        gcs[i] = GLMCopulaVCObs(Y, X, V, d, link)
    end
    gcm = GLMCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        V = [ones(di, di)]
        gcs[i] = NBCopulaVCObs(Y, X, V, d, link)
    end
    gcm = NBCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

"""
    VC_model(df, y, grouping, covariates, d, link)

Form the variance component model (VCM) for regression with a random intercept covariance matrix and the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `covariates`: Covariate names of interest as a vector of `Symbol`s.
    Each variable name must be present in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance components vector.
    One can put true (e.g. `penalized = true`) to add this penalty for numerical stability (default `penalized = false`).
"""
function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        V = [ones(di, di)]
        gcs[i] = GaussianCopulaVCObs(Y, X, V)
    end
    gcm = GaussianCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        V = [ones(di, di)]
        gcs[i] = GLMCopulaVCObs(Y, X, V, d, link)
    end
    gcm = GLMCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        V = [ones(di, di)]
        gcs[i] = NBCopulaVCObs(Y, X, V, d, link)
    end
    gcm = NBCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

"""
    VC_model(df, y, grouping, covariates, V, d, link)

Form the variance component model (VCM) for intercept only regression with the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `V`: Vector of Vector of Positive Semi-Definite (PSD) Covariance Matrices. `V` is of length n, where n is the number of groups/clusters.
    Each element of `V` is also a `Vector`, but of length m. Here m is the number of variance components.
    Each element of `V` is a `Vector` of d_i x d_i PSD covariance matrices under the VCM framework,
    where d_i is the cluster size of the ith cluster, which may vary for each cluster of observations i in [1, n].
    Each of these dimensions must match that specified in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance components vector.
    One can put true (e.g. `penalized = true`) to add this penalty for numerical stability (default `penalized = false`).
"""
function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    V::Vector{Vector{Matrix{Float64}}},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = GaussianCopulaVCObs(Y, X, V[i])
    end
    gcm = GaussianCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    V::Vector{Vector{Matrix{Float64}}},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = GLMCopulaVCObs(Y, X, V[i], d, link)
    end
    gcm = GLMCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    V::Vector{Vector{Matrix{Float64}}},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = NBCopulaVCObs(Y, X, V[i], d, link)
    end
    gcm = NBCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

"""
    VC_model(df, y, grouping, covariates, V, d, link)

Form the variance component model (VCM) for regression with the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `covariates`: Covariate names of interest as a vector of `Symbol`s.
    Each variable name must be present in `df`.
- `V`: Vector of Vector of Positive Semi-Definite (PSD) Covariance Matrices. `V` is of length n, where n is the number of groups/clusters.
    Each element of `V` is also a `Vector`, but of length m. Here m is the number of variance components.
    Each element of `V` is a `Vector` of d_i x d_i PSD covariance matrices under the VCM framework,
    where d_i is the cluster size of the ith cluster, which may vary for each cluster of observations i in [1, n].
    Each of these dimensions must match that specified in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance components vector.
    One can put true (e.g. `penalized = true`) to add this penalty for numerical stability (default `penalized = false`).
"""
function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    V::Vector{Vector{Matrix{Float64}}},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaVCObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = GaussianCopulaVCObs(Y, X, V[i])
    end
    gcm = GaussianCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    V::Vector{Vector{Matrix{Float64}}},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = GLMCopulaVCObs(Y, X, V[i], d, link)
    end
    gcm = GLMCopulaVCModel(gcs; penalized = penalized);
    return gcm
end

function VC_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    V::Vector{Vector{Matrix{Float64}}},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaVCObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = NBCopulaVCObs(Y, X, V[i], d, link)
    end
    gcm = NBCopulaVCModel(gcs; penalized = penalized);
    return gcm
end


text:
array_value(e::Identity) = Identity(e.group, array_value(e.p))

array_point(p) = ValidationMPoint(p)
array_point(p::ValidationMPoint) = p
array_point(e::Identity) = Identity(e.group, array_point(e.p))

function Base.inv(M::ValidationManifold, p; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    q = array_point(inv(M.manifold, array_value(p)))
    is_manifold_point(M, q, true; kwargs...)
    return q
end

function inv!(M::ValidationManifold, q, p; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    inv!(M.manifold, array_value(q), array_value(p))
    is_manifold_point(M, q, true; kwargs...)
    return q
end

function Base.identity(M::ValidationManifold, p; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    q = array_point(identity(M.manifold, array_value(p)))
    is_manifold_point(M, q, true; kwargs...)
    return q
end

function identity!(M::ValidationManifold, q, p; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    identity!(M.manifold, array_value(q), array_value(p))
    is_manifold_point(M, q, true; kwargs...)
    return q
end

function compose(M::ValidationManifold, p, q; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    x = array_point(compose(M.manifold, array_value(p), array_value(q)))
    is_manifold_point(M, x, true; kwargs...)
    return x
end

function compose!(M::ValidationManifold, x, p, q; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    compose!(M.manifold, array_value(x), array_value(p), array_value(q))
    is_manifold_point(M, x, true; kwargs...)
    return x
end

function translate(M::ValidationManifold, p, q, conv::ActionDirection; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    x = array_point(translate(M.manifold, array_value(p), array_value(q), conv))
    is_manifold_point(M, x, true; kwargs...)
    return x
end

function translate!(M::ValidationManifold, x, p, q, conv::ActionDirection; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    translate!(M.manifold, array_value(x), array_value(p), array_value(q), conv)
    is_manifold_point(M, x, true; kwargs...)
    return x
end

function inverse_translate(M::ValidationManifold, p, q, conv::ActionDirection; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    x = array_point(inverse_translate(M.manifold, array_value(p), array_value(q), conv))
    is_manifold_point(M, x, true; kwargs...)
    return x
end

function inverse_translate!(
    M::ValidationManifold,
    x,
    p,
    q,
    conv::ActionDirection;
    kwargs...,
)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    inverse_translate!(M.manifold, array_value(x), array_value(p), array_value(q), conv)
    is_manifold_point(M, x, true; kwargs...)
    return x
end

function translate_diff(M::ValidationManifold, p, q, X, conv::ActionDirection; kwargs...)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    is_tangent_vector(M, q, X, true; kwargs...)
    Y = ValidationTVector(
        translate_diff(M.manifold, array_value(p), array_value(q), array_value(X), conv),
    )
    pq = translate(M, p, q, conv)
    is_tangent_vector(M, pq, Y, true; kwargs...)
    return Y
end

function translate_diff!(
    M::ValidationManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection;
    kwargs...,
)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    is_tangent_vector(M, q, X, true; kwargs...)
    translate_diff!(
        M.manifold,
        array_value(Y),
        array_value(p),
        array_value(q),
        array_value(X),
        conv,
    )
    pq = translate(M, p, q, conv)
    is_tangent_vector(M, pq, Y, true; kwargs...)
    return Y
end

function inverse_translate_diff(
    M::ValidationManifold,
    p,
    q,
    X,
    conv::ActionDirection;
    kwargs...,
)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    is_tangent_vector(M, q, X, true; kwargs...)
    Y = ValidationTVector(
        inverse_translate_diff(
            M.manifold,
            array_value(p),
            array_value(q),
            array_value(X),
            conv,
        ),
    )
    pinvq = inverse_translate(M, p, q, conv)
    is_tangent_vector(M, pinvq, Y, true; kwargs...)
    return Y
end

function inverse_translate_diff!(
    M::ValidationManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection;
    kwargs...,
)
    is_manifold_point(M, p, true; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    is_tangent_vector(M, q, X, true; kwargs...)
    inverse_translate_diff!(
        M.manifold,
        array_value(Y),
        array_value(p),
        array_value(q),
        array_value(X),
        conv,
    )
    pinvq = inverse_translate(M, p, q, conv)
    is_tangent_vector(M, pinvq, Y, true; kwargs...)
    return Y
end

function group_exp(M::ValidationManifold, X; kwargs...)
    is_tangent_vector(
        M,
        make_identity(M.manifold, array_value(X)),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    q = array_point(group_exp(M.manifold, array_value(X)))
    is_manifold_point(M, q, true; kwargs...)
    return q
end

function group_exp!(M::ValidationManifold, q, X; kwargs...)
    is_tangent_vector(
        M,
        make_identity(M.manifold, array_value(X)),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    group_exp!(M.manifold, array_value(q), array_value(X))
    is_manifold_point(M, q, true; kwargs...)
    return q
end

function group_log(M::ValidationManifold, q; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    X = ValidationTVector(group_log(M.manifold, array_value(q)))
    is_tangent_vector(
        M,
        make_identity(M.manifold, array_value(X)),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    return X
end

function group_log!(M::ValidationManifold, X, q; kwargs...)
    is_manifold_point(M, q, true; kwargs...)
    group_log!(M.manifold, array_value(X), array_value(q))
    is_tangent_vector(
        M,
        make_identity(M.manifold, array_value(X)),
        array_value(X),
        true;
        check_base_point=false,
        kwargs...,
    )
    return X
end

ext: jl


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*------------------------------------------------------------------
type DoublePipe_Geometry
	DoublePipe_Geometry()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components",
				:Hidden=>true
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Pi Number",
				:Default=>3.14159265,
				:Symbol=>"\\pi",
				:Hidden=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Outside Diameter of Inner Pipe",
				:Lower=>1e-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inside Diameter of Inner Pipe",
				:Lower=>1e-10
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inside Diameter of Outer pipe",
				:Lower=>1e-10
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Effective Tube Length of one segment of Pipe",
				:Lower=>0.1,
				:Symbol=>"L_{pipe}"
			)),
			conductivity(Dict{Symbol,Any}(
				:Brief=>"Tube Wall Material Thermal Conductivity",
				:Default=>1.0,
				:Symbol=>"K_{wall}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Inside Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Default=>1e-6,
				:Lower=>0
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Outside Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Default=>1e-6,
				:Lower=>0
			)),
			[:PP,:NComp,:M,:Pi,:DoInner,:DiInner,:DiOuter,:Lpipe,:Kwall,:Rfi,:Rfo,],
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	Pi::constant
	DoInner::length
	DiInner::length
	DiOuter::length
	Lpipe::length
	Kwall::conductivity
	Rfi::positive
	Rfo::positive
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export DoublePipe_Geometry
function set(in::DoublePipe_Geometry)
	#"Component Molecular Weight"
	M = PP.MolecularWeight()
	 #"Pi Number"
	Pi = 3.14159265
end
function atributes(in::DoublePipe_Geometry,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="double pipe geometrical parameters."
	drive!(fields,_)
	return fields
end
DoublePipe_Geometry(_::Dict{Symbol,Any})=begin
	newModel=DoublePipe_Geometry()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
import Stochy: cps, desugar, striplineinfo

x = 0
cpsid(k,x) = k(x)

id(x) = x

# Atomic Expressions
@test (@cps 0) == 0
@test (@cps x) == x
@test (@cps true) == true
@test (@cps "test") == "test"
@test (@cps :test) == :test
@test eval(cps(:(:test), :identity)) == :test
@test (@cps (()->false)()) == false
@test (@cps (x->x)(false)) == false
@test (@cps ((x,y)->x+y)(1, 2)) == 3

# Function calls
@test (@cps cpsid(0)) == 0
@test (@cps ((a...)->a)(1,2)) == (1,2)

# Primitives
@test (@cps 1+1) == 2
@test (@cps x+1) == 1

# Conditionals
@test (@cps if true 1 else 0 end) == 1
@test (@cps if false 1 else 0 end) == 0

# Blocks
@test (@cps begin 1 end) == 1
@test (@cps begin 0; 1 end) == 1

# Function definitions
let
    @cps function f(); false; end
    @test (@cps f()) == :false
end

let
    @cps function f(x) x end
    @test (@cps f(false)) == false
end

let
    @cps function f(x,y) x+y end
    @test (@cps f(1, 2)) == 3
end

# Comparisons
@test (@cps 1==1) == true
@test (@cps x==x) == true
@test (@cps 2>1) == true
@test (@cps 1<2) == true
@test (@cps 1>=1) == true
@test (@cps 1<=1) == true

# Boolean operators
@test (@cps false || false) == false
@test (@cps true || false) == true
@test (@cps false || true) == true
@test (@cps true || true) == true
@test (@cps false && false) == false
@test (@cps true && false) == false
@test (@cps false && true) == false
@test (@cps true && true) == true

# Array literals
@test (@cps [x,1,2]) == [x,1,2]
@test (@cps [x,[1,[2]]]) == [x,1,2]

# Array indexing
@test (@cps [true,false][1]) == true

# Locals
@test (@cps begin; local a=1; a+1 end) == 2
@test (@cps begin; local a=1; local b=2; a+b end) == 3
@test (@cps begin; local a=cpsid(1); a; end) == 1
@test (@cps begin; local a=begin; 1; end; a; end) == 1

# Dot syntax
module M cpsid(k,x) = k(x) end
getmod(k) = k(M)
@test (@cps M.cpsid) == M.cpsid
@test (@cps M.cpsid(0)) == 0
@test (@cps getmod().cpsid(1)) == 1

let
    @cps function fact(n)
        n == 0 ? 1 : n * fact(n - 1)
    end
    @test (@cps fact(5)) == 120
end

# Ellipsis
@test (@cps tuple([1,2]...)) == (1,2)

# Desugar
expr_eq(e,f) = striplineinfo(e) == striplineinfo(f)

for expr in [:(1),
             :(x),
             :(:test),
             :(function f() false end),
             :(()->false),
             :(begin; end),
             :(f(x)),
             :(begin; f(x, y + z); end)]
    @test desugar(expr) == expr
end

@test expr_eq(desugar(:(begin; local a=1, b=2 end)),
                      :(begin; local a=1; local b=2 end))

@test expr_eq(desugar(:(begin; local a=1,b=2; begin; local x=1, y=2; end; end)),
                      :(begin; local a=1; local b=2; begin; local x=1; local y=2; end; end))

@test expr_eq(desugar(:(function f(); local x=1, y=2 end)),
                      :(function f(); local x=1; local y=2 end))

@test expr_eq(desugar(:(local x=1, y=2)),
                      :(begin; local x=1; local y=2 end))

@test expr_eq(desugar(:(~nothing)), :(sample(nothing)))


text:
using Schwarz: Constants, PowerLawPotential, OrbitLibrary
c = Constants

r0 = 10 * c.kpc
rho0 = 4e7 * c.Msun / c.kpc^3
alpha = -2.5
phi = PowerLawPotential(r0, rho0, alpha)

Nmax = 16
reps = 5
Nv = 4
rmax = 40 * c.kpc
multithreaded = true
n = nworkers()

sizes = []
means = []
stds = []

println("Pre-compiling OrbitLibrary")
lib = OrbitLibrary(1, 1, rmax, phi; multithreaded = true)
Nr = 1
lib, elapsed = @timed OrbitLibrary(Nr, Nv, rmax, phi; multithreaded = multithreaded)

for Nr in 1:Nmax
    println("Timing OrbitLibrary")
    println("    Nr = $Nr, Nv = $Nv")
    times = Array{Float64, 1}(reps)
    for i in 1:reps
        lib, elapsed = @timed OrbitLibrary(Nr, Nv, rmax, phi; multithreaded = multithreaded)
        times[i] = elapsed
    end
    push!(sizes, Nr)
    push!(means, mean(times))
    push!(stds, std(times))
end

writedlm("orbitLibrary_$n.txt", [sizes means stds])


text:
# Check to make sure types are imported properly
@test RadioMenu <: TerminalMenus.AbstractMenu

# Constructor
@test RadioMenu(["one", "two", "three"]).pagesize == 3
@test RadioMenu(string.(1:30), pagesize=-1).pagesize == 30
@test RadioMenu(string.(1:4), pagesize=10).pagesize == 4
@test RadioMenu(string.(1:100)).pagesize == 10

radio_menu = RadioMenu(string.(1:20))
@test TerminalMenus.options(radio_menu) == string.(1:20)
radio_menu.selected = 2
TerminalMenus.cancel(radio_menu)
@test radio_menu.selected == -1
@test TerminalMenus.header(radio_menu) == ""

# Output
TerminalMenus.config() # Use default chars
CONFIG = TerminalMenus.CONFIG

# Test writeLine function
term_width = 100
radio_menu = RadioMenu(string.(1:10))
buf = IOBuffer()
TerminalMenus.writeLine(buf, radio_menu, 1, true, term_width)
@test String(take!(buf)) == string(CONFIG[:cursor], " 1")
TerminalMenus.config(cursor='+')
TerminalMenus.writeLine(buf, radio_menu, 1, true, term_width)
@test String(take!(buf)) == "+ 1"
TerminalMenus.config(charset=:unicode)
TerminalMenus.writeLine(buf, radio_menu, 1, true, term_width)
@test String(take!(buf)) == string(CONFIG[:cursor], " 1")

# Test using STDIN
radio_menu = RadioMenu(string.(1:10))
@test simulateInput(3, radio_menu, :down, :down, :enter)
radio_menu = RadioMenu(["single option"])
@test simulateInput(1, radio_menu, :up, :up, :down, :up, :enter)
radio_menu = RadioMenu(string.(1:3), pagesize=1)
@test simulateInput(3, radio_menu, :down, :down, :down, :down, :enter)


text:
# Inspired from JuMP/test/solvers

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

# Load available solvers
grb = try_import(:Gurobi)
cpx = try_import(:CPLEX)
xpr = try_import(:Xpress)
mos = false && try_import(:MosekTools)
cbc = try_import(:Cbc)
if cbc; import Clp; end

import JuMP
import GLPK
glp = true
ipt = try_import(:Ipopt)
eco = try_import(:ECOS)

# Create LP solver list
lp_solvers = Module[]
grb && push!(lp_solvers, Gurobi)
cpx && push!(lp_solvers, CPLEX)
xpr && push!(lp_solvers, Xpress)
mos && push!(lp_solvers, Mosek)
cbc && push!(lp_solvers, Clp)
glp && push!(lp_solvers, GLPK)
ipt && push!(lp_solvers, Ipopt)
eco && push!(lp_solvers, ECOS)


text:
module Pipelines

using Dates
using UUIDs
using Logging

include("utils.jl")
export do_nothing, isok,
str, to_str, to_cmd,
replaceext, removeext,
redirect_to_files

include("CmdDependency.jl")
export CmdDependency,
check_dependency,
check_dependency_file,
check_dependency_dir,
exec

include("Program.jl")
export Program,
infer_outputs

include("CmdProgram.jl")
export CmdProgram,
prepare_cmd

include("JuliaProgram.jl")
export JuliaProgram

include("pretty_print.jl")

# binding documentations
@doc (@doc _run) run

function __init__()
    Pipelines.set_default_stdout()
    Pipelines.set_default_stderr()
end
end # module


text:
using CVCore
using CVImgProc
using CVImgCodecs
using Base.Test

lenapng = Pkg.dir("OpenCV", "test", "lena.png")
lenamat = imread(lenapng)
@test size(lenamat) == (512,512,3)
@assert all(!isnan(lenamat))

@testset "UMat{T,N}" begin
    umat = UMat{UInt8}(3, 3)
    mat = Mat(umat)
    @test isa(mat, Mat)
    @test size(umat) == (3, 3)

    # make sure resize doesn't change the type of mat
    umat_resized = resize(umat, size(umat))
    @test isa(umat_resized, UMat)
    @test size(umat) == size(umat_resized)
end

@testset "cv::resize" begin
    mat = copy(lenamat)
    rows,cols,cn = size(mat)

    # shape as a tuple, acts like column-major
    mat2 = resize(mat, (200,300))
    @test size(mat2) == (200,300,3)

    # shape as a cv::Size, acts like row-major
    mat3 = resize(mat, CVCore.cvSize(200,300))
    @test size(mat3) == (300,200,3)
end


text:
module SymmetrySnap
    using LinearAlgebra
    using Random
    using DelimitedFiles
    using MinkowskiReduction

    """
        getLatticeDiff(oldA, a, oldB, b, oldC, c, output)

    Save the change in lattice to output vector

    # Arguments
    - `oldA, oldB, oldC::Vector{Float64}`: Old lattice vectors
    - `a, b, c::Vector{Float64}`: New lattice vectors
    - `output::Vector{String}`: Output list
    """
    function getLatticeDiff(oldA, a, oldB, b, oldC, c, numSymmetries, output)
        push!(output, "Old Lattice -> New Lattice : L2 Norm, SymOps = $(numSymmetries)")
        push!(output, join([oldA, a, norm(oldA - a)]," -> "," : "))
        push!(output, join([oldB, b, norm(oldB - b)]," -> "," : "))
        push!(output, join([oldC, c, norm(oldC - c)]," -> "," : "))
        push!(output, "")
        return output
    end

    """
        getAtomDiff(oldAtomPos, atomPos, symmetries, output)

    Save the change in atom positions and number of symmetries found to output vector

    # Arguments
    - `oldAtomPos::Vector{Vector{Float64}}`: Old atom positions
    - `atomPos::Vector{Vector{Float64}}`: New atom positions
    - `symmetries::Int8`: Number of symmetries found
    - `output::Vector{String}`: Output list
    """
    function getAtomDiff(oldAtomPos, atomPos, symmetries, output)
        push!(output, "Old Atom Position -> New Atom Position : L2 Norm, SymOps")
        for i in 1:length(atomPos)
            push!(output, join([join([oldAtomPos[i], atomPos[i], norm(oldAtomPos[i] - atomPos[i])]," -> "," : "),symmetries[i]],", "))
        end
        return output
    end

    """
        mapAtomsToCell(a, b, c, latVecInv, atomPositions)

    Map given atom positions into positions inside the unit cell.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `latVecInv::Array{Float64, 2}`: Inverse of lattice vectors
    - `atomPositions::Vector{Vector{Float64}}`

    # Examples
    ```jldoctest
    julia> mapAtomsToCell([1,0,0], [0,1,0], [0,0,1], [1, 0, 0; 0, 1, 0; 0, 0, 1], [1.5,-0.5,3.5])
    [0.5, 0.5, 0.5]
    ```
    """
    function mapAtomsToCell(a, b, c, latVecInv, atomPositions)
        newPositions = []
        for atom in atomPositions
            addTo = true
            shift = latVecInv * atom
            for i in 1:3
                shift[i] = floor(shift[i])
            end
            newPos = atom - a*shift[1] - b*shift[2] - c*shift[3]
            for newAtom in newPositions
                if isapprox(newPos, newAtom)
                    addTo = false
                    break
                end
            end
            if addTo == true
                push!(newPositions, newPos)
            end
        end
        return newPositions
    end

    """
        toDirect(a, b, c, atomPositions)

    Convert atom positions in Cartesian coordinates to direct coordinates.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `atomPositions::Vector{Vector{Float64}}`: atom positions in Cartesian coordinates

    # Examples

    """
    function toDirect(a, b, c, atomPositions)
        latVecs = Array{Float64, 2}(undef, 3, 3)
        latVecs[1:end,1] = a
        latVecs[1:end,2] = b
        latVecs[1:end,3] = c

        newPositions = Vector{Vector{Float64}}(undef, length(atomPositions))
        for i in 1:length(atomPositions)
            newPositions[i] = \(latVecs, atomPositions[i])
        end
        return newPositions
    end

    """
        toCartesian(a, b, c, atomPositions)

        Convert atom positions in direct coordinates to Cartesian coordinates.

        # Arguments
        - `a,b,c::Vector{Float64}`: lattice vectors
        - `atomPositions::Vector{Vector{Float64}}`: atom positions in direct coordinates

        # Examples

    """
    function toCartesian(a, b, c, atomPositions)
        latVecs = Array{Float64, 2}(undef, 3, 3)
        latVecs[1:end,1] = a
        latVecs[1:end,2] = b
        latVecs[1:end,3] = c

        newPositions = Vector{Vector{Float64}}(undef, length(atomPositions))
        for i in 1:length(atomPositions)
            newPositions[i] = latVecs * atomPositions[i]
        end
        return newPositions
    end

    """
        test_atom(a, b, c, lat_vec_inv, symmetry_op, atom_1_pos, atom_1_index, atom_2_pos, atom_2_index, natoms[, epsilon=1e-2])

    Test if atom 1 maps to the position of atom 2 when the symmetry operator is applied. If valid, returns matrix A that represents the mapping from atom 1 to atom 2.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `latVecInv::Array{Float64, 2}`: inverse lattice vectors matrix
    - `symmetryOp::Array{Float64, 2}`: symmetry operator
    - `atom_1_index, atom_2_index::{Int}`: indexes of atoms being compared
    - `atom_1_pos, atom_2_pos::Vector{Float64}`: positions of atoms being compared
    - `nAtoms::Int8`: Total number of atoms in the system
    - `epsilon::Float64`: Largest normed difference between positions for a valid symmetry operator

    # Examples

    """
    function testAtom(a, b, c, latVecInv, symmetryOp, atom1Pos, atom1Index, atom2Pos, atom2Index, nAtoms, shift, epsilon=1e-2)
        # In the case of a different crystal structure, A would be the lattice vectors as rows
        m = symmetryOp * atom1Pos #what is m??

        if shift != nothing
            m = m .+ shift
        end

        compareCell = latVecInv * m
        for i in 1:3
            compareCell[i] = floor(compareCell[i])
        end

        transformedAtom = symmetryOp * atom1Pos .- a*compareCell[1] .- b*compareCell[2] .- c*compareCell[3]
        if shift == nothing
            shift = atom2Pos .- transformedAtom
        end

        transformedAtom = transformedAtom .+ shift

        if norm(atom2Pos .- symmetryOp * atom1Pos .- shift) < epsilon || norm(atom2Pos .- transformedAtom) < epsilon
            A = zeros(Float64, 3, 3*nAtoms+1)
            if norm(atom2Pos .- symmetryOp * atom1Pos .- shift) < epsilon
                compareCell = [0,0,0]
            end
            for i in 1:3
                for j in 1:3
                    if symmetryOp[i,j] != 0
                        A[i,3*(atom1Index-1)+j] += symmetryOp[i,j]
                    end
                    if -Matrix(I, 3, 3)[i,j] != 0
                        A[i,3*(atom2Index-1)+j] += -Matrix(I, 3, 3)[i,j]
                    end
                end
            end
            A[:,3*nAtoms+1] = a*compareCell[1] .+ b*compareCell[2] .+ c*compareCell[3] .- shift
            A[abs.(A) .< 1e-6] .= 0
            if isapprox(A, zero(A), atol = 1e-10)
                return nothing, shift, true
            end
            return A, shift, true
        else
            return nothing, shift, false
        end
    end

    """
        runSymmetry(a,b,c,latVecInv, symmetryOp, atomPositions[, epsilon=1e-2])

    Test the given symmetry operator on each atom in atomPositions to verify if it is a true symmetry, and return the list of symmetry equivalent linear algebra problems.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `latVecInv::Array{Float64, 2}`: inverse lattice vectors matrix
    - `symmetryOp::Array{Float64, 2}`: symmetry operator
    - `atomPositions::Vector{Vector{Float64}}`: atom positions
    - `epsilon::Float64`: Largest normed difference between positions for a valid symmetry operator

    # Examples

    """
    function runSymmetry(a, b, c, latVecInv, symmetryOp, atomPositions, epsilon=1e-2)
        # Which atom to shift by
        for i in 1:length(atomPositions)
            As = Vector{Array{Float64, 2}}(undef, 0)
            global shift = nothing
            indexFound = []
            # Iterate over all atoms for atom 1
            for j in 1:length(atomPositions)
                global symmFound = false
                # Iterate over all atoms for atom 2
                for k in 1:length(atomPositions)
                    if i+k-1 in indexFound
                        continue
                    end
                    A, shift, symmFound = testAtom(a, b, c, latVecInv, symmetryOp, atomPositions[j], j, atomPositions[(i+k-2)%length(atomPositions)+1], (i+k-2)%length(atomPositions)+1, length(atomPositions), shift, epsilon)
                    if A != nothing
                        push!(As, A)
                        push!(indexFound, i+k-1)
                    end
                    if symmFound
                        push!(indexFound, i+k-1)
                        break
                    end
                end
                if !symmFound
                    break
                end
            end
            if symmFound
                return As
            end
        end
        return nothing
    end

    """
        symmetrize(a, b, c, lat_vec_inv, atom_positions[, epsilon=1e-2])

    Using the valid symmetry operators found, solve the \\(Ax=b\\) problem to get the smallest correction in atom positions needed to align with the improved lattice vectors.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `latVecInv::Array{Float64, 2}`: inverse lattice vectors matrix
    - `atomPositions::Vector{Vector{Float64}}`: atom positions
    - `epsilon::Float64`: Largest normed difference between positions for a valid symmetry operator

    # Examples

    """
    function symmetrize(a, b, c, latVecInv, atomPositions, symmetries, epsilon=1e-2)
        As = Vector{Array{Float64, 2}}(undef, 0)
        symm_num = 0
        for i in 1:length(symmetries)
            newA = runSymmetry(a, b, c, latVecInv, symmetries[i], atomPositions, epsilon)
            if newA != nothing
                symm_num += 1
                check = false
                for new in newA
                    for A in As
                        if isapprox(A,new, atol = 1e-10)
                            check = true
                            break
                        end
                    end
                    if !check
                        push!(As,new)
                    end
                end
            end
        end
        if length(As) == 0
            return atomPositions, zeros(length(atomPositions))
        end

        A = zeros(Float64, 3*length(As), 3*length(atomPositions))
        B = zeros(Float64, 3*length(As))
        for i in 0:length(As)-1
            A[3*i+1:3*i+3,:] = As[i+1][1:end,1:end .!= length(As[i+1][1,1:end])]
            B[3*i+1:3*i+3] = As[i+1][1:end,end]
        end

        aSVD = svd(A)
        newA = zeros(3*length(As), 3*length(atomPositions))
        for i in 1:length(aSVD.S)
            uArray = aSVD.U[1:end, i]
            vArray = aSVD.Vt[i, 1:end]
            part = aSVD.S[i]*(aSVD.U[1:end, i]*transpose(aSVD.Vt[i, 1:end]))
            part[abs.(part) .< 1e-6] .= 0
            newA += part
        end
        A = newA

        xUser = vcat(atomPositions...)
        xDiff = xUser - \(A,B)
        null = transpose(LinearAlgebra.nullspace(A))

        for i in 1:length(null[1:end, 1])
            xDiff -= null[i,1:end]*dot(xDiff,null[i,1:end])/norm(null[i,1:end])
        end
        xCorrected = xUser - xDiff
        newXs = Array{Array{Float64, 1}, 1}(undef, length(atomPositions))
        for i in 1:length(atomPositions)
            newX = Array{Float64, 1}(undef, 3)
            for j in 1:3
                newX[j] = xCorrected[3*(i-1) + j]
            end
            newXs[i] = newX
        end

        return a, b, c, newXs, repeat([symm_num], length(atomPositions))
    end

    """
        getSymmetryOps(a,b,c,latVecInv)

    Finds the applicable symmetry operators for the system by solving the following equation,

    \\(R = A B^{-1}\\),

    where \\(A\\) is the original lattice vectors, \\(B\\) is the new lattice vectors, and \\(R\\) is the symmetry operator being tested.

    If \\(R R^T = I\\), \\(R\\) is a valid symmetry operator.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `latVecInv::Array{Float64, 2}`: inverse lattice vectors matrix

    # Examples

    """
    function getSymmetryOps(a, b, c, latVecInv, atomPositions, epsilon=1e-2)
        translations = [[0,0,1],[0,1,0],[1,0,0],[0,0,-1],[0,-1,0],[-1,0,0],[1,1,0],[1,-1,0],[-1,1,0],
                        [-1,-1,0],[1,0,1],[1,0,-1],[-1,0,1],[-1,0,-1],[0,1,1],[0,1,-1],[0,-1,1],[0,-1,-1],
                        [1,1,1],[1,1,-1],[1,-1,1],[-1,1,1],[1,-1,-1],[-1,1,-1],[-1,-1,1],[-1,-1,-1]]
        latticePoints = []
        for v in translations
            push!(latticePoints,v[1]*a+v[2]*b+v[3]*c)
        end

        oldLattice = Array{Float64, 2}(undef, 3, 3)
        oldLattice[1:end, 1] = a
        oldLattice[1:end, 2] = b
        oldLattice[1:end, 3] = c

        newLattices = Vector{Array{Float64, 2}}(undef, 15600)
        transformations = Vector{Array{Float64, 2}}(undef, 15600)
        index = 1
        for i in 1:length(latticePoints)
            for j in 1:length(latticePoints)
                for k in 1:length(latticePoints)
                    if i == k || j == k || i == j
                        continue
                    end
                    newLattice = Array{Float64, 2}(undef, 3, 3)
                    newLattice[1:end, 1] = latticePoints[i]
                    newLattice[1:end, 2] = latticePoints[j]
                    newLattice[1:end, 3] = latticePoints[k]
                    newLattices[index] = newLattice

                    transformLattice = Array{Float64, 2}(undef, 3, 3)
                    transformLattice[1, 1:end] = translations[i]
                    transformLattice[2, 1:end] = translations[j]
                    transformLattice[3, 1:end] = translations[k]
                    transformations[index] = transformLattice
                    index += 1
                end
            end
        end

        allOperators = []
        allErrors = []
        for i in 1:length(newLattices)
            operator = newLattices[i]*latVecInv
            operatorT = transpose(operator)
            push!(allOperators, operator)
            push!(allErrors, norm(operator*operatorT - Matrix(I,3,3)))
        end

        p = sortperm(allErrors)
        sortOperators = allOperators[p]
        sortErrors = allErrors[p]
        sortTransforms = transformations[p]

        # Symmetrix Matrix, so only 6 elements: xx, yy, zz, xy, xz, yz
        elementMap = Dict(
            (1, 1) => 1,
            (2, 2) => 2,
            (3, 3) => 3,
            (1, 2) => 4,
            (2, 1) => 4,
            (1, 3) => 5,
            (3, 1) => 5,
            (2, 3) => 6,
            (3, 2) => 6
            )
        numOps = [2, 3, 4, 6, 8, 12, 16, 24, 48]
        lattices = []
        errors = []
        opCounts = []
        for numOp in numOps
            A = zeros(6*numOp, 6)
            for opNum in 0:numOp-1
                for i in 1:3
                    for j in i:3
                        A[6*opNum + elementMap[(i, j)], elementMap[(i, j)]] -= 1
                        for k in 1:3
                            for l in 1:3
                                A[6*opNum + elementMap[(i, j)], elementMap[(k, l)]] += sortTransforms[opNum+1][i, k]*sortTransforms[opNum+1][j, l]
                            end
                        end
                    end
                end
            end

            xUser = transpose(oldLattice)*oldLattice

            xUser = [xUser[1, 1], xUser[2, 2], xUser[3, 3], xUser[1, 2], xUser[1, 3], xUser[2, 3]]
            diff = xUser

            null = transpose(nullspace(A))

            for i in 1:length(null[1:end, 1])
                diff -= null[i,1:end]*dot(diff,null[i,1:end])/norm(null[i,1:end])
            end

            xCorrected = xUser - diff

            xCorrFull = [xCorrected[1] xCorrected[4] xCorrected[5];xCorrected[4] xCorrected[2] xCorrected[6];xCorrected[5] xCorrected[6] xCorrected[3]]

            newSVD = svd(xCorrFull)

            newSMat = Diagonal(sqrt.(newSVD.S))
            newVtMat = newSVD.Vt

            findUSVD = nothing
            try
                findUSVD = svd(oldLattice*inv(newSMat*newVtMat))
            catch
                continue
            end
            newUMat = findUSVD.U * findUSVD.Vt

            push!(lattices, newUMat*newSMat*newVtMat)
            push!(errors, norm(oldLattice - lattices[end]))
            push!(opCounts, numOp)
        end

        lattice = Array{Float64, 2}(undef, 3, 3)
        numValidOps = 0
        for i in 1:length(errors)+1
            if i == length(errors)+1 || errors[i] > epsilon
                if i == 1
                    i == 2
                end
                lattice = lattices[i-1]
                numValidOps = opCounts[i-1]
                break
            end
        end

        transformation = lattice*latVecInv
        for i in 1:length(atomPositions)
            atomPositions[i] = transformation*atomPositions[i]
        end

        latVecInv = inv(lattice)
        a = lattice[1:end, 1]
        b = lattice[1:end, 2]
        c = lattice[1:end, 3]

        validOperators = []
        for i in 1:numValidOps
            newLattice = transpose(sortTransforms[i]*transpose(lattice))
            push!(validOperators, newLattice*latVecInv)
        end

        return a, b, c, latVecInv, atomPositions, validOperators
    end

    """
        symmSnap(a, b, c, atomPositions, atomTypes[, isDirect = false, fout = "symm_snap.out", epsilon=1e-2])

    Get corrected/reduced lattice vectors, and corrected atom positions to correspond with a Bravais lattice.

    # Arguments
    - `a,b,c::Vector{Float64}`: lattice vectors
    - `atomPositions::Vector{Vector{Float64}}`: atom positions
    - `atomTypes::_______`:
    - `isDirect::Bool`: indicate if atomPositions are already in direct coordinates
    - `fout::String`: file name to write out additional information to
    - `epsilon::Float64`: Largest normed difference between positions for a valid symmetry operator

    # Examples

    """
    function symmSnap(a, b, c, atomPositions, atomTypes, isDirect = false, fout = "symm_snap.out", epsilon=1e-2)
        if isDirect
            atomPositions = toCartesian(a, b, c, atomPositions)
        end

        lattice = Array{Float64, 2}(undef, 3, 3)
        lattice[1:end, 1] = a
        lattice[1:end, 2] = b
        lattice[1:end, 3] = c
        latVecInv = inv(lattice)
        atomPositions = mapAtomsToCell(a, b, c, latVecInv, atomPositions) #new atom positions, but mapped to the unit cell

        lattice = Array{Float64, 2}(undef, 3, 3)
        lattice[1:end, 1] = a
        lattice[1:end, 2] = b
        lattice[1:end, 3] = c
        latVecInv = inv(lattice)

        oldA = deepcopy(a)
        oldB = deepcopy(b)
        oldC = deepcopy(c)

        a, b, c, latVecInv, atomPositions, symmetries = getSymmetryOps(a, b, c, latVecInv, atomPositions)

        output = getLatticeDiff(oldA, a, oldB, b, oldC, c, length(symmetries), [])

        oldAtomPos = deepcopy(atomPositions)

        index = 1
        num_symm = Vector{Int8}(undef, length(atomPositions))

        for num in atomTypes
            a, b, c, atomPositions[index:index+num-1], num_symm[index:index+num-1] = symmetrize(a, b, c, latVecInv, atomPositions[index:index+num-1], symmetries, epsilon)
            index += num
        end

        if isDirect
            atomPositions = toDirect(a, b, c, atomPositions)
        end

        output = getAtomDiff(oldAtomPos, atomPositions, num_symm, output)

        writedlm(fout, output)

        if dot(a, cross(b, c)) < 0
            tempa = a
            a = b
            b = tempa
        end

        return a, b, c, atomPositions
    end

    function snapPoscar(in_file = "POSCAR", out_file = "POSCAR.fixed", log_file = "symm_snap.out", epsilon=1e-2)
        fin=readlines(in_file)
        a=[parse(Float64,i) for i in split(fin[3])]
        b=[parse(Float64,i) for i in split(fin[4])]
        c=[parse(Float64,i) for i in split(fin[5])]
        atom_types=[parse(Int8,i) for i in split(fin[6])]
        is_direct = lowercase(fin[7][1]) == "d" ? true : false
        atom_positions=[[parse(Float64,i) for i in split(fin[j])[1:3]] for j in 8:length(fin)]

        a, b, c, atom_positions = SymmetrySnap.symmSnap(a, b, c, atom_positions, atom_types, is_direct, log_file, epsilon)

        fin[3] = join(a, " ")
        fin[4] = join(b, " ")
        fin[5] = join(c, " ")

        writedlm(out_file, fin)
    end
end


text:
"""
    scatter(x, y)

Creates an `EChart` where (x,y) tuples are plotted as dots.

## Methods
```julia
scatter(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}})
scatter(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2})
scatter(df::AbstractDataFrame, x::Symbol, y::Symbol)
scatter(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol)
```

## Arguments
* `mark::Union{String, AbstractVector} = "scatter"` : how to display plotted points
* `legend::Bool` : display legend?
* `scale::Bool = false` : show full Y-axis or truncated
* `large::Bool = true` : minimize overplotting
* `largeThreshold::Int = 2000` : number of points before overplotting optimization occurs
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

## Examples
```julia
sc = scatter(rand(30), rand(30))
```
"""
function scatter(x::AbstractVector{<:Union{Missing, Real}},
				 y::AbstractVector{<:Union{Missing, Real}};
			mark::Union{String, AbstractVector} = "scatter",
			legend::Bool = false,
			scale::Bool = false,
			large::Bool = true,
			largeThreshold::Int = 2000,
			kwargs...)

	ec = xy_plot(x, arrayofarray(x,y); mark = mark, legend = legend, scale = scale, kwargs...)
	ec.xAxis[1]._type = "value"
	ec.xAxis[1].data = nothing #This is necessary to re-use xy_plot, deletes data since scatter uses arrayofarray

	#Enable optimization for lots of data
	[x.large = large for x in ec.series]
	[x.largeThreshold = largeThreshold for x in ec.series]

	return ec

end

#Single X, multiple Y. This is unlikier scenario
function scatter(x::AbstractVector{<:Union{Missing, Real}},
				 y::AbstractArray{<:Union{Missing, Real},2};
			mark::Union{String, AbstractVector} = "scatter",
			legend::Bool = true,
			scale::Bool = false,
			large::Bool = true,
			largeThreshold::Int = 2000,
			kwargs...)

	#Initialize with single series
	ec = scatter(x, y[:,1]; mark = mark, legend = legend, scale = scale, kwargs...)

	#Append remaining series
	for i in 2:size(y)[2]
		push!(ec.series, XYSeries(_type = "scatter", data = arrayofarray(x, y[:,i])))
	end

	seriesnames!(ec)

	legend == true ? legend!(ec) : nothing

	#Enable optimization for lots of data
	[x.large = large for x in ec.series]
	[x.largeThreshold = largeThreshold for x in ec.series]

	return ec

end

#dataframe, single series
function scatter(df::AbstractDataFrame, x::Symbol, y::Symbol;
			mark::Union{String, AbstractVector} = "scatter",
			legend::Bool = false,
			scale::Bool = false,
			large::Bool = true,
			largeThreshold::Int = 2000,
			kwargs...)

	#Initialize with single series
	ec = scatter(df[x], df[y]; mark = mark, legend = legend, scale = scale, kwargs...)

	#Add legend if desired
	legend == true ? legend!(ec) : nothing

	#Name axes since we know them
	xaxis!(ec, name = string(x))
	yaxis!(ec, name = string(y))

	#Enable optimization for lots of data
	[x.large = large for x in ec.series]
	[x.largeThreshold = largeThreshold for x in ec.series]

	return ec

end

#dataframe with group: uses single series df method
function scatter(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol;
			mark::Union{String, AbstractVector} = "scatter",
			legend::Bool = true,
			scale::Bool = false,
			large::Bool = true,
			largeThreshold::Int = 2000,
			kwargs...)

	#Create grouped df
	subdf = groupby(df, group)

	#Get number of groups
	numgroups = length(subdf)

	#Initialize with single series
	ec = scatter(subdf[1], x, y; mark = mark, legend = legend, scale = scale, kwargs...)

	#Append remaining series
	for i in 2:numgroups
		push!(ec.series, XYSeries(_type = "scatter", data = arrayofarray(subdf[i][x], subdf[i][y])))
	end

	#Add series names based on levels of grouped df
	seriesnames!(ec, [unique(subdf[x][group])[1] for x in 1:numgroups])

	#Add legend if desired
	legend == true ? legend!(ec) : nothing

	#Enable optimization for lots of data
	[x.large = large for x in ec.series]
	[x.largeThreshold = largeThreshold for x in ec.series]

	return ec

end


text:
"""
Subtypes of `AbstractRule` are types which represent the primitive derivative
propagation "rules" that can be composed to implement forward- and reverse-mode
automatic differentiation.

More specifically, a `rule::AbstractRule` is a callable Julia object generally
obtained via calling [`frule`](@ref) or [`rrule`](@ref). Such rules accept
differential values as input, evaluate the chain rule using internally stored/
computed partial derivatives to produce a single differential value, then
return that calculated differential value.

For example:

```julia-repl
julia> using ChainRulesCore: frule, rrule, AbstractRule

julia> x, y = rand(2);

julia> h, dh = frule(hypot, x, y);

julia> h == hypot(x, y)
true

julia> isa(dh, AbstractRule)
true

julia> Δx, Δy = rand(2);

julia> dh(Δx, Δy) == ((y / h) * Δx + (x / h) * Δy)
true

julia> h, (dx, dy) = rrule(hypot, x, y);

julia> h == hypot(x, y)
true

julia> isa(dx, AbstractRule) && isa(dy, AbstractRule)
true

julia> Δh = rand();

julia> dx(Δh) == (y / h) * Δh
true

julia> dy(Δh) == (x / h) * Δh
true
```

See also: [`frule`](@ref), [`rrule`](@ref), [`Rule`](@ref), [`DNERule`](@ref), [`WirtingerRule`](@ref)
"""
abstract type AbstractRule end

# this ensures that consumers don't have to special-case rule destructuring
Base.iterate(rule::AbstractRule) = (rule, nothing)
Base.iterate(::AbstractRule, ::Any) = nothing

# This ensures we don't need to check whether the result of `rrule`/`frule` is a tuple
# in order to get the `i`th rule (assuming it's 1)
Base.getindex(rule::AbstractRule, i::Integer) = i == 1 ? rule : throw(BoundsError())

"""
    accumulate(Δ, rule::AbstractRule, args...)

Return `Δ + rule(args...)` evaluated in a manner that supports ChainRulesCore'
various `AbstractDifferential` types.

This method intended to be customizable for specific rules/input types. For
example, here is pseudocode to overload `accumulate` w.r.t. a specific forward
differentiation rule for a given function `f`:

```
df(x) = # forward differentiation primitive implementation

frule(::typeof(f), x) = (f(x), Rule(df))

accumulate(Δ, rule::Rule{typeof(df)}, x) = # customized `accumulate` implementation
```

See also: [`accumulate!`](@ref), [`store!`](@ref), [`AbstractRule`](@ref)
"""
accumulate(Δ, rule::AbstractRule, args...) = add(Δ, rule(args...))

"""
    accumulate!(Δ, rule::AbstractRule, args...)

Similar to [`accumulate`](@ref), but compute `Δ + rule(args...)` in-place,
storing the result in `Δ`.

Note that this function internally calls `Base.Broadcast.materialize!(Δ, ...)`.

See also: [`accumulate`](@ref), [`store!`](@ref), [`AbstractRule`](@ref)
"""
function accumulate!(Δ, rule::AbstractRule, args...)
    return materialize!(Δ, broadcastable(add(cast(Δ), rule(args...))))
end

accumulate!(Δ::Number, rule::AbstractRule, args...) = accumulate(Δ, rule, args...)

"""
    store!(Δ, rule::AbstractRule, args...)

Compute `rule(args...)` and store the result in `Δ`, potentially avoiding
intermediate temporary allocations that might be necessary for alternative
approaches (e.g. `copyto!(Δ, extern(rule(args...)))`)

Note that this function internally calls `Base.Broadcast.materialize!(Δ, ...)`.

Like [`accumulate`](@ref) and [`accumulate!`](@ref), this function is intended
to be customizable for specific rules/input types.

See also: [`accumulate`](@ref), [`accumulate!`](@ref), [`AbstractRule`](@ref)
"""
store!(Δ, rule::AbstractRule, args...) = materialize!(Δ, broadcastable(rule(args...)))

#####
##### `Rule`
#####

Cassette.@context RuleContext

const RULE_CONTEXT = Cassette.disablehooks(RuleContext())

Cassette.overdub(::RuleContext, ::typeof(+), a, b) = add(a, b)
Cassette.overdub(::RuleContext, ::typeof(*), a, b) = mul(a, b)

Cassette.overdub(::RuleContext, ::typeof(add), a, b) = add(a, b)
Cassette.overdub(::RuleContext, ::typeof(mul), a, b) = mul(a, b)

"""
    Rule(propation_function[, updating_function])

Return a `Rule` that wraps the given `propation_function`. It is assumed that
`propation_function` is a callable object whose arguments are differential
values, and whose output is a single differential value calculated by applying
internally stored/computed partial derivatives to the input differential
values.

If an updating function is provided, it is assumed to have the signature `u(Δ, xs...)`
and to store the result of the propagation function applied to the arguments `xs` into
`Δ` in-place, returning `Δ`.

For example:

```
frule(::typeof(*), x, y) = x * y, Rule((Δx, Δy) -> Δx * y + x * Δy)

rrule(::typeof(*), x, y) = x * y, (Rule(ΔΩ -> ΔΩ * y'), Rule(ΔΩ -> x' * ΔΩ))
```

See also: [`frule`](@ref), [`rrule`](@ref), [`accumulate`](@ref), [`accumulate!`](@ref), [`store!`](@ref)
"""
struct Rule{F,U<:Union{Function,Nothing}} <: AbstractRule
    f::F
    u::U
end

# NOTE: Using `Core.Typeof` instead of `typeof` here so that if we define a rule for some
# constructor based on a `UnionAll`, we get `Rule{Type{Thing}}` instead of `Rule{UnionAll}`
Rule(f) = Rule{Core.Typeof(f),Nothing}(f, nothing)

(rule::Rule{F})(args...) where {F} = Cassette.overdub(RULE_CONTEXT, rule.f, args...)

# Specialized accumulation
# TODO: Does this need to be overdubbed in the rule context?
accumulate!(Δ, rule::Rule{F,U}, args...) where {F,U<:Function} = rule.u(Δ, args...)

#####
##### `DNERule`
#####

"""
    DNERule(args...)

Construct a `DNERule` object, which is an `AbstractRule` that signifies that the
current function is not differentiable with respect to a particular parameter.
**DNE** is an abbreviation for Does Not Exist.
"""
struct DNERule <: AbstractRule end

DNERule(args...) = DNE()

#####
##### `WirtingerRule`
#####

"""
TODO
"""
struct WirtingerRule{P<:AbstractRule,C<:AbstractRule} <: AbstractRule
    primal::P
    conjugate::C
end

function (rule::WirtingerRule)(args...)
    return Wirtinger(rule.primal(args...), rule.conjugate(args...))
end

#####
##### `frule`/`rrule`
#####

#=
In some weird ideal sense, the fallback for e.g. `frule` should actually be "get
the derivative via forward-mode AD". This is necessary to enable mixed-mode
rules, where e.g. `frule` is used within a `rrule` definition. For example,
broadcasted functions may not themselves be forward-mode *primitives*, but are
often forward-mode *differentiable*.

ChainRulesCore, by design, is decoupled from any specific AD implementation. How,
then, do we know which AD to fall back to when there isn't a primitive defined?

Well, if you're a greedy AD implementation, you can just overload `frule` and/or
`rrule` to use your AD directly. However, this won't play nice with other AD
packages doing the same thing, and thus could cause load-order-dependent
problems for downstream users.

It turns out, Cassette solves this problem nicely by allowing AD authors to
overload the fallbacks w.r.t. their own context. Example using ForwardDiff:

```
using ChainRulesCore, ForwardDiff, Cassette

Cassette.@context MyChainRuleCtx

# ForwardDiff, itself, can call `my_frule` instead of
# `frule` to utilize the ForwardDiff-injected ChainRulesCore
# infrastructure
my_frule(args...) = Cassette.overdub(MyChainRuleCtx(), frule, args...)

function Cassette.execute(::MyChainRuleCtx, ::typeof(frule), f, x::Number)
    r = frule(f, x)
    if isa(r, Nothing)
        fx, df = (f(x), Rule(Δx -> ForwardDiff.derivative(f, x) * Δx))
    else
        fx, df = r
    end
    return fx, df
end
```
=#

"""
    frule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (rule_for_ΔΩ₁::AbstractRule, rule_for_ΔΩ₂::AbstractRule, ...))

where each returned propagation rule `rule_for_ΔΩᵢ` can be invoked as

    rule_for_ΔΩᵢ(Δx₁, Δx₂, ...)

to yield `Ωᵢ`'s corresponding differential `ΔΩᵢ`. To illustrate, if all involved
values are real-valued scalars, this differential can be written as:

    ΔΩᵢ = ∂Ωᵢ_∂x₁ * Δx₁ + ∂Ωᵢ_∂x₂ * Δx₂ + ...

If no method matching `frule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```julia-repl
julia> x = rand();

julia> sinx, dsin = frule(sin, x);

julia> sinx == sin(x)
true

julia> dsin(1) == cos(x)
true
```

unary input, binary output scalar function:

```julia-repl
julia> x = rand();

julia> sincosx, (dsin, dcos) = frule(sincos, x);

julia> sincosx == sincos(x)
true

julia> dsin(1) == cos(x)
true

julia> dcos(1) == -sin(x)
true
```

See also: [`rrule`](@ref), [`AbstractRule`](@ref), [`@scalar_rule`](@ref)
"""
frule(::Any, ::Vararg{Any}; kwargs...) = nothing

"""
    rrule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (rule_for_Δx₁::AbstractRule, rule_for_Δx₂::AbstractRule, ...))

where each returned propagation rule `rule_for_Δxᵢ` can be invoked as

    rule_for_Δxᵢ(ΔΩ₁, ΔΩ₂, ...)

to yield `xᵢ`'s corresponding differential `Δxᵢ`. To illustrate, if all involved
values are real-valued scalars, this differential can be written as:

    Δxᵢ = ∂Ω₁_∂xᵢ * ΔΩ₁ + ∂Ω₂_∂xᵢ * ΔΩ₂ + ...

If no method matching `rrule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```julia-repl
julia> x = rand();

julia> sinx, dx = rrule(sin, x);

julia> sinx == sin(x)
true

julia> dx(1) == cos(x)
true
```

binary input, unary output scalar function:

```julia-repl
julia> x, y = rand(2);

julia> hypotxy, (dx, dy) = rrule(hypot, x, y);

julia> hypotxy == hypot(x, y)
true

julia> dx(1) == (y / hypot(x, y))
true

julia> dy(1) == (x / hypot(x, y))
true
```

See also: [`frule`](@ref), [`AbstractRule`](@ref), [`@scalar_rule`](@ref)
"""
rrule(::Any, ::Vararg{Any}; kwargs...) = nothing


text:
module BibTeX

using DataStructures
using BibInternal

export parse_string

"""
    mutable struct Accumulator

A structure to accumulate part of the input stream as a simple string. It differs from Position that store `line` and `column` information.
"""
mutable struct Accumulator
    from::Int
    to::Int
end

"""
    Content

Store the different BibTeX elements once succesfully parsed.

!!! warning "Note:"
    Free text, comments entries and preambles entries are currently ignored.
"""
struct Content
    # comments::Dict{Int, String}
    entries::OrderedDict{String,BibInternal.Entry}
    # free::Dict{Int, String}
    # preambles::Dict{Int, String}
    strings::Dict{String,String}
end

"""
    Content()

Create an (almost) empty content constructor for a new `Parser`. Some usual BibTeX content is added:
- List of months abbreviations in English (added automatically by some BibTeX class)

Feel free to contribute to more default content.
"""
function Content()
    # comments = Dict{Int, String}()
    entries = OrderedDict{String,BibInternal.Entry}()
    # free = Dict{Int, String}()
    # preambles = Dict{Int, String}()
    strings = Dict{String,String}([
        "jan" => "January",
        "feb" => "February",
        "mar" => "March",
        "apr" => "April",
        "may" => "May",
        "jun" => "June",
        "jul" => "July",
        "aug" => "August",
        "sep" => "September",
        "oct" => "October",
        "nov" => "November",
        "dec" => "December",
    ])
    return Content(entries, strings)
end

"""
    Field

A structure that can store BibTeX fields information.

# Arguments:
- `braces::Int`: counter to the number of opened braces in a field being parsed
- `name::String`: name of the field being parsed
- `quotes::Bool`: `true` if the field is delimited by quotes
- `value::String`: the value of the field
"""
mutable struct Field
    braces::Int
    name::String
    quotes::Bool
    value::String

    Field() = new(0, "", false, "")
end

"""
    Position

A structure pointing to a position in the input based on rows and columns.
"""
mutable struct Position
    row::Int
    col::Int
end

"""
    Position()

Initial position constructor
"""
Position() = Position(0, 0)

"""
    Storage

Store the content of an entry being parsed.

# Arguments:
- `delim::Union{Char, Nothing}`: the character delimiting the entry
- `fields::Vector{Field}`: a collection of the entry's fields
- `key::String`: the key of the entry
- `kind::String`: the kind/type of BibTeX entry
"""
mutable struct Storage
    delim::Union{Char,Nothing}
    fields::Vector{Field}
    key::String
    kind::String
end

"""
    Storage()
Empty storage constructor called when a new entry is being parsed.
"""
Storage() = Storage(nothing, Vector{Field}(), "", "")

"""
    make_entry(storage)

Make a `BibInternal.Entry` from a completed entry in a parser storage.
"""
function make_entry(storage)
    d = Dict("_type" => storage.kind)
    foreach(field -> push!(d, field.name => field.value), storage.fields)
    return d
end

"""
    BibTeXError

Description of a BibTeXError, including the relevant position in the input string.

# Arguments:
- `kind::Symbol`: the type of BibTeX error
- `input::String`: the BibTeX string being parsed
- `start::Position`: the row/col start position of the BibTeXError
- `stop::Position`: the row/col end position of the BibTeXError
"""
struct BibTeXError
    kind::Symbol
    input::String
    start::Position
    stop::Position

    BibTeXError(k, i, sta, sto) = new(k, "'" * strip(i) * "'", deepcopy(sta), deepcopy(sto))
end

function warn(error, ::Val)
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The entry kind is invalid from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

function warn(error, ::Val{:invalid_string})
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The string entry is invalid from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

function warn(error, ::Val{:incomplete_entry})
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The entry is incomplete and end from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

function warn(error, ::Val{:invalid_key})
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The entry key is invalid from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

function warn(error, ::Val{:invalid_field_name})
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The field name is invalid from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

function warn(error, ::Val{:invalid_field_number})
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The field value has an invalid format (number) from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

function warn(error, ::Val{:invalid_field_var})
    row_start, col_start = error.start.row, error.start.col
    row_stop, col_stop = error.stop.row, error.stop.col
    return str = "The field value has an invalid format (string variable) from (line $row_start, character $col_start) to (line $row_stop, character $col_stop): $(error.input)"
end

"""
    warn(error)

Dispatch a BibTeX error as a Julia warning based on the type of error.
"""
warn(error) = @warn warn(error, Val(error.kind))

"""
    Parser

A structure allowing to parse a BibTeX formatted string one character at a time.

# Arguments:
- `acc::Accumulator`: an accumulated string from last `dump` to `content`
- `content::Content`: the current content of a parser
- `errors::Vector{BibTeXError}`: a collection of BibTeX errors
- `field::Field`: temporary storage for an entry fields names
- `input::Vector{Char}`: the BibTeX string
- `pos_start::Position`: pointer to the raw/col start position
- `pos_end::Position`: pointer to the raw/col end position
- `storage::Storage`: temporary storage of the content of an entry being parsed
- `rules_checker::Bool`: indicate which level `:error`, `:warn`, or `:none` should be raised when an entry violates BibTeX rules
- `task::Symbol`: describe which part of the BibTeX gramma is being parsed
"""
mutable struct Parser
    acc::Accumulator
    content::Content
    errors::Vector{BibTeXError}
    field::Field
    input::Vector{Char}
    pos_start::Position
    pos_end::Position
    rules_checker::Symbol
    storage::Storage
    task::Symbol

    function Parser(
        input;
        acc=Accumulator(1, 0),
        content=Content(),
        errors=Vector{BibTeXError}(),
        field=Field(),
        pos_start=Position(1, 1),
        pos_end=Position(1, 0),
        rules_checker=:error,
        storage=Storage(),
        task=:free,
    )
        return new(
            acc,
            content,
            errors,
            field,
            collect(input),
            pos_start,
            pos_end,
            rules_checker,
            storage,
            task,
        )
    end
end

"""
    rev(char)

Return the closing character `)`/`}` matching either `(` or `{`.
"""
rev(char) = char == '(' ? ')' : '}'

"""
    get_entries(parser)

Retrieve the entries succesfully parsed by `parser`.
"""
get_entries(parser) = parser.content.entries

"""
    get_acc(parser; from = 1, to = 0)

Retrieve the `Accumulator` of the parser.

# Arguments:
- `parser`: a `Parser`
- `from`: a positive offset from the start of the `Accumulator`, default to `1`
- `to`:  a negative offset from the end of the `Accumulator`, default to `0`
"""
function get_acc(parser; from=1, to=0)
    a = from + parser.acc.from - 1
    b = parser.acc.to - to
    return prod(parser.input[a:b])
end

"""
    set_delim!(parser, char)

Set the delimiter for this section of the entry to parenthesis or braces.
"""
set_delim!(parser, char) = parser.storage.delim = char

"""
    set_entry_kind!(parser, kind) = begin

Set the kind of the entry being parsed in the `parser`'s storage.
"""
set_entry_kind!(parser, kind) = parser.storage.kind = kind

"""
    inc_col!(parser) = begin

Increment the `column` field of the pointer to the end position of the Parser.
"""
inc_col!(parser) = parser.pos_end.col += 1

"""
    inc_row!(parser)

Increment the `row` field of the pointer to the end position of the Parser.
"""
function inc_row!(parser)
    parser.pos_end.row += 1
    return parser.pos_end.col = 1
end

"""
    inc!(parser, char, dumped)

Increment the start/end position and accumulator of the parser.

# Arguments:
- `parser`: a `BibTeX.Parser` structure
- `char`: the character being parsed
- `dumped`: boolean value describing if the `Accumulator` content is being dumped into the parser `Content`
"""
function inc!(parser, char, dumped)
    char == '\n' ? inc_row!(parser) : inc_col!(parser)
    parser.acc.to += 1
    if dumped
        parser.pos_start = deepcopy(parser.pos_end)
        parser.acc.from = parser.acc.to
    end
end

is_dumped(::Parser, char, ::Val{:free}) = char == '@'
dump!(parser, char, ::Val) = char == '@' && (parser.task = :entry)

is_dumped(::Parser, char, ::Val{:entry}) = occursin(r"[@{\(\n]", char)
function dump!(parser, char, ::Val{:entry})
    if char == '\n'
        parser.task = :free
        e = BibTeXError(:invalid_kind, get_acc(parser), parser.pos_start, parser.pos_end)
        push!(parser.errors, e)
        parser.storage = Storage()
    elseif char ∈ ['{', '(']
        set_delim!(parser, char)
        acc = split(lowercase(get_acc(parser; from=2)), r"[\t ]+")
        if length(acc) ≤ 1 && !isempty(acc[1])
            set_entry_kind!(parser, acc[1])
            # parser.task = acc[1] ∈ ["comment", "preamble", "string"] ? Symbol(acc[1]) : :key
            parser.task = acc[1] ∈ ["string"] ? Symbol(acc[1]) : :key
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_kind, get_acc(parser), parser.pos_start, parser.pos_end
            )
            push!(parser.errors, e)
            parser.storage = Storage()
        end
    end
end

is_dumped(::Parser, char, ::Val{:key}) = char ∈ ['@', ',']
function dump!(parser, char, ::Val{:key})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == ','
        acc = split(get_acc(parser; from=2), r"[\t ]+")
        if length(acc) ≤ 1 && !isempty(acc[1])
            parser.storage.key = acc[1]
            parser.task = :field_name
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_key, get_acc(parser; from=2), parser.pos_start, parser.pos_end
            )
            push!(parser.errors, e)
            parser.storage = Storage()
        end
    end
end

is_dumped(::Parser, char, ::Val{:field_name}) = char ∈ ['=', '@']
function dump!(parser, char, ::Val{:field_name})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '='
        acc = split(get_acc(parser; from=2), r"[\t\r\n ]+"; keepempty=false)
        if length(acc) == 1
            parser.field.name = acc[1]
            parser.task = :field_in
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_field_name,
                get_acc(parser; from=2),
                parser.pos_start,
                parser.pos_end,
            )
            push!(parser.errors, e)
            parser.storage = Storage()
        end
    end
end

is_dumped(::Parser, char, ::Val{:field_in}) = occursin(r"[0-9@a-zA-Z\"{]", char)
function dump!(parser, char, ::Val{:field_in})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '"'
        parser.pos_start.col += 1
        parser.acc.from += 1
        parser.task = :field_inquote
    elseif char == '{'
        parser.task = :field_inbrace
        parser.field.braces += 1
    elseif occursin(r"[0-9]", char)
        parser.task = :field_number
    elseif occursin(r"[a-zA-Z]", char)
        parser.task = :field_var
    end
end

function is_dumped(parser, char, ::Val{:field_inbrace})
    if char == '{'
        parser.field.braces += 1
        return false
    elseif char == '}'
        parser.field.braces -= 1
        return (parser.field.braces == 0)
    elseif char == '@'
        return true
    else
        return false
    end
end
function dump!(parser, char, ::Val{:field_inbrace})
    if char == '@'
        parser.task = :entry
        parser.field.braces = 0
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '}'
        parser.field.value *= get_acc(parser; from=2)
        push!(parser.storage.fields, deepcopy(parser.field))
        parser.field.value = ""
        parser.task = :field_out
    end
end

is_dumped(::Parser, char, ::Val{:field_inquote}) = char ∈ ['{', '"', '}']
function dump!(parser, char, ::Val{:field_inquote})
    if char == '"' && parser.field.braces == 0
        parser.field.value *= get_acc(parser; from=2)
        parser.task = :field_outquote
    elseif char == '{'
        parser.field.value *= get_acc(parser; from=2) # * parser.field.value
        parser.field.braces += 1
    elseif char == '}'
        parser.field.value *= get_acc(parser; from=1, to=-1) # * parser.field.value
        parser.field.braces -= 1
    end
end

function is_dumped(parser, char, ::Val{:field_outquote})
    return char ∈ ['@', ',', '#', rev(parser.storage.delim)]
end
function dump!(parser, char, ::Val{:field_outquote})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '#'
        parser.task = :field_concat
    else
        push!(parser.storage.fields, deepcopy(parser.field))
        parser.field.value = ""
        if char == ','
            parser.task = :field_next
        elseif char == rev(parser.storage.delim)
            entry = make_entry(parser.storage)
            check = parser.rules_checker
            push!(
                parser.content.entries,
                parser.storage.key =>
                    BibInternal.make_bibtex_entry(parser.storage.key, entry; check),
            )
            parser.storage = Storage()
            parser.task = :free
        end
    end
end

is_dumped(::Parser, char, ::Val{:field_concat}) = occursin(r"[a-zA-Z\"@]", char)
function dump!(parser, char, ::Val{:field_concat})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '"'
        parser.pos_start.col += 1
        parser.acc.from += 1
        parser.task = :field_inquote
    elseif occursin(r"[a-zA-Z]", char)
        parser.task = :field_var
    end
end

function is_dumped(parser, char, ::Val{:field_var})
    return char ∈ ['@', ',', '#', rev(parser.storage.delim)]
end
function dump!(parser, char, ::Val{:field_var})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    else
        acc = split(get_acc(parser), r"[\t\r\n ]+"; keepempty=false)
        if length(acc) == 1
            parser.field.value *= parser.content.strings[acc[1]]
            if char == '#'
                parser.task = :field_concat
            else
                push!(parser.storage.fields, deepcopy(parser.field))
                parser.field.value = ""
                if char == ','
                    parser.task = :field_next
                elseif char == rev(parser.storage.delim)
                    entry = make_entry(parser.storage)
                    check = parser.rules_checker
                    push!(
                        parser.content.entries,
                        parser.storage.key =>
                            BibInternal.make_bibtex_entry(parser.storage.key, entry; check),
                    )
                    parser.storage = Storage()
                    parser.task = :free
                end
            end
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_field_var, get_acc(parser), parser.pos_start, parser.pos_end
            )
            push!(parser.errors, e)
            parser.storage = Storage()
        end
    end
end

is_dumped(parser, char, ::Val{:field_number}) = char ∈ ['@', ',', rev(parser.storage.delim)]
function dump!(parser, char, ::Val{:field_number})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    else
        acc = split(get_acc(parser), r"[\t\r\n ]+"; keepempty=false)
        # @show acc
        if length(acc) == 1
            parser.field.value = acc[1]
            push!(parser.storage.fields, deepcopy(parser.field))
            parser.field.value = ""
            if char == ','
                parser.task = :field_next
            elseif char == rev(parser.storage.delim)
                entry = make_entry(parser.storage)
                check = parser.rules_checker
                push!(
                    parser.content.entries,
                    parser.storage.key =>
                        BibInternal.make_bibtex_entry(parser.storage.key, entry; check),
                )
                parser.task = :free
            end
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_field_number, get_acc(parser), parser.pos_start, parser.pos_end
            )
            push!(parser.errors, e)
        end
    end
end

is_dumped(parser, char, ::Val{:field_out}) = char ∈ ['@', ',', rev(parser.storage.delim)]
function dump!(parser, char, ::Val{:field_out})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == ','
        parser.task = :field_next
    elseif char == rev(parser.storage.delim)
        entry = make_entry(parser.storage)
        check = parser.rules_checker
        push!(
            parser.content.entries,
            parser.storage.key => BibInternal.make_bibtex_entry(parser.storage.key, entry; check),
        )
        parser.storage = Storage()
        parser.task = :free
    end
end

is_dumped(parser, char, ::Val{:field_next}) = char ∈ ['=', '@', rev(parser.storage.delim)]
function dump!(parser, char, ::Val{:field_next})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '='
        acc = split(get_acc(parser; from=2), r"[\t\r\n ]+"; keepempty=false)
        # @show acc
        if length(acc) == 1
            parser.field.name = acc[1]
            parser.task = :field_in
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_field_name, get_acc(parser), parser.pos_start, parser.pos_end
            )
            push!(parser.errors, e)
            parser.storage = Storage()
        end
    elseif char == rev(parser.storage.delim)
        entry = make_entry(parser.storage)
        check = parser.rules_checker
        push!(
            parser.content.entries,
            parser.storage.key => BibInternal.make_bibtex_entry(parser.storage.key, entry; check),
        )
        parser.storage = Storage()
        parser.task = :free
    end
end

is_dumped(::Parser, char, ::Val{:string}) = char ∈ ['=', '@']
function dump!(parser, char, ::Val{:string})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '='
        acc = split(get_acc(parser; from=2), r"[\t\r\n ]+"; keepempty=false)
        if length(acc) == 1
            parser.field.name = acc[1]
            parser.task = :string_inquote
        else
            parser.task = :free
            e = BibTeXError(
                :invalid_string,
                get_acc(parser; from=2, to=-1),
                parser.pos_start,
                parser.pos_end,
            )
            push!(parser.errors, e)
            parser.storage = Storage()
        end
    end
end

is_dumped(::Parser, char, ::Val{:string_inquote}) = char ∈ ['"', '@']
function dump!(parser, char, ::Val{:string_inquote})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == '"'
        parser.task = :string_value
    end
end

is_dumped(::Parser, char, ::Val{:string_value}) = char ∈ ['"']
function dump!(parser, char, ::Val{:string_value})
    if char == '"'
        parser.field.value = get_acc(parser; from=2)
        parser.task = :string_outquote
    end
end

is_dumped(parser, char, ::Val{:string_outquote}) = char ∈ [rev(parser.storage.delim), '@']
function dump!(parser, char, ::Val{:string_outquote})
    if char == '@'
        parser.task = :entry
        e = BibTeXError(
            :incomplete_entry, get_acc(parser), parser.pos_start, parser.pos_end
        )
        push!(parser.errors, e)
    elseif char == rev(parser.storage.delim)
        parser.content.strings[parser.field.name] = parser.field.value
        parser.task = :free
        parser.storage = Storage()
    end
end

"""
    is_dumped(parser, char)

Check if an `Accumulator` needs to be `dump!`ed. Dispatch to the appropriate `is_dumped` method according to the state of the parser.
"""
is_dumped(parser, char) = is_dumped(parser, char, Val(parser.task))

"""
    dump!(parser, char = ' ')

Dump the content of the parser `Accumulator` into the parser `Content`. Dispatch to the appropriate `dump!` method according to the state of the parser.
"""
function dump!(parser, char=' ')
    dump!(parser, char, Val(parser.task))
    parser.pos_start = deepcopy(parser.pos_end)
    return parser.acc.from = parser.acc.to
end

"""
    parse!(parser, char)

Parse a single character of a BibTeX string. Modify the `Parser` accordingly.
"""
function parse!(parser, char)
    dumped = is_dumped(parser, char)
    dumped && dump!(parser, char)
    return inc!(parser, char, dumped)
end

"""
    parse_string(str)

Parse a BibTeX string of entries. Raise a detailed warning for each invalid entry.
"""
function parse_string(str; check=:error)
    parser = Parser(str; rules_checker=check)
    foreach(char -> parse!(parser, char), parser.input)
    foreach(error -> warn(error), parser.errors)
    return get_entries(parser)
end

"""
    parse_file(path)

Parse a BibTeX file located at `path`. Raise a detailed warning for each invalid entry.
"""
parse_file(path; check=:error) = parse_string(read(path, String); check)

end # module


text:
function template()
    tpl = mt"""#!/bin/sh
     # embedded options to qsub - start with #PBS
     # -- Name of the job ---
     #BSUB -J AmrajJob
     # �~@~S- specify queue --
     #BSUB -q {{cluster}}
     ### -- ask for number of cores (default: 1) --
     #BSUB -n {{cores}}
     ### -- specify that the cores must be on the same host --
     #BSUB -R "span[hosts=1]"
     ### -- specify that we need 2GB of memory per core/slot --
     #BSUB -R "rusage[mem=8GB]"
     ### -- specify that we want the job to get killed if it exceeds 3 GB per core/slot --
     #BSUB -M 8GB
     ### -- set walltime limit: hh:mm --
     #BSUB -W 8:00
     ### -- set the email address --
     # please uncomment the following line and put in your e-mail address,
     # if you want to receive e-mail notifications on a non-default address
     #BSUB -u amraj@dtu.dk
     ### -- send notification at start --
     ##BSUB -B
     ### -- send notification at completion --
     ##BSUB -N
     ### -- Specify the output and error file. %J is the job-id --
     ### -- -o and -e mean append, -oo and -eo mean overwrite --
     #BSUB -o {{{workspace}}}/Output.out
     #BSUB -e {{{workspace}}}/Error.err

     export JULIA_NUM_THREADS={{cores}}
     export SERVER_ID=$SERVER_ID

     {{{command}}}
     """
    return tpl
end

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
# closest points using a particular metric, using a K-D tree implemented by NearestNeighbors
function closest_points(P,Q; metric=Euclidean())
    kdtree = KDTree(Q, metric)
    nearestidx, dists = nn(kdtree, P)
    return [(i,nearestidx[i]) for i=1:size(P,2)]
end

# generate matrices for Kabsch from a list of correspondences
matched_points(P::AbstractMatrix, Q::AbstractMatrix, matches) = (hcat([P[:,i] for (i,j) in matches]...), hcat([Q[:,j] for (i,j) in matches]...))
matched_points(P, Q, matches) = ([P[i] for (i,j) in matches], [Q[j] for (i,j) in matches])
matched_points(P,Q; correspondence=closest_points, kwargs...) = matched_points(P, Q, correspondence(P,Q; kwargs...))


text:
"""
    Revise.WatchList

A struct for holding files that live inside a directory.
Some platforms (OSX) have trouble watching too many files. So we
watch parent directories, and keep track of which files in them
should be tracked.

Fields:
- `timestamp`: mtime of last update
- `trackedfiles`: Set of filenames, generally expressed as a relative path
"""
mutable struct WatchList
    timestamp::Float64         # unix time of last revision
    trackedfiles::Set{String}
end

const DefMapValue = Union{RelocatableExpr,Tuple{Vector{Any},Int}}  # ([sigt1,...], lineoffset)

"""
    DefMap

Maps `def`=>nothing or `def=>([sigt1,...], lineoffset)`, where:

- `def` is an expression
- the value is `nothing` if `def` does not define a method
- if it does define a method, `sigt1...` are the signature-types and `lineoffset` is the
  difference between the line number when the method was compiled and the current state
  of the source file.

See the documentation page [How Revise works](@ref) for more information.
"""
const DefMap = OrderedDict{RelocatableExpr,Union{DefMapValue,Nothing}}

"""
    SigtMap

Maps `sigt=>def`, where `sigt` is the signature-type of a method and `def` the expression
defining the method.

See the documentation page [How Revise works](@ref) for more information.
"""
const SigtMap = IdDict{Any,RelocatableExpr}   # sigt=>def

"""
    FMMaps

`defmap` (`source=>sigtypes`) and `sigtmap` (`sigtypes=>source`) mappings for a
particular file/module combination.
See the documentation page [How Revise works](@ref) for more information.
"""
struct FMMaps
    defmap::DefMap
    sigtmap::SigtMap
end
FMMaps() = FMMaps(DefMap(), SigtMap())

Base.isempty(fmm::FMMaps) = isempty(fmm.defmap)

function Base.show(io::IO, fmm::FMMaps)
    limit = get(io, :limit, true)
    if limit
        print(io, "FMMaps(<$(length(fmm.defmap)) expressions>, <$(length(fmm.sigtmap)) signatures>)")
    else
        println(io, "FMMaps with the following expressions:")
        for def in keys(fmm.defmap)
            print(io, "  ")
            Base.show_unquoted(io, def, 2)
            print(io, '\n')
        end
    end
end

"""
    FileModules

For a particular source file, the corresponding `FileModules` is an
`OrderedDict(mod1=>fmm1, mod2=>fmm2)`,
mapping the collection of modules "active" in the file (the parent module and any
submodules it defines) to their corresponding [`FMMaps`](@ref).

The first key is guaranteed to be the module into which this file was `include`d.

To create a `FileModules` from a source file, see [`parse_source`](@ref).
"""
const FileModules = OrderedDict{Module,FMMaps}

"""
    fm = FileModules(mod::Module)

Initialize an empty `FileModules` for a file that is `include`d into `mod`.
"""
FileModules(mod::Module) = FileModules(mod=>FMMaps())

Base.isempty(fm::FileModules) = length(fm) == 1 && isempty(first(values(fm)))

"""
    FileInfo(fm::FileModules, cachefile="")

Structure to hold the per-module expressions found when parsing a
single file.
`fm` holds the [`FileModules`](@ref) for the file.

Optionally, a `FileInfo` can also record the path to a cache file holding the original source code.
This is applicable only for precompiled modules and `Base`.
(This cache file is distinct from the original source file that might be edited by the
developer, and it will always hold the state
of the code when the package was precompiled or Julia's `Base` was built.)
When a cache is available, `fm` will be empty until the file gets edited:
the original source code gets parsed only when a revision needs to be made.

Source cache files greatly reduce the overhead of using Revise.
"""
struct FileInfo
    fm::FileModules
    cachefile::String
end
FileInfo(fm::FileModules) = FileInfo(fm, "")

"""
    FileInfo(mod::Module, cachefile="")

Initialze an empty FileInfo for a file that is `include`d into `mod`.
"""
FileInfo(mod::Module, cachefile::AbstractString="") = FileInfo(FileModules(mod), cachefile)

FileInfo(fm::FileModules, fi::FileInfo) = FileInfo(fm, fi.cachefile)

"""
    PkgData(id, path, fileinfos::Dict{String,FileInfo})

A structure holding the data required to handle a particular package.
`path` is the top-level directory defining the package,
and `fileinfos` holds the [`FileInfo`](@ref) for each file defining the package.

For the `PkgData` associated with `Main` (e.g., for files loaded with [`includet`](@ref)),
the corresponding `path` entry will be empty.
"""
mutable struct PkgData
    id::PkgId
    path::String
    fileinfos::Dict{String,FileInfo}
end

PkgData(id::PkgId, path) = PkgData(id, path, Dict{String,FileInfo}())
PkgData(id::PkgId, ::Nothing) = PkgData(id, "")
PkgData(id::PkgId) = PkgData(id, normpath(basepath(id)))

struct GitRepoException <: Exception
    filename::String
end

function Base.showerror(io::IO, ex::GitRepoException)
    print(io, "no repository at ", ex.filename, " to track stdlibs you must build Julia from source")
end

"""
    Rescheduler(f, args)

To facilitate precompilation and reduce latency, we replace

```julia
function watch_manifest(mfile)
    wait_changed(mfile)
    # stuff
    @async watch_manifest(mfile)
end

@async watch_manifest(mfile)
```

with a rescheduling type:

```julia
fresched = Rescheduler(watch_manifest, (mfile,))
schedule(Task(fresched))
```

where now `watch_manifest(mfile)` should return `true` if the task
should be rescheduled after completion, and `false` otherwise.
"""
struct Rescheduler{F,A}
    f::F
    args::A
end

function (thunk::Rescheduler{F,A})() where {F,A}
    if thunk.f(thunk.args...)::Bool
        schedule(Task(thunk))
    end
end


text:
using BasisFunctions: hasderivative, hasantiderivative, resize

P = 80

b = PBSpline(2, 5)

@testset "$(rpad("Periodic BSplines",P))" begin

    BasisFunctions.Test.test_generic_dict_interface(PBSpline(3, 20))

    @testset "$(rpad("Basic functionality",P))" begin
        ξ = nodes(b)
        @test ξ == MoreBasisFunctions.get_pbspline_nodes(Float64, 5)
        @test nnodes(b) == 5
        @test degree(b) == 2
        @test support(b) == PBSplineInterval{eltype(ξ)}()
        @test nodes(b) == nodes(similar(b, eltype(ξ), 5))
        @test nodes(b) == nodes(similar(b, eltype(ξ), ξ))
        @test nodes(resize(b, 5)) == ξ
        @test nodes(resize(b, 3)) == [0.0, 0.5, 1.0]
    end

    @testset "$(rpad("Basis evaluation",P))" begin
        @test b[1](0.2) == 0.5
        @test b[1](0.3) == 0.75
        @test b[1](0.6) == 0.0

        @test b[5](0.0) == 0.5
        @test b[5](1.0) == 0.5
        @test b[5](0.9) == 0.125
    end

    @testset "$(rpad("Basis derivative",P))" begin
        @test eval_element_derivative(b, 1, 0.1) == 2.5
        @test eval_element_derivative(b, 1, 0.3) == 0.0
    end

end


text:
export cart_to_mat, mat_to_vec, remove_background, gaussf

"""
    cart_to_mat(cart)

converts a list of cartesian indices to a matrix with the first dimension corresponding to the dimension direction.

"""
function cart_to_mat(cart)
    all_max_vec = zeros(length(cart[1]), length(cart))
    n=1
    for idx in cart
        all_max_vec[:,n] = [Tuple(idx)...]
        n += 1
    end
    return all_max_vec
end

"""
    mat_to_vec(mat)

converts a 2D matrix to a vector of vectors, one vector for each column.

"""
function mat_to_vec(mat)
    [mat[:,n] for n=1:size(mat,2)]
end

"""
    gaussf(img::AbstractArray, sigma=1.0; dtype=Float32)
performs Gaussian filtering using FFTs. `sigma` denotes the kernel sizes along the various dimensions.
#arguments
+ `img`: the image to apply the Gaussian filter to
+ ``
"""
function gaussf(img::AbstractArray{T}, sigma=1.0; dtype=Float32) where (T <: Real)
    dims=1:ndims(img)
    shiftdims = dims[2:end]
    img = dtype.(img)
    f = rfft(img, dims)
    return irfft(f .* ifftshift(gaussian(eltype(img), size(f), offset=CtrRFT, sigma=size(img) ./(2π .*sigma)), shiftdims), size(img, dims[1]), dims)
end

function gaussf(img::AbstractArray, sigma=1.0)
    return (ifft(fft(img) .* ifftshift(gaussian(eltype(img), size(img), sigma=size(img) ./(2π*sigma)))))
end

"""
    remove_background(img, sigma=5.0)

estimates the background by applying a Gaussian filter to the data and measuring for the minimum of this filtered data.
#arguments
+ `img`: the image to remove the background from
+ `sigma`: the kernelsize for the Gaussian filter
#returns
a tuple of the background-subtracted image and the background value.
Note that no zero-clipping is performed. The result type is by defaul Float32, if the image itself is not FLoat64 or Complex.
"""
function remove_background(img, sigma=5.0; dtype=Float32)
    bg = minimum(gaussf(img,sigma; dtype=dtype))
    return img .- bg, bg
end


text:
using Test

include("common.jl")
include("byte_stream_test_harness.jl")

for test_file in readdir(@__DIR__)
    if test_file != "runtests.jl" && endswith(test_file, ".jl")
        include(test_file)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge