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

@doc raw"""Adapts a Secret into a volume.  The contents of the target Secret's Data field will be presented in a volume as files using the keys in the Data field as the file names. Secret volumes support ownership management and SELinux relabeling.

    IoK8sApiCoreV1SecretVolumeSource(;
        defaultMode=nothing,
        items=nothing,
        optional=nothing,
        secretName=nothing,
    )

    - defaultMode::Int32 : Optional: mode bits to use on created files by default. Must be a value between 0 and 0777. Defaults to 0644. Directories within the path are not affected by this setting. This might be in conflict with other options that affect the file mode, like fsGroup, and the result can be other mode bits set.
    - items::Vector{IoK8sApiCoreV1KeyToPath} : If unspecified, each key-value pair in the Data field of the referenced Secret will be projected into the volume as a file whose name is the key and content is the value. If specified, the listed keys will be projected into the specified paths, and unlisted keys will not be present. If a key is specified which is not present in the Secret, the volume setup will error unless it is marked optional. Paths must be relative and may not contain the '..' path or start with '..'.
    - optional::Bool : Specify whether the Secret or its keys must be defined
    - secretName::String : Name of the secret in the pod's namespace to use. More info: https://kubernetes.io/docs/concepts/storage/volumes#secret
"""
mutable struct IoK8sApiCoreV1SecretVolumeSource <: SwaggerModel
    defaultMode::Any # spec type: Union{ Nothing, Int32 } # spec name: defaultMode
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1KeyToPath} } # spec name: items
    optional::Any # spec type: Union{ Nothing, Bool } # spec name: optional
    secretName::Any # spec type: Union{ Nothing, String } # spec name: secretName

    function IoK8sApiCoreV1SecretVolumeSource(;defaultMode=nothing, items=nothing, optional=nothing, secretName=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecretVolumeSource, Symbol("defaultMode"), defaultMode)
        setfield!(o, Symbol("defaultMode"), defaultMode)
        validate_property(IoK8sApiCoreV1SecretVolumeSource, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiCoreV1SecretVolumeSource, Symbol("optional"), optional)
        setfield!(o, Symbol("optional"), optional)
        validate_property(IoK8sApiCoreV1SecretVolumeSource, Symbol("secretName"), secretName)
        setfield!(o, Symbol("secretName"), secretName)
        o
    end
end # type IoK8sApiCoreV1SecretVolumeSource

const _property_map_IoK8sApiCoreV1SecretVolumeSource = Dict{Symbol,Symbol}(Symbol("defaultMode")=>Symbol("defaultMode"), Symbol("items")=>Symbol("items"), Symbol("optional")=>Symbol("optional"), Symbol("secretName")=>Symbol("secretName"))
const _property_types_IoK8sApiCoreV1SecretVolumeSource = Dict{Symbol,String}(Symbol("defaultMode")=>"Int32", Symbol("items")=>"Vector{IoK8sApiCoreV1KeyToPath}", Symbol("optional")=>"Bool", Symbol("secretName")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1SecretVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1SecretVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecretVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecretVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecretVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecretVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1SecretVolumeSource)
    true
end

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


text:
# Discrete entropy estimators

###########################################################
#
#   Entropy estimator types
#
###########################################################

type MaximumLikelihoodEntropyEstimator <: EntropyEstimator end
typealias MLEntropyEstimator MaximumLikelihoodEntropyEstimator
typealias MLEE MaximumLikelihoodEntropyEstimator

type MillerMaddowBiasCorrectedEntropyEstimator <: EntropyEstimator end
typealias MillerMaddowBiasCorrected MillerMaddowBiasCorrectedEntropyEstimator
typealias MMBCEE MillerMaddowBiasCorrectedEntropyEstimator

###########################################################
#
#   Entropy estimator implementations
#
###########################################################

# Maximum Likelihood

function entropy(e::MLEE, f::Frequencies)
    p = f.f / sum(f.f)
    p::Vector{Float64}
    StatsBase.entropy(p)
end

entropy_ml(o::Observations) = entropy(MLEE(), o)
entropy_ml(o::Observations, base::Real) = entropy(MLEE(), o, base)

# Millor-Maddow Corrected

function entropy{T}(e::MillerMaddowBiasCorrected, f::Frequencies{T})
    n = sum(f.f) # number of samples
    nz = sum(f.f .!= zero(T)) # number of non-zero categories
    entropy(MLEE(), f) + (nz - 1) / (2 * n)
end

ext: jl


text:
function pk_qualify(pk::Vector{Float64})
    pk_good = sum(pk .<= 0.5)
    pk_ok = length(pk[pk .<= 0.7]) - pk_good
    pk_bad = length(pk[pk .<= 1]) - pk_good - pk_ok
    (good=pk_good, ok=pk_ok, bad=pk_bad, very_bad=sum(pk .> 1))
end

@userplot Pk_Plot

@recipe function f(p::Pk_Plot)
    if length(p.args) != 1 || !(p.args[1] isa AbstractVector{<:Real})
        error("PSIS diagnostic plots should be given one vector. Got: ", typeof(p.args))
    end
    pk = p.args[1]

    # default plotting attributes
    title --> "PSIS diagnostic plot"
    legend --> :topleft
    xguide --> "Datapoint"
    yguide --> "Pareto shape k"

    @series begin
        seriestype := :scatter
        marker --> 2.5
        label := "Pk points"
        pk
    end

    @series begin
        seriestype := :hline
        label := "pk = 0.5"
        [0.5]
    end

    @series begin
        seriestype := :hline
        label := "pk = 0.7"
        [0.7]
    end

    @series begin
        seriestype := :hline
        label := "pk = 1.0"
        [1.0]
    end
end


text:
using Statistics
using ProgressMeter
using MATLAB
using MAT
using Dates
include("atn_coupled_module.jl")

N = 10_000          # number of nodes
m = 5               # number of active links per active node
tmax = 1000         # maximum duration of the simulation
Ni = Int(0.01 * N)  # initial number of infects
gamma = -2.1        # exponent of the power law

ActMin = 0.001      # minimum activation probability
ActMax = 1
agent_act_1 =
    (
        (ActMax .^ (gamma + 1) - ActMin .^ (gamma + 1)) .* rand(N) .+
        ActMin .^ (gamma + 1)
    ) .^ (1.0 / (gamma + 1))

agent_act_2 =
    (
        (ActMax .^ (gamma + 1) - ActMin .^ (gamma + 1)) .* rand(N) .+
        ActMin .^ (gamma + 1)
    ) .^ (1.0 / (gamma + 1))

mu_1 = 0.1
mu_2 = 0.01

lambda_1 = 0.2
lambda_2 = 0.3

beta_1 = 0.9
beta_2 = beta_1

Iend_1, Iend_2 = atn_coupled_module.infected_totals(
    tmax,
    agent_act_1,
    agent_act_2,
    m,
    mu_1,
    mu_2,
    lambda_1,
    lambda_2,
    beta_1,
    beta_2,
)

location_mat = "/mnt/storage/MAT/"
right_now = replace(
    replace(replace(string(Dates.now()), "." => "_"), ":" => "_"),
    "-" => "_",
)

file_name = "practice" * right_now
mat_name = location_mat * file_name * ".mat"

# matwrite(
#     mat_name,
#     Dict("AvgInfectedEtaLambda" => AvgInfectedEtaLambda);
#     compress = true,
# )
time = Float64.(1:tmax)
mat"figure(); hold on; plot($time,$Iend_1,'.-');plot($time,$Iend_2,'.-')"


text:
function scomponents_test()
    A = load_matrix_network("cores_example")
    cc = scomponents(MatrixNetwork(A))

    sizes = [15;1;5]
    components = 3
    if cc.sizes != sizes
        error("scomponents failed")
    end
    if length(cc.map) != 21
        error("scomponents failed")
    end
    if cc.number != components
        error("scomponents failed")
    end
    return true
end


text:
module SurNoR_2020

using CSVFiles
using CSV
using DataFrames
using DataFramesMeta
using Statistics
using SpecialFunctions
using Random
using PyPlot
using HypothesisTests

include("Structs_form.jl")
include("Functions_for_model.jl")
include("Functions_for_data.jl")
include("Functions_for_behavioral_landmarks.jl")
include("Functions_for_simulation.jl")
include("Functions_for_EEG.jl")
include("Functions_for_model_evaluation.jl")

end # module


text:
using ImagePartitioning
using Test, TestImages

@testset "ImagePartitioning.jl" begin
    include("partition.jl")
end


text:
immutable EnumerableSelect{T, S, Q} <: Enumerable{T}
    source::S
    f::Q
end

function select{TS}(source::Enumerable{TS}, f::Function, f_expr::Expr)
    T = Base.return_types(f, (TS,))[1]
    S = typeof(source)
    return EnumerableSelect{T,S,FunctionWrapper{T,Tuple{TS}}}(source, f)
end

function start{T,S,Q}(iter::EnumerableSelect{T,S,Q})
    s = start(iter.source)
    return s
end

function next{T,S,Q}(iter::EnumerableSelect{T,S,Q}, s)
    x = next(iter.source, s)
    v = x[1]
    s_new = x[2]
    v_new = iter.f(v)::T
    return v_new, s_new
end

function done{T,S,Q}(iter::EnumerableSelect{T,S,Q}, state)
    return done(iter.source, state)
end


text:
using OrderedCollections

include("graph.jl")
include("spatial.jl")
include("temporal.jl")

INT64_MIN = -(2^63)
INT64_MAX = (2^63) - 1

mutable struct PackStreamHydrator
    graph
    hydration_functions
    use_julia_dates

    function PackStreamHydrator(protocol_version, use_julia_dates=true)
        graph = Graph()
        hydration_functions = Dict()
        merge!(hydration_functions, graph_hydration_functions(graph))
        if protocol_version >= 2
            merge!(hydration_functions, spatial_hydration_functions())
            merge!(hydration_functions, temporal_hydration_functions())
        end
        return new(graph, hydration_functions, use_julia_dates)
    end
end

function hydrate(psh::PackStreamHydrator, values)
    function hydrate_(value)
        if isa(value, Structure)
            if haskey(psh.hydration_functions, value.tag)
                return psh.hydration_functions[value.tag](psh, map(hydrate_, value.fields))
            else
                return value
            end
        elseif isa(value, Vector)
            return map(hydrate_, value)
        elseif isa(value, Dict)
            r = Dict()
            for (k, v) in value
                r[k] = hydrate_(v)
            end
            return r
        else
            return value
        end
    end

    return map(hydrate_, values)
end

function hydrate_records(psh::PackStreamHydrator, keys, record_values)
    return [OrderedDict{String, Any}(zip(keys, hydrate(psh, value))) for value in record_values]
end

mutable struct PackStreamDehydrator
    supports_bytes
    protocol_version

    function PackStreamDehydrator(protocol_version; supports_bytes=false)
        return new(supports_bytes, protocol_version)
    end
end

function dehydrate(psd::PackStreamDehydrator, values::Vector)

    function dehydrate__(value::Vector)
        return map(dehydrate_, value)
    end

    function dehydrate__(value::Dict)
        r = Dict()
        for (k, v) in value
            r[k] = dehydrate_(v)
        end
        return r
    end

    function dehydrate_(value)
        if (isa(value, Nothing) || isa(value, Bool) || isa(value, AbstractFloat) || isa(value, String))
            return value
        elseif isa(value, Integer)
            if INT64_MIN <= value <= INT64_MAX
                return value
            else
                throw(ErrorException("Integer out of range"))
            end
        elseif isa(value, Vector{UInt8})
            if psd.supports_bytes
                return value
            else
                throw(TypeError())
            end
        elseif isa(value, Vector) || isa(value, Dict)
            dehydrate__(value)
        elseif psd.protocol_version >= 2
            dehydrate(value)
        else
            throw(TypeError())
        end
    end

    return map(dehydrate_, values)
end


text:
function getCompleteVariableList(
    mb      :: ModelBlock,
    vartype :: Symbol,
)

    d = nothing

    if vartype == :dynamic
        d = OrderedDict(

            "TEMP"               => ( mb.fi.sv[:TEMP], :T, :mask ),
            "SALT"               => ( mb.fi.sv[:SALT], :T, :mask ),
            "UVEL"               => ( mb.fi.sv[:UVEL], :U, nothing ),
            "VVEL"               => ( mb.fi.sv[:VVEL], :V, nothing ),
            "WVEL"               => ( mb.fi.sv[:WVEL], :W, nothing ),
            "TAUX"               => ( mb.fi.TAUX,        :sT, :mask ),
            "TAUY"               => ( mb.fi.TAUY,        :sT, :mask ),
            "TAUX_east"          => ( mb.fi.TAUX_east,   :sT, :mask ),
            "TAUY_north"         => ( mb.fi.TAUY_north,  :sT, :mask ),
            "SWFLX"              => ( mb.fi.SWFLX,       :sT, :mask ),
            "NSWFLX"             => ( mb.fi.NSWFLX,      :sT, :mask ),
            "VSFLX"              => ( mb.fi.VSFLX,       :sT, :mask ),
            "ADVT"               => ( mb.fi.sv[:ADVT],     :T, :mask ),
            "ADVS"               => ( mb.fi.sv[:ADVS],     :T, :mask ),
            "HMXL"               => ( mb.fi.HMXL,          :sT, :mask ),

            "WKRSTT"       => ( mb.fi.sv[:WKRSTT], :T, :mask),
            "WKRSTS"       => ( mb.fi.sv[:WKRSTS], :T, :mask),

            "VDIFFT"       => ( mb.fi.sv[:VDIFFT], :T, :mask),
            "VDIFFS"       => ( mb.fi.sv[:VDIFFS], :T, :mask),

            "Q_FRZMLTPOT"        => ( mb.fi.Q_FRZMLTPOT,  :sT, :mask ),
            "Q_FRZMLTPOT_NEG"    => ( mb.fi.Q_FRZMLTPOT_NEG, :sT, :mask ),
            "Q_FRZHEAT"          => ( mb.fi.Q_FRZHEAT,    :sT, :mask ),
            "Q_LOST"             => ( mb.fi.Q_LOST,       :sT, :mask ),
            "CHKTEMP"            => ( mb.tmpfi.sv[:CHKTEMP],  :sT, :mask ),
            "CHKSALT"            => ( mb.tmpfi.sv[:CHKSALT],  :sT, :mask ),

            "INTMTEMP"           => ( mb.tmpfi.sv[:INTMTEMP], :T, :mask),
            "INTMSALT"           => ( mb.tmpfi.sv[:INTMSALT], :T, :mask),

            "WKRST_TARGET_TEMP"         => nothing,
            "WKRST_TARGET_SALT"         => nothing,
            "QFLXT"          => nothing,
            "QFLXS"          => nothing,

            "Ks_H_U"    => ( mb.tmpfi.check_usage[:Ks_H_U], :U, nothing),
            "Ks_H_V"    => ( mb.tmpfi.check_usage[:Ks_H_V], :V, nothing),

        )

        if mb.tmpfi.datastream != nothing

            if mb.ev.config["weak_restoring"] == "on"
                d["WKRST_TARGET_TEMP"] = ( mb.tmpfi.datastream["TEMP"], :T , :mask)
                d["WKRST_TARGET_SALT"] = ( mb.tmpfi.datastream["SALT"], :T , :mask)
            end

        end

        d["QFLXT"] = ( mb.fi.sv[:QFLXT], :T, :mask )
        d["QFLXS"] = ( mb.fi.sv[:QFLXS], :T, :mask )

    elseif vartype == :static

        d = Dict(
            # COORDINATE
            "deepmask_T"=> ( mb.ev.topo.deepmask_T, :T, nothing),
            "topoz_sT"  => ( mb.ev.topo.topoz_sT, :sT, nothing),
            "Nz_bot_sT" => ( mb.ev.topo.Nz_bot_sT, :sT, nothing),
            "area_sT"   => ( mb.ev.gd_slab.Δx_T .* mb.ev.gd_slab.Δy_T, :sT, nothing),
            "mask_sT"   => ( mb.ev.topo.sfcmask_sT, :sT, nothing),
            "z_cW"      => ( reshape(mb.ev.config["z_w"], :, 1, 1), :cW, nothing),
            "dz_cT"     => ( mb.ev.gd.Δz_T[:, 1:1, 1:1], :cT, nothing),
            "lon_sT"    => ( rad2deg.(mb.ev.gd_slab.λ_T), :sT, nothing),
            "lat_sT"    => ( rad2deg.(mb.ev.gd_slab.ϕ_T), :sT, nothing),
        )

    else
        throw(ErrorException("Unknown vartype: " * string(vartype)))
    end

    return d
end

function getVarDesc(varname)
    return haskey(HOOM.var_desc, varname) ? HOOM.var_desc[varname] : Dict()
end

function getDynamicVariableList(
    mb :: ModelBlock;
    varnames :: AbstractArray{String} = Array{String}(undef,0),
    varsets  :: AbstractArray{String} = Array{String}(undef,0),
)

    all_varlist = getCompleteVariableList(mb, :dynamic)

    output_varnames = []

    for varname in varnames
        if haskey(all_varlist, varname)
            push!(output_varnames, varname)
        else
            throw(ErrorException(format("Varname {:s} does not exist.", varname)))
        end
    end

    for varset in varsets

        if varset == "ALL"

            append!(output_varnames, keys(all_varlist))

        elseif varset == "ESSENTIAL"

            append!(output_varnames, [

                "TEMP",
                "SALT",
                "UVEL",
                "VVEL",
                "WVEL",
                "TAUX",
                "TAUY",
                "TAUX_east",
                "TAUY_north",
                "SWFLX",
                "NSWFLX",
                "VSFLX",
                "ADVT",
                "ADVS",
                "HMXL",

                "WKRSTT",
                "WKRSTS",

                "VDIFFT",
                "VDIFFS",

                "Q_FRZMLTPOT",
                "Q_LOST",
                "CHKTEMP",
                "CHKSALT",

            ])

        else

            throw(ErrorException("Unknown varset: " * string(varset)))

        end
    end

    function makeSubset(dict, keys)
        subset_dict = Dict()
        for k in keys
            if dict[k] != nothing
                subset_dict[k] = dict[k]
            else
                println("$(k) is nothing")
            end
        end
        return subset_dict
    end

    #println("NOW AGAIN: ", output_varnames)

    output_varlist = makeSubset(all_varlist, output_varnames)

    return output_varlist
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Nexss PROGRAMMER 2.0.0 - Julia
# Default template for JSON Data
# STDIN
NexssStdin = readline();

# JSON in Julia installation
import Pkg;
if !in("JSON", keys(Pkg.installed()))
    # write(stderr,"pkgJSON is not intalled");
    @warn("pkgJSON is not intalled")
    Pkg.add("JSON");
end
import JSON;
# END JSON part

# Modify Data
parsedJson = JSON.parse(NexssStdin);
parsedJson["juliaOutput"] = "Hello from Julia! $VERSION";
# parsedJson["test"] = "test";
NexssStdout = JSON.json(parsedJson)

# STDOUT
write(stdout,NexssStdout);


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

# I want to remove this dependence
using Primes

import Base: ==, ∩, +
import Primes: factor

include("cyclotomic.jl")
include("linearcode.jl")

abstract type AbstractCyclicCode <: AbstractLinearCode end
abstract type AbstractBCHCode <: AbstractCyclicCode end
abstract type AbstractReedSolomonCode <: AbstractBCHCode end

# need to write == functions
struct CyclicCode <: AbstractCyclicCode
    F::FqNmodFiniteField # base field
    E::FqNmodFiniteField # splitting field
    R::FqNmodPolyRing # polynomial ring of generator polynomial
    β::fq_nmod # n-th root of primitive element of splitting field
    n::Integer
    k::Integer
    d::Union{Integer, Missing}
    b::Integer
    δ::Integer
    qcosets::Vector{Vector{Int64}}
    qcosetsreps::Vector{Int64}
    defset::Vector{Int64}
    g::fq_nmod_poly
    h::fq_nmod_poly
    e::fq_nmod_poly
    G::fq_nmod_mat
    Gorig::Union{fq_nmod_mat, Missing}
    H::fq_nmod_mat
    Horig::Union{fq_nmod_mat, Missing}
    Gstand::fq_nmod_mat
    Hstand::fq_nmod_mat
    weightenum::Union{WeightEnumerator, Missing}
end

struct BCHCode <: AbstractBCHCode
    F::FqNmodFiniteField # base field
    E::FqNmodFiniteField # splitting field
    R::FqNmodPolyRing # polynomial ring of generator polynomial
    β::fq_nmod # n-th root of primitive element of splitting field
    n::Integer
    k::Integer
    d::Union{Integer, Missing}
    b::Integer
    δ::Integer
    qcosets::Vector{Vector{Int64}}
    qcosetsreps::Vector{Int64}
    defset::Vector{Int64}
    g::fq_nmod_poly
    h::fq_nmod_poly
    e::fq_nmod_poly
    G::fq_nmod_mat
    Gorig::Union{fq_nmod_mat, Missing}
    H::fq_nmod_mat
    Horig::Union{fq_nmod_mat, Missing}
    Gstand::fq_nmod_mat
    Hstand::fq_nmod_mat
    weightenum::Union{WeightEnumerator, Missing}
end

struct ReedSolomonCode <: AbstractReedSolomonCode
    F::FqNmodFiniteField # base field
    E::FqNmodFiniteField # splitting field
    R::FqNmodPolyRing # polynomial ring of generator polynomial
    β::fq_nmod # n-th root of primitive element of splitting field
    n::Integer
    k::Integer
    d::Union{Integer, Missing}
    b::Integer
    δ::Integer
    qcosets::Vector{Vector{Int64}}
    qcosetsreps::Vector{Int64}
    defset::Vector{Int64}
    g::fq_nmod_poly
    h::fq_nmod_poly
    e::fq_nmod_poly
    G::fq_nmod_mat
    Gorig::Union{fq_nmod_mat, Missing}
    H::fq_nmod_mat
    Horig::Union{fq_nmod_mat, Missing}
    Gstand::fq_nmod_mat
    Hstand::fq_nmod_mat
    weightenum::Union{WeightEnumerator, Missing}
end

function _generatorpolynomial(R::FqNmodPolyRing, β::fq_nmod, Z::Vector{Int64})
    # from_roots(R, [β^i for i in Z]) - R has wrong type for this
    g = one(R)
    for i in Z
        g *= (gen(R) - β^i)
    end

    return g
end
_generatorpolynomial(R::FqNmodPolyRing, β::fq_nmod, qcosets::Vector{Vector{Int64}}) = _generatorpolynomial(R, β, vcat(qcosets...))

function _generatormatrix(F::FqNmodFiniteField, n::Integer, k::Integer, g::fq_nmod_poly)
    # if g = x^10 + α^2*x^9 + x^8 + α*x^7 + x^3 + α^2*x^2 + x + α
    # g.coeffs = [α  1  α^2  1  0  0  0  α  1  α^2  1]
    coeffs = collect(coefficients(g))
    len = length(coeffs)
    k + len - 1 <= n || error("Too many coefficients for $k shifts in _generatormatrix.")

    G = zero_matrix(F, k, n)
    for i in 1:k
        G[i, i:i + len - 1] = coeffs
    end
    return G
end

function definingset(nums::Vector{Int64}, q::Integer, n::Integer, flat::Bool=true)
    arr = Vector{Vector{Int64}}()
    arrflat = Vector{Int64}()
    for x in nums
        Cx = cyclotomiccoset(x, q, n)
        if Cx[1] ∉ arrflat
            arrflat = [arrflat; Cx]
            push!(arr, Cx)
        end
    end

    !flat || return sort!(vcat(arr...))
    return arr
end

function _idempotent(g::fq_nmod_poly, h::fq_nmod_poly, n::Integer)
    # solve 1 = a(x) g(x) + b(x) h(x) for a(x) then e(x) = a(x) g(x) mod x^n - 1
    d, a, b = gcdx(g, h)
    return mod(g * a, gen(parent(g))^n - 1)
end

# MattsonSolomontransform(f, n)
# inverseMattsonSolomontransform

basefield(C::AbstractCyclicCode) = C.F
splittingfield(C::AbstractCyclicCode) = C.E
polynomialring(C::AbstractCyclicCode) = C.R
primitiveroot(C::AbstractCyclicCode) = C.β
offset(C::AbstractCyclicCode) = C.b
designdistance(C::AbstractCyclicCode) = C.δ
qcosets(C::AbstractCyclicCode) = C.qcosets
qcosetsreps(C::AbstractCyclicCode) = C.qcosetsreps
definingset(C::AbstractCyclicCode) = C.defset
generatorpolynomial(C::AbstractCyclicCode) = C.g
paritycheckpolynomial(C::AbstractCyclicCode) = C.h
idempotent(C::AbstractCyclicCode) = C.e
isprimitive(C::AbstractCyclicCode) = length(n) == order(basefield(C)) - 1
isnarrowsense(C::AbstractCyclicCode) = iszero(offset(C)) # should we define this as b = 1 instead?
isreversible(C::AbstractCyclicCode) = [length(C) - i for i in defset] ⊆ defset

function Base.show(io::IO, C::AbstractCyclicCode)
    if get(io, :compact, false)
        # to use type "show(IOContext(stdout, :compact=>true), C)" instead
        if typeof(C) <: ReedSolomonCode
            println(io, "[$(length(C)), $(dimension(C)), ≥$(designdistance(C)); $(offset(C))]_$(order(field(C))) Reed Solomon code.")
        elseif typeof(C) <: BCHCode
            println(io, "[$(length(C)), $(dimension(C)), ≥$(designdistance(C)); $(offset(C))]_$(order(field(C))) BCH code over splitting field GF($(order(splittingfield(C)))).")
        else
            println(io, "[$(length(C)), $(dimension(C)), ≥$(designdistance(C)); $(offset(C))]_$(order(field(C))) cyclic code over splitting field GF($(order(splittingfield(C)))).")
        end
    else
        if typeof(C) <: ReedSolomonCode
            println(io, "[$(length(C)), $(dimension(C)), ≥$(designdistance(C)); $(offset(C))]_$(order(field(C))) Reed Solomon code.")
        elseif typeof(C) <: BCHCode
            println(io, "[$(length(C)), $(dimension(C)), ≥$(designdistance(C)); $(offset(C))]_$(order(field(C))) BCH code over splitting field GF($(order(splittingfield(C)))).")
        else
            println(io, "[$(length(C)), $(dimension(C)), ≥$(designdistance(C)); $(offset(C))]_$(order(field(C))) cyclic code over splitting field GF($(order(splittingfield(C)))).")
        end
        println(io, "$(order(field(C)))-Cyclotomic cosets: ")
        for (i, x) in enumerate(qcosetsreps(C))
            if i == 1
                print(io, "\tC_$x ∪ ")
            elseif i == 1 && i == length(qcosetsreps(C))
                println(io, "\tC_$x")
            elseif i != length(qcosetsreps(C))
                print(io, "C_$x ∪ ")
            else
                println(io, "C_$x")
            end
        end
        println(io, "Generator polynomial:")
        println(io, "\t", generatorpolynomial(C))
        println(io, "Generator matrix: $(dimension(C)) × $(length(C))")
        for i in 1:dimension(C)
            print(io, "\t")
            for j in 1:length(C)
                if j != length(C)
                    print(io, "$(C.G[i, j]) ")
                elseif j == length(C) && i != dimension(C)
                    println(io, "$(C.G[i, j])")
                else
                    print(io, "$(C.G[i, j])")
                end
            end
        end
    end
end

function finddelta(n::Integer, cosets::Vector{Vector{Int64}})
    # println(cosets)
    defset = sort!(vcat(cosets...))
    runs = Vector{Vector{Int64}}()
    for x in defset
        useddefset = Vector{Int64}()
        reps = Vector{Int64}()
        cosetnum = 0
        for i in 1:length(cosets)
            if x ∈ cosets[i]
                cosetnum = i
                append!(useddefset, cosets[i])
                append!(reps, x)
                break
            end
        end

        y = x + 1
        while y ∈ defset
            if y ∈ useddefset
                append!(reps, y)
            else
                cosetnum = 0
                for i in 1:length(cosets)
                    if y ∈ cosets[i]
                        cosetnum = i
                        append!(useddefset, cosets[i])
                        append!(reps, y)
                        break
                    end
                end
            end
            y += 1
        end
        push!(runs, reps)
    end

    runlens = [length(i) for i in runs]
    (consec, ind) = findmax(runlens)
    # there are δ - 1 consecutive numbers for designed distance δ
    δ = consec + 1
    # start of run
    offset = runs[ind][1]
    # BCH Bound is thus d≥δ

    # moving to Hartmann-Tzeng Bound refinement
    currbound = δ
    if consec > 1
        for A in runs
            if length(A) == consec
                for b in 1:(n - 1)
                    if gcd(b, n) ≤ δ
                        for s in 0:(δ - 2)
                            B = [mod(j * b, n) for j in 0:s]
                            AB = [x + y for x in A for y in B]
                            if AB ⊆ defset
                                if currbound < δ + s
                                    currbound = δ + s
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    return δ, offset, currbound
end

function largestconsecrun(arr::Vector{Int64})
    n = length(arr)
    maxlen = 1
    for i = 1:n
        mn = arr[i]
        mx = arr[i]

        for j = (i + 1):n
            mn = min(mn, arr[j])
            mx = max(mx, arr[j])

            if (mx - mn) == (j - i)
                maxlen = max(maxlen, mx - mn + 1)
            end
        end
    end

    return maxlen
end

function dualdefiningset(arr::Vector{Int64}, n::Integer)
    full = [i for i in 0:(n - 1)]
    temp = Vector{Int64}()
    for i in full
        if i ∉ arr
            append!(temp, i)
        end
    end

    return sort!([mod(n - i, n) for i in temp])
end

function CyclicCode(q::Integer, n::Integer, cosets::Vector{Vector{Int64}}, verify::Bool=true)
    !(q <= 1 || n <= 1) ||error("Invalid parameters past to CyclicCode constructor: q = $q, n = $n.")

    if !isprime(q)
        # this used to work with just AbstractAlgebra
        factors = Primes.factor(q)
        if length(factors) != 1
            error("There is no finite field of order $(prod(factors)).")
        end
        (p, t), = factors
    else
        p = q
        t = 1
    end

    F, _ = FiniteField(p, t, "α") # changed to keep ReedSolomonCodes printing α's'
    deg = ord(n, q)
    E, α = FiniteField(p, t * deg, "α")
    R, _ = PolynomialRing(E, "x")
    β = α^(div(q^deg - 1, n))
    # println("here so far")

    defset = sort!(vcat(cosets...))
    k = n - length(defset)
    # println(k)
    comcosets = complementqcosets(q, n, cosets)
    # println("here 2")
    g = _generatorpolynomial(R, β, defset)
    h = _generatorpolynomial(R, β, vcat(comcosets...))
    e = _idempotent(g, h, n)
    G = _generatormatrix(F, n, k, g)
    H = _generatormatrix(F, n, n - k, reverse(h))
    # println("here 3")
    # println(G)
    Gstand, Hstand = _standardform(G)
    # println("here 4")
    δ, b, HT = finddelta(n, cosets)
    # println("here 5")

    if verify
        # println("above")
        flag, htest = divides(gen(R)^n - 1, g)
        # println(flag)
        # println(htest)
        flag || error("Incorrect generator polynomial, does not divide x^$n - 1.")
        # println(htest, parent(h))
        # println(h, parent(h))
        htest == h || error("Division of x^$n - 1 by the generator polynomial does not yield the constructed parity check polynomial.")
        # if htest != h
        #     println("sucks")
        #     # error("test")
        # end
        if size(H) == (n - k, k)
            H = deepcopy(H')
        end
        !(!iszero(G * H') || !iszero(H * G')) || error("Generator and parity check matrices are not transpose orthogonal.")
        for r in 1:size(Gstand, 1)
            iszero(Gstand[r, :] * H') || error("Column swap appeared in _standardform.")
        end

        # check e=e^2
    end

    if δ >= 2 && defset == definingset([i for i = b:(b + δ - 2)], q, n, true)
        if deg == 1 && n == q - 1
            return ReedSolomonCode(F, E, R, β, n, k, n - k + 1, b, HT, cosets,
                sort!([arr[1] for arr in cosets]), defset, g, h, e, G, missing,
                H, missing, Gstand, Hstand, missing)
        end

        return BCHCode(F, E, R, β, n, k, missing, b, HT, cosets,
            sort!([arr[1] for arr in cosets]), defset, g, h, e, G, missing, H,
            missing, Gstand, Hstand, missing)
    end

    return CyclicCode(F, E, R, β, n, k, missing, b, HT, cosets,
        sort!([arr[1] for arr in cosets]), defset, g, h, e, G, missing, H, missing,
        Gstand, Hstand, missing)
end

# currently untested - not fully fixed yet
function CyclicCode(q::Integer, n::Integer, g::fq_nmod_poly, verify::Bool=true)
    flag, htest = divides(gen(R)^n - 1, g)
    flag || error("Given polynomial does not divide x^$n - 1.")

    if !isprime(q)
        factors = factor(q)
        if length(factors) != 1
            error("There is no finite field of order $(prod(factors)).")
        end
        (q, n), = factors
    end

    R = parent(g)
    F = base_ring(R)
    α = gen(F)
    t = ord(n, q)
    β = α^(div(q^t - 1, n))

    dic = Dict{fq_nmod, Int64}()
    for i in 1:n
        dic[β^i] = i
    end
    rts = roots(g)
    defset = sort!([dic[rt] for rt in rts])
    qcosets = definingset(defset, q, n, false)
    k = n - length(defset)
    comcosets = complementqcosets(q, n, qcosets)
    e = _idempotent(g, h)
    G = _generatormatrix(n, k, g)
    H = _generatormatrix(n, n - k, reverse(htest))
    Gstand, Hstand = _standardform(G)
    _, _, HT = finddelta(n, qcosets)

    if verify
        h, _, _, _ = _generatorpolynomial(q, n, vcat(comcosets...))
        htest == h || error("Division of x^$n - 1 by the generator polynomial does not yield the constructed parity check polynomial.")
        if size(H) == (n - k, k)
            H = deepcopy(H')
        end
        !(!iszero(G * H') || !iszero(H * G')) || error("Generator and parity check matrices are not transpose orthogonal.")
        for r in 1:size(Gstand, 1)
            iszero(Gstand[r, :] * H') || error("Column swap appeared in _standardform.")
        end
    end

    return CyclicCode(F, E, R, β, n, k, missing, b, HT, qcosets,
        sort!([arr[1] for arr in qcosets]), defset, g, h, e, G, H, Gstand, Hstand,
        missing)
end

# self orthogonal cyclic codes are even-like
# does this require them too have even minimum distance?
# self orthogonal code must contain all of its self orthogonal q-cosets and at least one of every q-coset pair
function BCHCode(q::Integer, n::Integer, δ::Integer, b::Integer=0, verify::Bool=true)
    δ >= 2 || error("BCH codes require δ ≥ 2 but the constructor was given δ = $δ.")
    !(q <= 1 || n <= 1) || error("Invalid parameters past to BCHCode constructor: q = $q, n = $n.")

    if !isprime(q)
        factors = Primes.factor(q)
        if length(factors) != 1
            error("There is no finite field of order $(prod(factors)).")
        end
        (p, t), = factors
    else
        p = q
        t = 1
    end

    F, _ = FiniteField(p, t, "α") # changed to keep ReedSolomonCodes printing α's'
    deg = ord(n, q)
    E, α = FiniteField(p, t * deg, "α")
    R, _ = PolynomialRing(E, "x")
    β = α^(div(q^deg - 1, n))

    cosets = definingset([i for i = b:(b + δ - 2)], q, n, false)
    defset = sort!(vcat(cosets...))
    k = n - length(defset)
    comcosets = complementqcosets(q, n, cosets)
    g = _generatorpolynomial(R, β, defset)
    h = _generatorpolynomial(R, β, vcat(comcosets...))
    e = _idempotent(g, h, n)
    G = _generatormatrix(F, n, k, g)
    H = _generatormatrix(F, n, n - k, reverse(h))
    Gstand, Hstand = _standardform(G)
    δ, b, HT = finddelta(n, cosets)

    if verify
        flag, htest = divides(gen(R)^n - 1, g)
        flag || error("Incorrect generator polynomial, does not divide x^$n - 1.")
        htest == h || error("Division of x^$n - 1 by the generator polynomial does not yield the constructed parity check polynomial.")
        if size(H) == (n - k, k)
            H = deepcopy(H')
        end
        !(!iszero(G * H') || !iszero(H * G')) || error("Generator and parity check matrices are not transpose orthogonal.")
        for r in 1:size(Gstand, 1)
            iszero(Gstand[r, :] * H') || error("Column swap appeared in _standardform.")
        end
    end

    if deg == 1 && n == q - 1
        return ReedSolomonCode(F, E, R, β, n, k, n - k + 1, b, HT, cosets,
            sort!([arr[1] for arr in cosets]), defset, g, h, e, G, missing, H,
            missing, Gstand, Hstand, missing)
    end

    return BCHCode(F, E, R, β, n, k, n - k + 1, b, HT, cosets,
        sort!([arr[1] for arr in cosets]), defset, g, h, e, G, missing, H,
        missing, Gstand, Hstand, missing)
end

function ReedSolomonCode(q::Integer, δ::Integer, b::Integer=0, verify::Bool=true)
    δ >= 2 || error("Reed Solomon codes require δ ≥ 2 but the constructor was given δ = $δ.")
    q > 4 || error("Invalid or too small parameters past to ReedSolomonCode constructor: q = $q.")

    # n = q - 1
    # if ord(n, q) != 1
    #     error("Reed Solomon codes require n = q - 1.")
    # end

    if !isprime(q)
        factors = factor(q)
        if length(factors) != 1
            error("There is no finite field of order $(prod(factors)).")
        end
        (p, t), = factors
    else
        p = q
        t = 1
    end

    F, α = FiniteField(p, t, "α") # changed to keep ReedSolomonCodes printing α's'
    R, _ = PolynomialRing(F, "x")

    n = q - 1
    cosets = definingset([i for i = b:(b + δ - 2)], q, n, false)
    defset = sort!(vcat(cosets...))
    k = n - length(defset)
    comcosets = complementqcosets(q, n, cosets)
    g = _generatorpolynomial(R, α, defset)
    h = _generatorpolynomial(R, α, vcat(comcosets...))
    e = _idempotent(g, h, n)
    G = _generatormatrix(F, n, k, g)
    H = _generatormatrix(F, n, n - k, reverse(h))
    Gstand, Hstand = _standardform(G)
    δ, b, HT = finddelta(n, cosets)

    if verify
        flag, htest = divides(gen(R)^n - 1, g)
        flag || error("Incorrect generator polynomial, does not divide x^$n - 1.")
        htest == h || error("Division of x^$n - 1 by the generator polynomial does not yield the constructed parity check polynomial.")
        if size(H) == (n - k, k)
            H = deepcopy(H')
        end
        !(!iszero(G * H') || !iszero(H * G')) || error("Generator and parity check matrices are not transpose orthogonal.")
        for r in 1:size(Gstand, 1)
            iszero(Gstand[r, :] * H) || error("Column swap appeared in _standardform.")
        end
    end

    return ReedSolomonCode(F, F, R, α, n, k, n - k + 1, b, HT, cosets,
        sort!([arr[1] for arr in cosets]), defset, g, h, e, G, missing, H,
        missing, Gstand, Hstand, missing)
end

function complement(C::AbstractCyclicCode, verify::Bool=true)
    D = CyclicCode(Int64(order(field(C))), length(C),
        complementqcosets(Int64(order(field(C))), length(C), qcosets(C)))
    if verify
        if paritycheckpolynomial(C) != generatorpolynomial(D) || idempotent(D) != (1 - idempotent(C))
            error("Error constructing the complement cyclic code.")
        end
    end

    return D
end

# C1 ⊆ C2 iff g_2(x) | g_1(x) iff T_2 ⊆ T_1
⊆(C1::AbstractCyclicCode, C2::AbstractCyclicCode) = C2.defset ⊆ C1.defset
issubcode(C1::AbstractCyclicCode, C2::AbstractCyclicCode) = C1 ⊆ C2

function ==(C1::AbstractCyclicCode, C2::AbstractCyclicCode)
    # should also check primitive root but so far the user is not given a choice here
    return field(C1) == field(C2) && length(C1) == length(C2) && definingset(C1) == definingset(C2)
end

function dual(C::AbstractCyclicCode)
    # one is even-like and the other is odd-like
    return CyclicCode(Int64(order(field(C))), length(C),
        dualqcosets(Int64(order(field(C))), length(C), qcosets(C)))
end

# this checks def set, need to rewrite == for linear first
isselfdual(C::AbstractCyclicCode) = C == dual(C)

# don't think this is necessary in order to invoke the ⊆ for CyclicCode
# function isselforthogonal(C::AbstractCyclicCode)
#     # A code is self-orthogonal if it is a subcode of its dual.
#     return C ⊆ dual(C)
# end

# function μa(C::CyclicCode)
#     # check gcd(a, n) = 1
#     # technically changes g(x) and e(x) but the q-cosets are the same?
# end

function ∩(C1::AbstractCyclicCode, C2::AbstractCyclicCode)
    # has generator polynomial lcm(g_1(x), g_2(x))
    # has generator idempotent e_1(x) e_2(x)

    if field(C1) == field(C2) && length(C1) == length(C2)
        return CyclicCode(Int64(order(field(C1))), length(C1),
            definingset(definingset(C1) ∪ definingset(C2), Int64(order(field(C1))), length(C1), false))
    else
        error("Cannot intersect two codes over different base fields or lengths.")
    end
end

function +(C1::AbstractCyclicCode, C2::AbstractCyclicCode)
    # has generator polynomial gcd(g_1(x), g_2(x))
    # has generator idempotent e_1(x) + e_2(x) - e_1(x) e_2(x)
    if field(C1) == field(C2) && length(C1) == length(C2)
        defset = definingset(C1) ∩ definingset(C2)
        if length(defset) != 0
            return CyclicCode(Int64(order(field(C1))), length(C1),
                definingset(defset, Int64(order(field(C1))), length(C1), false))
        else
            error("Addition of codes has empty defining set.")
        end
    else
        error("Cannot add two codes over different base fields or lengths.")
    end
end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
using Primes, Printf

println("Primes ≤ $hi written in common bases.")
@printf("%8s%8s%8s%8s", "bin", "oct", "dec", "hex")
for i in primes(50)
    @printf("%8s%8s%8s%8s\n", bin(i), oct(i), dec(i), hex(i))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/julia
# -*- Julia 0.7 -*-
#
# Calculate Julia Set with Julia.
#
# The resulting image will be saved in PNM format to file result.ppm
# This ppm file can be viewed directly, or you can convert it into
# another format like this:
#   $ ffmpeg -i result.pnm -vcodec png output.png
# This command provides performance information.
#   $ sudo /usr/bin/time perf stat julia main.jl
#
# https://en.wikipedia.org/wiki/Julia_set      for "julia set"
# http://mathworld.wolfram.com/JuliaSet.html   for "julia set"
# https://en.wikipedia.org/wiki/Netpbm_format  for "pnm picture format"
#
# Copyright © Zhou Mo <cdluminate AT gmail.com>
# MIT License

using SharedArrays
using Printf

# configure
#C = -0.62772 -0.42193im
C = 0.285 + 0.0im
#C = randn() + randn() * im
xmin = -1.8
xmax = 1.8
ymin = -1.8
ymax = 1.8
samples = 1920
@printf("%s", "Dump configuration
 -> C       = $(C)
 -> xrange  = [ $(xmin), $(xmax) ]
 -> yrange  = [ $(ymin), $(ymax) ]
 -> samples = $(samples)
")

# helper functions
function getPNMHeader(width, height, magic="P2")
   # Magic "P2" for gray PNM picture
   return "$(magic) $(width) $(height) 255
   "
end

function getCount(z, c, maxiter=255, threshold=2)
   n = 0
   while abs(z)<threshold && n < maxiter
      z = z*z + c
      n += 1
   end
   return n
end

function normalize(mat)
   return (mat - minimum(mat)) / (maximum(mat) - minimum(mat))
end

function fakecolorize(mat)
   I, J = size(mat)
   ret = zeros(UInt8, I, J*3)
   for i = 1:I
      for j = 1:J
         # ret[i, (j-1)*3+1:j*3] = mat[i,j]
         needle = round(UInt8, mat[i,j])
         value = [ 0, 0, 0 ]
         if needle <= 63
            value = [ 0, 4*needle, 255 ]
         elseif needle <= 127
            value = [ 0, 255, 511-4*needle ]
         elseif needle <= 191
            value = [ 4*needle-511, 255, 0 ]
         else
            value = [ 255, 1023- 4*needle, 0 ]
         end
         ret[i, (j-1)*3+1:j*3] = value
      end
   end
   return ret
end

# calculate
output = SharedArray{Int}(samples, samples)
@printf("Calculating Julia Set\n")
@time for (i,x) in enumerate(linspace(xmin, xmax, samples))
   @printf("\0337 -> progress: %.2f%%\0338", i*100/samples)
   for (j,y) in enumerate(linspace(ymin, ymax, samples))
      # FIXME: can't use @parallel?
      # FIXME: why should I rotate the resulting image like this?
      output[j,i] = getCount(x + y * im, C)
   end
end
println("")

# post process
@printf("Post processing\n")
output = round.(UInt8, 255*normalize(output))

# save picture
#@printf("Saving gray result to result.gray.pnm\n")
#f = open("result.gray.pnm", "w+")
#write(f, getPNMHeader(samples, samples))
#writedlm(f, output, " ")
#close(f)

println("Saving color result to result.$(C).color.pnm")
fc = open("result.$(C).color.pnm", "w+")
write(fc, getPNMHeader(samples, samples, "P3"))
writedlm(fc, fakecolorize(output), " ")
close(fc)

@printf("Done\n")

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using NCDatasets,NetCDF,Dates,MeshArrays,Printf

"""
    NCvar

Data structure containing information needed to write a NetCDF file. This
includes a list of filenames (see `Bindata`) if the data is not loaded into memory.
"""
struct NCvar
    name::String
    units::String
    dims
    values
    atts::Union{Dict,Nothing}
    backend::Module
end

function NCvar(name::String, units::String, dims::NCvar, values, atts::Union{Dict,Nothing}, backend::Module)
    return NCvar(name, units, [dims], values, atts, backend)
end

"""
    BinData

Data structure containing a string or an array of strings (NetCDF
    file names) as well as metadata needed to read a file.
"""
struct BinData # Pointer to data stored in binary files- contains info needed to read in
    fnames::Union{Array{String},String}
    precision::Type
    iosize::Tuple
    fldidx::Int
end

"""
    BinData(fnames::Union{Array{String},String},precision::Type,iosize::Tuple)

Construct a BinData struct for files that contain one field.
"""
function BinData(fnames::Union{Array{String},String},precision::Type,iosize::Tuple)
    return BinData(fnames,precision,iosize,1)
end

function BinData(fnames::Union{Array{String},String},precision::Type,iosize::Tuple,meta::Dict,fldname::String)
    return BinData(fnames,precision,iosize,findfirst(meta["fldList"] .== fldname)[2])
end

"""
    NCData

Data structure containing a string or an array of strings (file names) of
    NetCDF files as well as information needed to read a file.
"""
struct NCData
    fname::AbstractString
    varname::AbstractString
    backend::Module
    precision::Type
end

function getnsteps(d)
    if isa(d,BinData)
        if isa(d.fnames,Array)
            nsteps = length(d.fnames)
        else
            nsteps = 1
        end
    elseif isa(d,NCData)
        ds = Dataset(d.fname)
        tdim = dimnames(ds[d.varname])[end]
        timeUnits = ["minutes","seconds","hours","days","minute","second","hour","day"]
        if any(occursin.(timeUnits,Ref(lowercase(ds[tdim].attrib["units"]))))
            nsteps = length(ds[tdim])
        else
            nsteps = 1
        end
    elseif isa(d,TileData)
        nsteps = getnsteps(d.vals)
    elseif isa(d,Array) && isa(d[1],Array)
        nsteps = length(d)
    else
        nsteps = 1
    end

    return nsteps
end

"""
    TileData{T}

Data structure containing either a `MeshArray` struct or `BinData` struct (see `vals`),
    `MeshArray` structs describing the tile layout (`tileinfo`), and other information for
    reading/writing tile data.
"""
struct TileData{T}
    vals::T
    tileinfo::Dict
    tilesize::Tuple
    precision::Type
    numtiles::Int
end

"""
    TileData(vals,tilesize::Tuple)

Construct a TileData struct. First generate the tileinfo, precision, and numtiles attributes.
"""
function TileData(vals,tilesize::Tuple,grid::gcmgrid)
    ni,nj = tilesize
    tileinfo = findtiles(ni,nj,grid)
    if isa(vals,BinData)
        prec = vals.precision
    else
        prec = Float32
    end
    return TileData(vals,tileinfo,tilesize,Float32,Int(maximum(tileinfo["tileNo"])))
end
TileData(vals,tilesize::Tuple,grid::String="LLC90") = TileData(vals,tilesize::Tuple,GCMGridSpec(grid))

"""
    findidx(A,val)

Helper function for getting the indices for tiles. A is a gcmfaces struct, val is a
    numeric value to get the indices of. Returns a Dict of the indices for each face.
    Not currently exported. Maybe move to MeshArrays?
"""
function findidx(A,val)
    tileidx = Dict()
    if isa(A,MeshArrays.gcmfaces)
        nfaces = A.nFaces
    elseif isa(A,MeshArray)
        nfaces = A.grid.nFaces
    end
    for i = 1:nfaces
        idx = findall(x->x==val,A.f[i])
        tileidx[i] = [x for x in idx]
    end
    return tileidx
end

"""
    gettile(fldvals,tileinfo,tilesize,tilenum::Int)

Helper function for retrieving a tile from a gcmfaces struct as a numeric Array. Not
    currently exported.
"""
function gettile(fldvals,tileinfo,tilesize,tilenum::Int)
    tilidx = findidx(tileinfo["tileNo"],tilenum)

    if isa(fldvals,MeshArrays.gcmfaces)
        is3D = length(size(fldvals)) == 3
        if is3D; n3 = size(fldvals)[3]; end
    else
        is3D = length(size(fldvals)) == 2 && size(fldvals)[2] > 1
        if is3D; n3 = size(fldvals)[2]; end
    end

    if is3D # has depth
        tilfld = Array{Any,3}(nothing,0,tilesize[2],n3)
    else
        tilfld = Array{Any,2}(nothing,0,tilesize[2])
    end
    if isa(fldvals,MeshArrays.gcmfaces)
        nfaces = fldvals.nFaces
    elseif isa(fldvals,MeshArray)
        nfaces = fldvals.grid.nFaces
    end
    for iF = 1:nfaces
        if ~isempty(tilidx[iF])
            imin = minimum(tilidx[iF])[1]; imax = maximum(tilidx[iF])[1]
            jmin = minimum(tilidx[iF])[2]; jmax = maximum(tilidx[iF])[2]
            if is3D
                if isa(fldvals,MeshArray)
                    tilfld = [tilfld; cat([fldvals.f[iF,d][imin:imax,jmin:jmax] for d in 1:n3]...,dims=3)]
                else
                    tilfld = [tilfld; fldvals.f[iF][imin:imax,jmin:jmax,:]]
                end
            else
                tilfld = [tilfld; fldvals.f[iF][imin:imax,jmin:jmax]]
            end
        end
    end

    if is3D
        tilfld = reshape(tilfld,tilesize[1],tilesize[2],:)
    else
        tilfld = reshape(tilfld,tilesize[1],tilesize[2])
    end

    return tilfld
end

"""
    gettiles(fldvals,tilenum::Int)

Helper function for retrieving a tile from a gcmfaces struct as a numeric Array along
    with associated latitude and longitude. Not currently exported.
"""
function gettiles(tilfld,tilenum::Int)
    tilesize = tilfld.tilesize

    tilfld = gettile(tilfld.vals,tilfld.tileinfo,tilesize,tilenum)
    tillat = gettile(tilfld.tileinfo["XC"],tilfld.tileinfo,tilfld.tilesize,tilenum)
    tillon = gettile(tilfld.tileinfo["YC"],tilfld.tileinfo,tilfld.tilesize,tilenum)

    return tilfld,tillat,tillon
end

"""
    getindex(v::NCvar,i::Integer)

Retrive an NCvar at time step i.

Currently only works when values is a BinData, NCdata, or array of time steps.
"""
function getindex(v::NCvar,i::Integer)
    NCvar(v.name,v.units,v.dims,v.values[i],v.atts,v.backend)
end

"""
    getindex(d::BinData,i::Integer)

 Return a new BinData object with the ith file, this being time index i.
"""
function getindex(d::BinData,i::Integer) # Gets file at time index i
    if isa(d.fnames,Array)
        newfnames = [d.fnames[i]]
    elseif i == 1
        newfnames = d.fnames
    end

    BinData(newfnames,d.precision,d.iosize)
end

"""
    getindex(d::NCData,i)

Retrieve the ith time step from the NetCDF file.
"""
function getindex(d::NCData,i) # Gets data at time index i
    readncdata(d,i)
end

"""
    readbin(fname::String,prec::Type,iosize::Tuple,fldidx=1)

Read in a binary file as an array as previously done via `MeshArrays.read_bin`
"""
function readbin(fname::String,prec::Type,iosize::Tuple,fldidx=1)
    n3 = 1
    if length(iosize) == 3
        n1,n2,n3 = iosize
    else
        n1,n2 = iosize
    end

    if prec == Float64
        reclen = 8
    else
        reclen = 4
    end

    #if isempty(n3)
    #    n3=Int64(stat(fname).size/n1/n2/reclen)
    #end

    field = Array{prec,1}(undef,(n1*n2*n3))
    fid = open(fname)
    if fldidx != 1
        seek(fid,n1*n2*n3*(fldidx-1)*reclen)
    end
    read!(fid,field)
    close(fid)
    field = hton.(field)

    n3>1 ? s=(n1,n2,n3) : s=(n1,n2)
    return reshape(field,s)

end

"""
    readbin(flddata::BinData,tidx=1)

Read in a binary file as an array as previously done via `MeshArrays.read_bin`
"""
readbin(flddata::BinData,tidx=1) = readbin(flddata.fnames[tidx],flddata.precision,flddata.iosize,flddata.fldidx)
"""
    readncdata(var::NCData,i::Union{Colon,Integer}=:)

Read netcdf file as specified in `NCData` argument. Optional
argument `i` can be used to read a specific records / times.
"""
function readncdata(var::NCData,i::Union{Colon,Integer}=:)
    ds = Dataset(var.fname)
    ndims = length(size(ds[var.varname]))

    if ndims == 1
        values = ds[var.varname][i]
    elseif ndims == 2
        values = ds[var.varname][:,i]
    elseif ndims == 3
        values = ds[var.varname][:,:,i]
    elseif ndims == 4
        values = ds[var.varname][:,:,:,i]
    else
        println("Variables greater than four dimensions not currently supported")
    end
    close(ds)
    return values
end

function readdata(flddata,tidx=1)
    if isa(flddata,BinData)
        res = readbin(flddata,tidx)#readbin(flddata.fnames[tidx],flddata.precision,flddata.iosize,flddata.fldidx)
    elseif isa(flddata,NCData)
        res = readncdata(flddata,tidx)
    end
    return res
end

"""
    addDim(ds::NCDatasets.Dataset,dimvar::NCvar) # NCDatasets

Add a dimension to a `NCDatasets.Dataset`
"""
function addDim(ds::NCDatasets.Dataset,dimvar::NCvar)
    defDim(ds,dimvar.name,dimvar.dims[1])
end

"""
    addDim(dimvar::NCvar)

Add a dimension to a NetCDF file using `NetCDF.jl`
"""
function addDim(dimvar::NCvar) #NetCDF
    NcDim(dimvar.name,collect(dimvar.values),
    atts = merge(Dict(("units" =>dimvar.units)),dimvar.atts),
    unlimited = dimvar.dims==Inf)
end

"""
    addVar(ds::NCDatasets.Dataset,field::NCvar)

Add a variable to a NetCDF file using `NCDatasets.jl`
"""
function addVar(ds::NCDatasets.Dataset,field::NCvar)
    if ~isempty(field.units)
        atts = merge(Dict(("units" =>field.units)),field.atts)
    else
        atts = field.atts
    end
    if isa(field.values,Array)
        if isa(field.values[1],Array)
            prec = typeof(field.values[1][1])
        else
            prec = typeof(field.values[1])
        end
        fieldvar = defVar(ds,field.name,prec,
                            tuple([f.name for f in field.dims]...),attrib=atts)
    else
        fieldvar = defVar(ds,field.name,field.values.precision,
                            tuple([f.name for f in field.dims]...),attrib=atts)
    end
end

"""
    addVar(field::NCvar,dimlist::Array{NetCDF.NcDim})

Add a variable with dimensions dimlist to a NetCDF file using `NetCDF.jl`
"""
function addVar(field::NCvar,dimlist::Array{NetCDF.NcDim})
    if ~isempty(field.units)
        atts = merge(Dict(("units" =>field.units)),field.atts)
    else
        atts = field.atts
    end
    fieldvar = NcVar(field.name,dimlist,
                        atts = atts,
                        t = field.values.precision)
end

"""
    addVar(field::NCvar})

Add a variable and its dimensions to a NetCDF file using `NetCDF.jl`
"""
function addVar(field::NCvar)
    dimlist = addDim.(field.dims)
    if ~isempty(field.units)
        attributes = merge(Dict(("units" =>field.units)),field.atts)
    else
        attributes = field.atts
    end
    fieldvar = NcVar(field.name,dimlist,
        atts = merge(Dict(("units" =>field.units)),field.atts))#,
        #t = field.values.precision)
end

"""
    addData(v::Union{NCDatasets.CFVariable,NetCDF.NcVar},var::NCvar)

Fill variable with data in netcdf file using either `NCDatasets.jl`
or `NetCDF.jl`
"""
function addData(v::Union{NCDatasets.CFVariable,NetCDF.NcVar,Array},var::NCvar,startidx=1)
    isBinData = isa(var.values,BinData)
    isNCData = isa(var.values,NCData)
    isTileData = isa(var.values,TileData)
    if isTileData
        isBinData = isa(var.values.vals,BinData)
    end
    nsteps = getnsteps(var.values)
    ndims = length(var.dims)
    if hastimedim(var)
        ndims = ndims-1
    end
    if isBinData || isNCData || isTileData || isa(var.values[1],Array) # Binary files or array of timesteps

        if isBinData && ~ isTileData
            if isa(var.values.fnames,Array)
                fnames = var.values.fnames
            else
                fnames = [var.values.fnames]
            end
        end

        for i = startidx:nsteps
            if isTileData
                writetiles.(v,Ref(var),1:var.values.numtiles,Ref(i))
            else
                if isBinData || isNCData
                    v0 = readdata(var.values,i)
                else
                    v0 = var.values[i]
                end

                if ndims == 1
                    v[:,i] = v0
                elseif ndims == 2
                    v[:,:,i] = v0
                elseif ndims == 3
                    v[:,:,:,i] = v0
                end
            end
        end

    elseif isa(var.values[1],Number) # Single array of data- just insert it
        ndims = length(size(var.values))
        if ndims == 1
            v[:] = var.values
        elseif ndims == 2
            v[:,:] = var.values
        elseif ndims == 3
            v[:,:,:] = var.values
        elseif ndims == 4
            v[:,:,:,:] = var.values
        end
    else
        print("Unrecognized values")
    end

end

"""
    writetiles(v,var,tilenum,timeidx=1)

Helper function for writing a tile to a NetCDF file.
"""
function writetiles(v,var,tilenum,timeidx=1)
    if isa(v,Array)
        v = v[findfirst(isequal(var.name),name.(v))]
    end
    tileinfo = var.values.tileinfo; tilesize = var.values.tilesize; grid = tileinfo["XC"].grid
    if isa(var.values.vals,BinData)
        iosize = var.values.vals.iosize
        prec = var.values.precision
        if isa(var.values.vals.fnames,Array)
            fnames = var.values.vals.fnames
        else
            fnames = [var.values.vals.fnames]
        end
        if length(iosize) == 3
            f = Array{Array{prec,2},2}(undef,grid.nFaces,iosize[3])
            exarray = MeshArray(grid,f)
        else
            exarray = tileinfo["XC"]
        end
        v0 = read(readbin(fnames[timeidx],prec,iosize,
                                        var.values.vals.fldidx),exarray)
    else
        if isa(var.values.vals,MeshArray) || isa(var.values.vals,MeshArrays.gcmfaces)
            v0 = var.values.vals
        else
            v0 = var.values.vals[timeidx]
        end
    end

    v0 = gettile(v0,tileinfo,tilesize,tilenum)
    numdims = length(size(v0))
    if numdims == 1
        v[:,timeidx] = v0
    elseif numdims == 2
        v[:,:,timeidx] = v0
    elseif numdims == 3
        v[:,:,:,timeidx] = v0
    end

end

"""
    addDimData(v::Union{NCDatasets.CFVariable,NetCDF.NcVar,Array},var::NCvar)

Add dimension data to predefined dimensions in a NetCDF file.
"""
function addDimData(ds,dimvar::NCvar)
    atts = merge(Dict(("units" =>dimvar.units)),dimvar.atts)
    defVar(ds,dimvar.name,dimvar.values,(dimvar.name,),attrib=atts)
end

"""
    createfile(filename, field::Union{NCvar,Dict{String,NCvar}}, README;
                fillval=NaN, missval=NaN, itile=1, ntile=1)

Create NetCDF file and add variable + dimension definitions
using either `NCDatasets.jl` or `NetCDF.jl`
"""
function createfile(filename, field::Union{NCvar,Dict}, README;
                    fillval=NaN, missval=NaN, itile=1, ntile=1)

    if isa(field,Dict)

        dims = unique(vcat([field[v].dims for v in keys(field)]...))
        dims = filter( d -> isa(d,NCvar),dims)
        field = collect(values(field))

        fieldnames = getfield.(field[findall(hastimedim.(field))],:name)
        backend = field[1].backend
    else
        dims = field.dims
        fieldnames = field.name
        backend = field.backend
        if backend == NCDatasets
            field = [field]
        end
    end

    if isa(fieldnames,Array)
        fieldnamestring = join(fieldnames,",")
    else
        fieldnamestring = fieldnames
    end

    file_atts = vcat(["date" => Dates.format(today(),"dd-u-yyyy"),
    "Conventions" => "CF-1.6",
    "description" => fieldnamestring*" -- "*README[1]],
    [string(Char(65+(i-2))) => README[i] for i in 2:length(README)],
    [#string(Char(65+length(README)-1)) => "file created using NCTiles.jl",
    "_FillValue" => fillval,
    "missing_value" => missval,
    "itile" => itile,
    "ntile" => ntile])

    if backend == NCDatasets
        ds =  Dataset(filename,"c",attrib=file_atts)
        dimlist = addDim.(Ref(ds),dims)
        fieldvar = addVar.(Ref(ds), field)
        if length(fieldvar) == 1
            fieldvar = fieldvar[1]
        end
        return ds,fieldvar,dimlist
    elseif backend == NetCDF # Not supporting multiple variables right now
        #= for i = 1:length(field)
            f = field[1];
            if isa(f.values,BinData) || isa(f.values,NCData)
                prec = f.values.precision
            else
                prec = eltype(f.values)
            end
            nccreate(filename,f.name,Tuple(vcat([[d.name; [collect(d.values)]; d.atts] for d in f.dims]...))...,atts = f.atts,t = prec)
        end
        #ncclose(filename)
        ds = NetCDF.open(filename)
        if length(field) == 1
            fieldvar = ds[field[1].name]
        else
            fieldvar = getkey.(Ref(ds),[f.name for f in field])
        end

        return ds,fieldvar,collect(values(ds.dim)) =#
        #fieldvar = Array{NcVar,1}
        dimlist = addDim.(dims)
        #for f in field
        #    fieldvar = vcat(fieldvar,addVar(f))
        #end
        field = field
        fieldvar = addVar(field)#[addVar(f) for f in field]
        file_atts = Dict(file_atts)
        #nccreate(filename,)
        return NetCDF.create(filename,fieldvar, gatts = file_atts),fieldvar,dimlist
    end

end

"""
    readncfile(fname,backend::Module=NCDatasets)

Read in a NetCDF file and return variables/dimensions as `NCvar` structs, and
    file attributes as `Dict`. Large variables/dimensions are not loaded into
    memory. This can use either `NCDatasets.jl` or `NetCDF.jl`
"""
function readncfile(fname,backend::Module=NCDatasets)

    ds = Dataset(fname)

    dims = Dict{AbstractString,NCvar}()
    vars = Dict{AbstractString,NCvar}()

    for k in keys(ds)
        if ~haskey(dims,k)
            k_units = get(ds[k].attrib,"units","")
            k_atts = Dict(ds[k].attrib)
            if length(dimnames(ds[k])) == 1 && k == dimnames(ds[k])[1] # this variable is also a dimension
                k_dims = length(ds[k])
                if length(size(ds[k])) == 1
                    k_values = ds[k][:]
                elseif length(size(ds[k])) == 2
                    k_values = ds[k][:,:]
                else
                    k_values = NCData(fname,k,backend,typeof(ds[k].var[:][1])) # Pointer to this variable in the file
                end
                if !isa(k_values,NCData) && !any(isa.(k_values,Ref(Missing)))
                    k_values = typeof(k_values[1]).(k_values)
                end
                dims[k] = NCvar(k,k_units,k_dims,k_values,k_atts,backend)
            else # Not a dimension
                k_dims = Array{NCvar,1}()
                for d in dimnames(ds[k])
                    #global k_dims,dims
                    if !haskey(dims,d) # Add the dimension
                        d_units = get(ds[d].attrib,"units","")
                        d_dims = length(ds[d])
                        d_values = NCData(fname,d,backend,typeof(ds[d][1])) # Pointer to this variable in the file
                        d_atts = Dict(ds[d].attrib)
                        dims[d] = NCvar(d,d_units,d_dims,d_values,d_atts,backend)
                    end
                    k_dims = cat(k_dims,dims[d],dims=1)
                end
                hasTimeDim = hastimedim(k_dims)
                if  ~hasTimeDim && length(k_dims) == 1# extra variable, not main field, just load the data now
                    k_values = ds[k][:]
                elseif ~hasTimeDim && length(k_dims) == 2
                    k_values = ds[k][:,:]
                elseif ~hasTimeDim && length(k_dims) == 3
                    k_values = ds[k][:,:,:]
                elseif hasTimeDim
                    k_values = NCData(fname,k,backend,eltype(ds[k].var[:])) # Pointer to this variable in the file
                end

                if !hasTimeDim && !any(isa.(k_values,Ref(Missing)))
                    k_values = typeof(k_values[1]).(k_values)
                end

                vars[k] = NCvar(k,k_units,k_dims,k_values,k_atts,backend)
            end
        end
    end

    atts = ["_FillValue","missing_value","itile","ntile"]
    fileatts = Dict()
    for a in atts
        if haskey(ds.attrib,a)
            fileatts[a] = ds.attrib[a]
        end
    end
    close(ds)
    return vars,dims,fileatts
end

"""
    istimedim(d::NCvar)

Helper function: determines whether d is a time dimension. Not exported.
"""
function istimedim(d::NCvar)

    timeUnits = ["minutes","seconds","hours","days","minute","second","hour","day"]
    return any(occursin.(timeUnits,Ref(lowercase(d.units))))

end

"""
    findtimedim(v::NCvar)

Helper function: finds which dimension is a time dimension, if any. Not exported.
"""
function findtimedim(v::NCvar)
    return findall(istimedim.(v.dims))[1]
end

"""
    hastimedim(v::NCvar)

Helper function: determines whether a variable has a time dimension. Not exported.
"""
function hastimedim(v::NCvar)
    ncvardim = isa.(v.dims,NCvar)
    return any(ncvardim) && any(istimedim.(v.dims[ncvardim]))
end

"""
    hastimedimdims::Array{NCvar})

Helper function: determines whether an array of dimensions has a time dimension. Not
    exported.
"""
function hastimedim(dims::Array{NCvar})
    return any(istimedim.(dims))
end

"""
    parsemeta(metafile)

Parse out an `MITgcm` metadata file and return a `Dict` of fields in the file.
"""
function parsemeta(metafile)

    meta = read(metafile,String)
    meta = split(meta,";\n")
    meta = meta[isempty.(meta).==false]
    meta = replace.(meta,Ref(",\n"=>";"))
    meta = replace.(meta,Ref("\n"=>""))
    meta = replace.(meta,Ref("}"=>"]"))
    meta = replace.(meta,Ref("{"=>"["))
    #meta = replace.(meta,Ref(" "=>""))
    meta = replace.(meta,Ref("'"=>"\""))
    meta = replace.(meta,Ref(";]"=>"]"))
    meta = replace.(meta,Ref(","=>" "))

    meta = split.(meta,"=")
    meta = [[replace(x[1]," "=>"") x[2]] for x in meta]

    metaDict = Dict{String,Any}(m[1] => m[2] for m in meta)

    for k in keys(metaDict)
        val = eval(Meta.parse(metaDict[k]))
        if isa(val[1],String)
            val = replace.(val,Ref(" "=>""))
        end
        if length(val) == 1
            val = val[1]
        end
        metaDict[k] = val
    end
    metaDict["dataprec"] = titlecase(metaDict["dataprec"])
    return metaDict

end

"""
    readAvailDiagnosticsLog(fname,fldname)

Get the information for a particular field from the `available_diagnostics.log`
    file (`MITgcm` output).
"""
function readAvailDiagnosticsLog(fname,fldname)
    availdiags = readlines(fname)
    line = availdiags[findall(occursin.(@sprintf("%-8s",fldname),availdiags))[1]]

    line = split(line,'|')
    line = lstrip.(rstrip.(line))

    diagInfo = Dict([
    "diagNum" => parse(Int,line[1]),
    "fldname" => line[2],
    "levs" => parse(Int,line[3]),
    "mate" => line[4],
    "code" => line[5],
    "units" => line[6],
    "title" => line[7]
    ])

end


text:
using StaticArrays
using Plots
using BenchmarkTools

function f(u, p)
    α, β, γ, δ = p

    du1 = α*u[1] - β*u[1]*u[2]
    du2 = -γ*u[2] + δ*u[1]*u[2]

    @SVector [du1, du2]
end

function dopri5!(u_cache, u, p, t_span, f)
    # Values from the Dormand Prince Tableau
    a21  =  (1.0/5.0)
    a31  = (3.0/40.0)
    a32  = (9.0/40.0)
    a41  = (44.0/45.0)
    a42  = (-56.0/15.0)
    a43  = (32.0/9.0)
    a51  = (19372.0/6561.0)
    a52  = (-25360.0/2187.0)
    a53  = (64448.0/6561.0)
    a54  = (-212.0/729.0)
    a61  = (9017.0/3168.0)
    a62  = (-355.0/33.0)
    a63  = (46732.0/5247.0)
    a64  = (49.0/176.0)
    a65  = (-5103.0/18656.0)
    a71  = (35.0/384.0)
    a72  = (0.0)
    a73  = (500.0/1113.0)
    a74  = (125.0/192.0)
    a75  = (-2187.0/6784.0)
    a76  = (11.0/84.0)

    b1   = (35.0/384.0)
    b2   = (0.0)
    b3   = (500.0/1113.0)
    b4   = (125.0/192.0)
    b5   = (-2187.0/6784.0)
    b6   = (11.0/84.0)
    b7   = (0.0)

    dt = 0.25

    @inbounds for i in 2:Int64(t_span/dt)
        k1 = f(u, p)
        k2 = f(u + dt.*(a21.*k1), p)
        k3 = f(u + dt.*(a31.*k1 + a32.*k2), p)
        k4 = f(u + dt.*(a41.*k1 + a42.*k2 + a43.*k3), p)
        k5 = f(u + dt.*(a51.*k1 + a52.*k2 + a53.*k3 + a54.*k4), p)
        k6 = f(u + dt.*(a61.*k1 + a62.*k2 + a63.*k3 + a64.*k4 + a65.*k5), p)
        u = u + dt.*(b1.*k1 + b3.*k3 + b4.*k4 + b5.*k5 + b6.*k6)
        u_cache[:, i] = u
    end
end

p = (1.5, 1.0, 3.0, 1.0)
u₀ = [1.0, 1.0]
u_cache = Array{Float64}(undef, 2, 40)
u_cache[:, 1] = u₀
t_span = 10.0
u_cache

@btime dopri5!(u_cache, u₀, p, t_span, f)

dopri5!(u_cache, u₀, p, t_span, f)

u_cache

plot(0.25:0.25:t_span, u_cache[1, :])
plot!(0.25:0.25:t_span, u_cache[2, :])


text:
#   This file is part of AbstractTensors.jl. It is licensed under the AGPL license
#   Grassmann Copyright (C) 2019 Michael Reed

module AbstractTensors

# universal root Tensor type

"""
    TensorAlgebra{V} <: Number

Universal root tensor type with `Manifold` instance parameter `V`.
"""
abstract type TensorAlgebra{V} <: Number end
Base.@pure istensor(t::T) where T<:TensorAlgebra = true
Base.@pure istensor(t) = false

## Manifold{n}

"""
    Manifold{n} <: TensorAlgebra

Basis parametrization locally homeomorphic to `ℝ^n` product topology.
"""
abstract type Manifold{n} <: TensorAlgebra{n} end
Base.@pure ismanifold(t::T) where T<:Manifold = true
Base.@pure ismanifold(t) = false

"""
    TensorGraded{V,n} <: Manifold{n} <: TensorAlgebra

Graded elements of a `TensorAlgebra` in a `Manifold` topology.
"""
abstract type TensorGraded{V,G} <: Manifold{G} end
const TAG = (:TensorAlgebra,:TensorGraded)
Base.@pure isgraded(t::T) where T<:TensorGraded = true
Base.@pure isgraded(t) = false

# parameters accessible from anywhere

Base.@pure Manifold(::T) where T<:TensorAlgebra{V} where V = V
Base.@pure Manifold(::T) where T<:TensorGraded{V} where V = V
Base.@pure Manifold(V::T) where T<:Manifold = V
Base.@pure Base.parent(V::T) where T<:TensorAlgebra = Manifold(V)

import LinearAlgebra
import LinearAlgebra: UniformScaling, I, rank
import AbstractLattices: ∧, ∨

"""
    rank(::Manifold{n})

Dimensionality `n` of the `Manifold{n}` subspace representation.
"""
Base.@pure LinearAlgebra.rank(::M) where M<:Manifold{n} where n = n

"""
    ndims(t::TensorAlgebra)

Dimensionality of the psuedoscalar, `rank(Manifold(t))` of an element.
"""
Base.@pure Base.ndims(M::T) where T<:TensorAlgebra = rank(Manifold(M))

for (part,G) ∈ ((:scalar,0),(:vector,1),(:bivector,2),(:trivector,3))
    ispart = Symbol(:is,part)
    str = """
    $part(::TensorAlgebra)

Return the $part (rank $G) part of any `TensorAlgebra` element.
    """
    @eval begin
        @doc $str $part
        @inline $part(t::T) where T<:TensorGraded{V} where V = zero(V)
        @inline $part(t::T) where T<:TensorGraded{V,$G} where V = t
        @inline $ispart(t::T) where T<:TensorGraded = rank(t) == $G || iszero(t)
    end
end

"""
    pseudoscalar(::TensorAlgebra)

Return the pseudoscalar (full rank) part of any `TensorAlgebra` element.
"""
@inline volume(t::T) where T<:Manifold = t
@inline volume(t::T) where T<:TensorGraded{V,G} where {V,G} = G == ndims(V) ? t : zero(V)
@inline isvolume(t::T) where T<:TensorGraded = rank(t) == ndims(t) || iszero(t)

"""
    value(::TensorAlgebra)

Returns the internal representation of a `TensorAlgebra` element value.
"""
value(t::T) where T<:Number = t
value(t::T) where T<:AbstractArray = t

"""
    valuetype(t::TensorAlgebra)

Returns type of a `TensorAlgebra` element value's internal representation.
"""
valuetype(::T) where T<:Number = T

function Base.isapprox(a::S,b::T) where {S<:TensorAlgebra,T<:TensorAlgebra}
    rtol = Base.rtoldefault(valuetype(a), valuetype(b), 0)
    LinearAlgebra.norm(Base.:-(a,b))≤rtol*max(LinearAlgebra.norm(a),LinearAlgebra.norm(b))
end
function Base.isapprox(a::S,b::T) where {S<:TensorGraded,T<:TensorGraded}
    Manifold(a)==Manifold(b) && (rank(a)==rank(b) ? AbstractTensors.:≈(norm(a),norm(b)) : (isnull(a) && isnull(b)))
end

# universal vector space interopability, abstract tensor form evaluation, contraction

for X ∈ TAG, Y ∈ TAG
    @eval begin
        @inline interop(op::Function,a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = op(a,b)
        @inline interform(a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = a(b)
        # ^^ identity ^^ | vv union vv #
        @inline function interop(op::Function,a::A,b::B) where {A<:$X,B<:$Y}
            M = Manifold(a) ∪ Manifold(b)
            return op(M(a),M(b))
        end
        @inline function interform(a::A,b::B) where {A<:$X,B<:$Y}
            M = Manifold(a) ∪ Manifold(b)
            return M(a)(M(b))
        end
        @inline ∗(a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = (~a)*b
        @inline ⊛(a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = scalar(contraction(a,b))
        @inline ⨼(a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = contraction(b,a)
        @inline Base.:<(a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = contraction(b,a)
    end
    for op ∈ (:⨽,:(Base.:>),:(Base.:|),:(LinearAlgebra.dot))
        @eval @inline $op(a::A,b::B) where {A<:$X{V},B<:$Y{V}} where V = contraction(a,b)
    end
end

# extended compatibility interface

export TensorAlgebra, Manifold, TensorGraded, Distribution
export istensor, isgraded, isdistribution, ismanifold, rank
export scalar, isscalar, vector, isvector, bivector, isbivector, volume, isvolume
export value, valuetype, interop, interform, involute, unit, even, odd, contraction
export ⊘, ⊖, ⊗, ⊛, ⊙, ⊠, ×, ⨼, ⨽, ⋆, ∗, ⁻¹, ǂ, ₊, ₋, ˣ

# some shared presets

for op ∈ (:(Base.:+),:(Base.:-),:(Base.:*),:⊘,:⊗,:⊛,:∗,:⨼,:⨽,:contraction,:(LinearAlgebra.dot),:(Base.:|),:(Base.:(==)),:(Base.:<),:(Base.:>),:(Base.:<<),:(Base.:>>),:(Base.:>>>),:(Base.div),:(Base.rem),:(Base.:&))
    @eval begin
        @inline $op(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = interop($op,a,b)
        @inline $op(a::A,b::UniformScaling) where A<:TensorAlgebra = $op(a,Manifold(a)(b))
        @inline $op(a::UniformScaling,b::B) where B<:TensorAlgebra = $op(Manifold(b)(a),b)
    end
end

const ⊖ = *
@inline ⋆(t::UniformScaling{T}) where T = T<:Bool ? (t.λ ? 1 : -1) : t.λ
@inline Base.:/(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = a*Base.inv(b)
@inline Base.:/(a::UniformScaling,b::B) where B<:TensorAlgebra = Manifold(b)(a)*Base.inv(b)
@inline Base.:/(a::A,b::UniformScaling) where A<:TensorAlgebra = a*Base.inv(Manifold(a)(b))
@inline Base.:\(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = Base.inv(a)*b
@inline Base.:\(a::UniformScaling,b::B) where B<:TensorAlgebra = Base.inv(Manifold(b)(a))*b
@inline Base.:\(a::A,b::UniformScaling) where A<:TensorAlgebra = Base.inv(a)*Manifold(a)(b)

for op ∈ (:(Base.:+),:(Base.:*))
    @eval $op(t::T) where T<:TensorAlgebra = t
end
for op ∈ (:|,:!), T ∈ (TensorAlgebra,UniformScaling)
    @eval Base.$op(t::$T) = ⋆(t)
end
for op ∈ (:⊙,:⊠,:¬,:⋆)
    @eval function $op end
end
for op ∈ (:scalar,:involute,:even)
    @eval $op(t::T) where T<:Real = t
end
odd(::T) where T<:Real = 0
LinearAlgebra.cross(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = ⋆(∧(a,b))

@inline Base.exp(t::T) where T<:TensorAlgebra = 1+Base.expm1(t)
@inline Base.log(b,t::T) where T<:TensorAlgebra = Base.log(t)/Base.log(b)
@inline Base.:^(b::S,t::T) where {S<:Number,T<:TensorAlgebra} = Base.exp(t*Base.log(b))
@inline Base.:^(a::A,b::UniformScaling) where A<:TensorAlgebra = Base.:^(a,Manifold(a)(b))
@inline Base.:^(a::UniformScaling,b::B) where B<:TensorAlgebra = Base.:^(Manifold(b)(a),b)

for base ∈ (2,10)
    fl,fe = (Symbol(:log,base),Symbol(:exp,base))
    @eval Base.$fl(t::T) where T<:TensorAlgebra = Base.$fl(ℯ)*Base.log(t)
    @eval Base.$fe(t::T) where T<:TensorAlgebra = Base.exp(Base.log($base)*t)
end

@inline Base.cos(t::T) where T<:TensorAlgebra = Base.cosh(Manifold(t)(I)*t)
@inline Base.sin(t::T) where T<:TensorAlgebra = (i=Manifold(t)(I);Base.:/(Base.sinh(i*t),i))
@inline Base.tan(t::T) where T<:TensorAlgebra = Base.:/(Base.sin(t),Base.cos(t))
@inline Base.cot(t::T) where T<:TensorAlgebra = Base.:/(Base.cos(t),Base.sin(t))
@inline Base.sec(t::T) where T<:TensorAlgebra = Base.inv(Base.cos(t))
@inline Base.csc(t::T) where T<:TensorAlgebra = Base.inv(Base.sin(t))
@inline Base.asec(t::T) where T<:TensorAlgebra = Base.acos(Base.inv(t))
@inline Base.acsc(t::T) where T<:TensorAlgebra = Base.asin(Base.inv(t))
@inline Base.sech(t::T) where T<:TensorAlgebra = Base.inv(Base.cosh(t))
@inline Base.csch(t::T) where T<:TensorAlgebra = Base.inv(Base.sinh(t))
@inline Base.asech(t::T) where T<:TensorAlgebra = Base.acosh(Base.inv(t))
@inline Base.acsch(t::T) where T<:TensorAlgebra = Base.asinh(Base.inv(t))
@inline Base.tanh(t::T) where T<:TensorAlgebra = Base.:/(Base.sinh(t),Base.cosh(t))
@inline Base.coth(t::T) where T<:TensorAlgebra = Base.:/(Base.cosh(t),Base.sinh(t))
@inline Base.asinh(t::T) where T<:TensorAlgebra = Base.log(t+Base.sqrt(1+Base.:^(t,2)))
@inline Base.acosh(t::T) where T<:TensorAlgebra = Base.log(t+Base.sqrt(Base.:-(Base.:^(t,2),1)))
@inline Base.atanh(t::T) where T<:TensorAlgebra = Base.:/(Base.:-(Base.log(1+t),Base.log(Base.:-(1,t))),2)
@inline Base.acoth(t::T) where T<:TensorAlgebra = Base.:/(Base.:-(Base.log(t+1),Base.log(Base.:-(t,1))),2)
Base.asin(t::T) where T<:TensorAlgebra = (i=Manifold(t)(I);Base.:-(i)*Base.log(i*t+Base.sqrt(Base.:-(1,Base.:^(t,2)))))
Base.acos(t::T) where T<:TensorAlgebra = (i=Manifold(t)(I);Base.:-(i)*Base.log(t+i*Base.sqrt(Base.:-(1,Base.:^(t,2)))))
Base.atan(t::T) where T<:TensorAlgebra = (i=Manifold(t)(I);Base.:/(Base.:-(i),2)*Base.:-(Base.log(1+i*t),Base.log(Base.:-(1,i*t))))
Base.acot(t::T) where T<:TensorAlgebra = (i=Manifold(t)(I);Base.:/(Base.:-(i),2)*Base.:-(Base.log(Base.:-(t,i)),Base.log(t+i)))
Base.sinc(t::T) where T<:TensorAlgebra = iszero(t) ? one(Manifold(t)) : (x=(1π)*t;Base.:/(Base.sin(x),x))
Base.cosc(t::T) where T<:TensorAlgebra = iszero(t) ? zero(Manifold(t)) : (x=(1π)*t; Base.:-(Base.:/(Base.cos(x),t), Base.:/(sin(x),(x*t))))

# absolute value norm

@inline Base.abs(t::T) where T<:TensorAlgebra = Base.sqrt(Base.abs2(t))
@inline Base.abs2(t::T) where T<:TensorAlgebra = (a=t⋅t; isscalar(a) ? scalar(a) : a)
@inline Base.abs2(t::T) where T<:TensorGraded = contraction(t,t)
@inline norm(z) = LinearAlgebra.norm(z)
@inline LinearAlgebra.norm(t::T) where T<:TensorAlgebra = norm(value(t))
@inline unit(t::T) where T<:TensorAlgebra = Base.:/(t,Base.abs(t))
@inline Base.iszero(t::T) where T<:TensorAlgebra = LinearAlgebra.norm(t) ≈ 0
@inline Base.isone(t::T) where T<:TensorAlgebra = LinearAlgebra.norm(t)≈value(scalar(t))≈1
@inline contraction(a::A,b::B) where {A<:TensorGraded,B<:TensorGraded} = contraction(a,b)

# identity elements

for id ∈ (:zero,:one)
    @eval begin
        @inline Base.$id(t::T) where T<:TensorAlgebra = $id(Manifold(t))
        @inline Base.$id(::Type{T}) where T<:TensorAlgebra{V} where V = $id(V)
        @inline Base.$id(::Type{T}) where T<:TensorGraded{V} where V = $id(V)
    end
end

# postfix operators

struct Postfix{Op} end
@inline Base.:*(t,op::P) where P<:Postfix = op(t)
for op ∈ (:⁻¹,:ǂ,:₊,:₋,:ˣ)
    @eval const $op = $(Postfix{op}())
end
@inline (::Postfix{:⁻¹})(t) = inv(t)
@inline (::Postfix{:ǂ})(t) = conj(t)
@inline (::Postfix{:₊})(t) = even(t)
@inline (::Postfix{:₋})(t) = odd(t)
@inline (::Postfix{:ˣ})(t) = involute(t)

# dispatch

@inline norm(z::Expr) = abs(z)
@inline norm(z::Symbol) = z
Base.@pure isnull(::Expr) = false
Base.@pure isnull(::Symbol) = false
isnull(n) = iszero(n)
signbit(x::Symbol) = false
signbit(x::Expr) = x.head == :call && x.args[1] == :-
-(x) = Base.:-(x)
-(x::Symbol) = :(-$x)

for op ∈ (:conj,:inv,:sqrt,:abs,:exp,:expm1,:log,:log1p,:sin,:cos,:sinh,:cosh,:signbit)
    @eval @inline $op(z) = Base.$op(z)
end

for op ∈ (:/,:-,:^,:≈)
    @eval @inline $op(a,b) = Base.$op(a,b)
end

for T ∈ (Expr,Symbol)
    @eval begin
        ≈(a::$T,b::$T) = a == b
        ≈(a::$T,b) = false
        ≈(a,b::$T) = false
    end
end

for (OP,op) ∈ ((:∏,:*),(:∑,:+))
    @eval begin
        @inline $OP(x...) = Base.$op(x...)
        @inline $OP(x::AbstractVector{T}) where T<:Any = $op(x...)
    end
end

const PROD,SUM,SUB = ∏,∑,-

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RectangleCoverer

A method for covering spatial objects with a minimum axis-aligned
bounding rectangle (or bounding box).
"""
struct RectangleCoverer <: AbstractCoverer end

function cover(domain::AbstractDomain{T,N}, coverer::RectangleCoverer) where {N,T}
  lowerleft  = MVector(ntuple(i->typemax(T), N))
  upperright = MVector(ntuple(i->typemin(T), N))

  x = MVector{N,T}(undef)
  for l in 1:npoints(domain)
    coordinates!(x, domain, l)
    for d in 1:N
      x[d] < lowerleft[d]  && (lowerleft[d]  = x[d])
      x[d] > upperright[d] && (upperright[d] = x[d])
    end
  end

  Rectangle(Tuple(lowerleft), Tuple(upperright))
end

function cover(grid::RegularGrid, coverer::RectangleCoverer)
  lowerleft  = origin(grid)
  upperright = origin(grid) .+ (size(grid) .- 1) .* spacing(grid)

  Rectangle(lowerleft, upperright)
end

"""
    boundbox(object)

Return the minimum axis-aligned bounding rectangle of the spatial `object`.

### Notes

Equivalent to `cover(object, RectangleCoverer())`
"""
boundbox(object::AbstractSpatialObject) = cover(object, RectangleCoverer())


text:
# Autogenerated wrapper script for GTK3_jll for i686-w64-mingw32
export libgailutil3, libgdk3, libgtk3

using Glib_jll
using Cairo_jll
using Pango_jll
using FriBidi_jll
using FreeType2_jll
using gdk_pixbuf_jll
using Libepoxy_jll
using ATK_jll
using HarfBuzz_jll
using xkbcommon_jll
using iso_codes_jll
using Wayland_jll
using Xorg_libXrandr_jll
using Xorg_libX11_jll
using Xorg_libXrender_jll
using Xorg_libXi_jll
using Xorg_libXext_jll
using Xorg_libXcursor_jll
using Xorg_libXdamage_jll
using Xorg_libXfixes_jll
using Xorg_libXcomposite_jll
using Xorg_libXinerama_jll
using Fontconfig_jll
using at_spi2_atk_jll
JLLWrappers.@generate_wrapper_header("GTK3")
JLLWrappers.@declare_library_product(libgailutil3, "libgailutil-3-0.dll")
JLLWrappers.@declare_library_product(libgdk3, "libgdk-3-0.dll")
JLLWrappers.@declare_library_product(libgtk3, "libgtk-3-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll, Cairo_jll, Pango_jll, FriBidi_jll, FreeType2_jll, gdk_pixbuf_jll, Libepoxy_jll, ATK_jll, HarfBuzz_jll, xkbcommon_jll, iso_codes_jll, Wayland_jll, Xorg_libXrandr_jll, Xorg_libX11_jll, Xorg_libXrender_jll, Xorg_libXi_jll, Xorg_libXext_jll, Xorg_libXcursor_jll, Xorg_libXdamage_jll, Xorg_libXfixes_jll, Xorg_libXcomposite_jll, Xorg_libXinerama_jll, Fontconfig_jll, at_spi2_atk_jll)
    JLLWrappers.@init_library_product(
        libgailutil3,
        "bin\\libgailutil-3-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgdk3,
        "bin\\libgdk-3-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgtk3,
        "bin\\libgtk-3-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BinaryBuilder, Pkg
import Pkg.Types: VersionSpec

name = "ModularForms"
version = v"0.1.0"

sources = [GitSource("https://gitlab.com/mraum/ModularForms_jll_src.git",
                     "4d81d41e625fde3e36dc86df5061313a8447cacd")]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ModularForms_jll_src
if [[ ${target} == *musl* ]]; then
   export CFLAGS=-D_GNU_SOURCE
fi
./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix --with-arb=$prefix
make -j${nproc}
make install LIBDIR=$(basename ${libdir})
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=Sys.iswindows)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmodularforms_jll", :libmodularforms_jll)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Arb_jll"), compat = "~200.2200.0"),
    Dependency(PackageSpec(name="FLINT_jll"), compat = "^200.800.401")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="^1.6", preferred_gcc_version=v"9")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
# -------------------------------------------------------------------
#* Distillation Column model with:
#*
#*	- NumberOfTrays like tray;
#*	- a vessel in the bottom of column;
#*	- a splitter who separate the bottom product and the stream to reboiler;
#*	- steady state reboiler (thermosyphon);
#*	- a dynamic condenser without subcooling;
#*	- a splitter which separate reflux and distillate;
#*	- a pump in reflux stream.
#*
#* ------------------------------------------------------------------
type User_Distillation_thermosyphon_cond
	User_Distillation_thermosyphon_cond()=begin
		new(
			User_Section_ColumnBasic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["on", "off"],
				:Default=>"on",
				:Hidden=>true
			)),
			condenser(),
			splitter(),
			pump(),
			tank(),
			splitter(),
			reboilerSteady(),
			DanaReal(),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Condenser",
				:Hidden=>true
			)),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Reboiler",
				:Hidden=>true
			)),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Bottom Vessel",
				:Hidden=>true
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream From Top Condenser",
				:PosX=>0.73,
				:PosY=>0
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Vapour outlet stream From Top Condenser",
				:Hidden=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream From Top Splitter",
				:PosX=>1,
				:PosY=>0.45
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Liquid outlet stream From Top Splitter",
				:Hidden=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream From Bottom Splitter",
				:PosX=>1,
				:PosY=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Liquid outlet stream From Reboiler",
				:Hidden=>true
			)),
			[
				:(ConnectorCondenserVout.T = VapourDistillate.T),
				:(ConnectorCondenserVout.P = VapourDistillate.P),
				:(ConnectorCondenserVout.F = VapourDistillate.F),
				:(ConnectorCondenserVout.z = VapourDistillate.z),
				:(ConnectorSplitterTop.T = LiquidDistillate.T),
				:(ConnectorSplitterTop.P = LiquidDistillate.P),
				:(ConnectorSplitterTop.F = LiquidDistillate.F),
				:(ConnectorSplitterTop.z = LiquidDistillate.z),
				:(ConnectorSplitterBottom.T = BottomProduct.T),
				:(ConnectorSplitterBottom.P = BottomProduct.P),
				:(ConnectorSplitterBottom.F = BottomProduct.F),
				:(ConnectorSplitterBottom.z = BottomProduct.z),
				:(CondenserUnity.InletV.F*trays(1).vV = alfaTopo * _base_1.Ah * sqrt(2*(trays(1).OutletV.P - CondenserUnity.OutletL.P + 1e-8 * "atm") / (_base_1.alfa*trays(1).rhoV))),
				:(CondenserUnity.InletV.F = 0 * "mol/s"),
			],
			[
				"","","","","","","","","","","","","","",
			],
			[:CondenserVapourFlow,],
			[:CondenserUnity,:SplitterTop,:PumpUnity,:BottomVessel,:SplitterBottom,:ReboilerUnity,:alfaTopo,:HeatToCondenser,:HeatToReboiler,:HeatToBottomVessel,:VapourDistillate,:ConnectorCondenserVout,:LiquidDistillate,:ConnectorSplitterTop,:BottomProduct,:ConnectorSplitterBottom,]
		)
	end
	_base_1::User_Section_ColumnBasic
	CondenserVapourFlow::DanaSwitcher
	CondenserUnity::condenser
	SplitterTop::splitter
	PumpUnity::pump
	BottomVessel::tank
	SplitterBottom::splitter
	ReboilerUnity::reboilerSteady
	alfaTopo::DanaReal
	HeatToCondenser::energy_stream
	HeatToReboiler::energy_stream
	HeatToBottomVessel::energy_stream
	VapourDistillate::vapour_stream
	ConnectorCondenserVout::stream
	LiquidDistillate::liquid_stream
	ConnectorSplitterTop::stream
	BottomProduct::liquid_stream
	ConnectorSplitterBottom::stream
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export User_Distillation_thermosyphon_cond
function setEquationFlow(in::User_Distillation_thermosyphon_cond)
	# Condenser Connector Equations
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	# Top Splitter Connector Equations
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	# Bottom Splitter Connector Equations
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	let switch=CondenserVapourFlow
		if CondenserVapourFlow==CondenserUnity.InletV.F < 1e-6 * "kmol/h"
			set(switch,"off")
		end
		if CondenserVapourFlow==trays(1).OutletV.P > CondenserUnity.OutletL.P + 1e-1 * "atm"
			set(switch,"on")
		end
		if switch=="on"
			addEquation(13)
		elseif switch=="off"
			addEquation(14)
		end
	end
end
function atributes(in::User_Distillation_thermosyphon_cond,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/DistillationThermosyphonCond"
	fields[:Brief]="Model of a distillation column with dynamic condenser and steady reboiler."
	fields[:Info]="== Specify ==
* the feed stream of each tray (Inlet);
* the Murphree eficiency for each tray Emv;
* the pump head;
* the condenser vapor outlet flow (OutletV.F);
* the heat supllied in bottom tank;
* the heat supllied in condenser and reboiler;
* the Outlet1 flow in the bottom splitter (spbottom.Outlet1.F) that corresponds to the bottom product;

== Initial Conditions ==
* the trays temperature (OutletL.T);
* the trays liquid level (Level) OR the trays liquid flow (OutletL.F);
* (NoComps - 1) OutletL (OR OutletV) compositions for each tray;

* the condenser temperature (OutletL.T);
* the condenser liquid level (Level);
* (NoComps - 1) OutletL (OR OutletV) compositions;

* the bottom tank temperature (OutletL.T);
* the bottom tank liquid level (Level);
* (NoComps - 1) OutletL (OR OutletV) compositions.
"
	drive!(fields,_)
	return fields
end
User_Distillation_thermosyphon_cond(_::Dict{Symbol,Any})=begin
	newModel=User_Distillation_thermosyphon_cond()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
    ReplitDatabaseCore

The core functions which access the database.
This implementation accesses the database using `curl` via shell calls.
The `ReplitDB` interface type depends on these functions.

## Exports

- `get(key; url)`
- `set!(key, value; url)`
- `list([prefix]; url)`
- `delete!(key; url)`
"""
module ReplitDatabaseCore
export delete!, get, list, set!

using URIs: escapeuri, unescapeuri

"""
    delete!(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Deletes the value associated with the given key in the Repl.it database.
"""
function delete!(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])
    success(`curl -s -XDELETE $(url * "/" * escapeuri(key))`)
end

"""
    get(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Get the value associated with the given key in the Repl.it database.
"""
function get(key::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])
    readchomp(`curl -s $(url * "/" * escapeuri(key))`)
end

"""
    list(; url::AbstractString=ENV["REPLIT_DB_URL"])
    list(prefix::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Lists all of the keys in the Repl.it database, or all of the keys starting with `prefix` if specified.
"""
function list(prefix=""::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])::Vector{String}
    escapedurl = "$(url)?prefix=$(escapeuri(prefix))&encode=true"
    result = readchomp(`curl -s "$(escapedurl)"`)
    isempty(result) ? String[] : unescapeuri.(split(result, '\n'))
end

"""
    set!(key::AbstractString, value::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])

Set the given key/value pair in the Repl.it database.
"""
function set!(key::AbstractString, value::AbstractString; url::AbstractString=ENV["REPLIT_DB_URL"])
    success(`curl -s $(url) -d $(key * "=" * value)`)
end

end # module


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_volumedatarepresentation

"""
    vtk_volumedatarepresentation(;kwargs...)

A VolumeDataRepresentation component.
VolumeDataRepresentation expect the following set of properties
  - dimensions: [10, 20, 5]
  - spacing: [1, 1, 1]
  - origin: [0, 0, 0]
  - rgb: [...]
  - rgba: [...]
  - scalars: [...]
  - scalarsType: Float32Array
Keyword arguments:
- `id` (String; optional): The ID used to identify this component.
- `colorDataRange` (Array of Reals | String; optional): Data range use for the colorMap
- `colorMapPreset` (String; optional): Preset name for the lookup table color map
- `controllerSize` (Array of Reals; optional): Controller size in pixels
- `dimensions` (Array of Reals; optional): Number of points along x, y, z
- `mapper` (Dict; optional): Properties to set to the mapper
- `origin` (Array of Reals; optional): World coordinate of the lower left corner of your vtkImageData (i=0, j=0, k=0).
- `property` (Dict; optional): Properties to set to the volume.property
- `rescaleColorMap` (Bool; optional): Use opacity range to rescale color map
- `rgb` (Array of Reals; optional): Use RGB values to attach to the points/vertex
- `rgba` (Array of Reals; optional): Use RGBA values to attach to the points/vertex
- `scalars` (Array of Reals; optional): Field values to attach to the points
- `scalarsType` (String; optional): Types of numbers provided in scalars
- `spacing` (Array of Reals; optional): Spacing along x, y, z between points in world coordinates
- `volume` (Dict; optional): Properties to set to the volume
- `volumeController` (Bool; optional): Show volumeController
"""
function vtk_volumedatarepresentation(; kwargs...)
        available_props = Symbol[:id, :colorDataRange, :colorMapPreset, :controllerSize, :dimensions, :mapper, :origin, :property, :rescaleColorMap, :rgb, :rgba, :scalars, :scalarsType, :spacing, :volume, :volumeController]
        wild_props = Symbol[]
        return Component("vtk_volumedatarepresentation", "VolumeDataRepresentation", "dash_vtk", available_props, wild_props; kwargs...)
end

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libmicrohttpd"
version = v"0.9.73"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/libmicrohttpd/libmicrohttpd-0.9.73.tar.gz", "a37b2f1b88fd1bfe74109586be463a434d34e773530fc2a74364cfcf734c032e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libmicrohttpd-0.9.73/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmicrohttpd", :libmicrohttpd)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GnuTLS_jll", uuid="0951126a-58fd-58f1-b5b3-b08c7c4a876d"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

# methods with no rounding required
-(a::AbstractFloat, ::RoundingMode) = -a
abs(a::Float64, ::RoundingMode) = abs(a)

# tightly rounded methods from FastRounding.jl
for (f, ff) in zip((:+, :-, :*, :/), (:add, :sub, :mul, :div))
    ff = Symbol(ff, "_round")
    for T in (Float32, Float64)
        @eval $f(a::$T, b::$T, r::RoundingMode) = $ff(a, b, r)
    end
end

for f in (:sqrt, :inv), T in (Float32, Float64)
    ff = Symbol(f, "_round")
    @eval $f(a::$T, r::RoundingMode) = $ff(a, r)
end

# for these functions use CRlibm
for f in (:exp, :expm1, :log, :log1p, :log2, :log10, :sin, :cos, :tan, :asin,
          :acos, :atan, :sinh, :cosh, :sinpi, :cospi, :tanpi, :atanpi)
    @eval $f(a::Float64, r::RoundingMode) = CRlibm.$f(a, r)
end

atan_inner(y, x, r) = atan(/(y, x, r), r)

function atan(y::Float64, x::Float64, r::RoundingMode)
    if x > 0
        return atan_inner(y, x, r)
    elseif x < 0
        if ifelse(r == RoundUp, y >= 0, y > 0)
            p = pi_interval(Float64)
            return atan_inner(y, x, r) + ifelse(r == RoundUp, p.hi, p.lo)
        else
            p = pi_interval(Float64)
            return atan_inner(y, x, r) + ifelse(r == RoundUp, -p.lo, -p.hi)
        end
    else
        if y > 0
            p2 = half_pi(Float64)
            return ifelse(r == RoundUp, p2.hi, p2.lo)
        elseif y < 0
            p2 = half_pi(Float64)
            return ifelse(r == RoundUp, -p2.lo, -p2.hi)
        else
            return ifelse(r == RoundUp, zero(Float64), -zero(Float64))
        end
    end
end

function tan(x::Float64, ::RoundingMode{:FromZero})
    val_u = tan(x, RoundUp)
    val_d = tan(x, RoundDown)
    return ifelse(abs(val_u) < abs(val_d), val_d, val_u)
end

cot(x::Float64, r::RoundingMode{:Up}) = inv(tan(x, RoundToZero), r)
cot(x::Float64, r::RoundingMode{:Down}) = inv(tan(x, RoundFromZero), r)

# operations for BigFloat
for f in (:exp, :expm1, :log, :log1p, :log2, :log10, :sin, :cos, :tan, :asin,
          :acos, :atan, :sinh, :cosh, :sinpi, :cospi, :tanpi, :atanpi, :inv,
          :sqrt, :-, :tanh, :asinh, :acosh, :atanh, :cot, :cbrt, :exp2, :exp10)
    @eval function $f(a::BigFloat, r::RoundingMode)
        setrounding(BigFloat, r) do
            $f(a)
        end
    end
end

for f in (:+, :-, :*, :/, :atan, :^)
    @eval function $f(a::BigFloat, b::Real, r::RoundingMode)
        setrounding(BigFloat, r) do
            $f(a, b)
        end
    end
end