text:
mutable struct AtomicSpeciesVars
    # species files path
    sppath::String
    # species filenames
    spfname::Vector{String}
    # species name
    spname::Vector{String}
    # species symbol
    spsymb::Vector{String}
    # species nuclear charge
    spzn::Vector{Float64}
    # ptnucl is .true. if the nuclei are to be treated as point charges, if .false.
    # the nuclei have a finite spherical distribution
    ptnucl::Bool
    # nuclear radius
    rnucl::Vector{Float64}
    # nuclear volume
    volnucl::Vector{Float64}
    # number of radial mesh points to nuclear radius
    nrnucl::Vector{Int64}
    # number of coarse radial mesh points to nuclear radius
    nrcnucl::Vector{Int64}
    # nuclear Coulomb potential
    vcln::Array{Float64,2}
    # species electronic charge
    spze::Vector{Float64}
    # species mass
    spmass::Vector{Float64}
    # smallest radial point for each species
    rminsp::Vector{Float64}
    # effective infinity for species
    rmaxsp::Vector{Float64}
    # number of radial points to effective infinity for each species
    nrsp::Vector{Int64}
    # maximum nrsp over all the species
    nrspmax::Int64
    # maximum allowed states for each species
    maxstsp::Int64 # parameter: 40
    # number of states for each species
    nstsp::Vector{Int64}
    # maximum nstsp over all the species
    nstspmax::Int64
    # core-valence cut-off energy for species file generation
    ecvcut::Float64
    # semi-core-valence cut-off energy for species file generation
    esccut::Float64
    # state principle quantum number for each species
    nsp::Array{Int64,2} # (maxstsp,maxspecies)
    # state l value for each species
    lsp::Array{Int64,2} # (maxstsp,maxspecies)
    # state k value for each species
    ksp::Array{Int64,2} #(maxstsp,maxspecies)
    # spcore is .true. if species state is core
    spcore::Array{Bool,2} #(maxstsp,maxspecies)
    # total number of core states
    nstcr::Int64
    # state eigenvalue for each species
    evalsp::Array{Float64,2} # (maxstsp,maxspecies)
    # state occupancy for each species
    occsp::Array{Float64,2} # (maxstsp,maxspecies)
    # species radial mesh to effective infinity
    rsp::Array{Float64,2}
    # r^l on radial mesh to muffin-tin radius
    rlsp::Array{Float64,3}
    # species charge density
    rhosp::Array{Float64,2}
    # species self-consistent potential
    vrsp::Array{Float64,2}
    # exchange-correlation type for atomic species (the converged ground-state of
    # the crystal does not depend on this choice)
    xctsp::Tuple{Int64,Int64,Int64}
end

function AtomicSpeciesVars( Nspecies::Int64 )

    sppath = ""
    spfname = Vector{String}(undef, Nspecies)
    spname = Vector{String}(undef, Nspecies)
    spsymb = Vector{String}(undef, Nspecies)

    spzn = zeros(Float64, Nspecies)

    ptnucl = true

    rnucl = zeros(Float64, Nspecies)

    volnucl = zeros(Float64, Nspecies)

    nrnucl = zeros(Int64, Nspecies)
    nrcnucl = zeros(Int64, Nspecies)

    vcln = zeros(Float64,1,1) # XXX will be properly setup later

    spze = zeros(Float64, Nspecies)
    spmass = zeros(Float64, Nspecies)

    rminsp = zeros(Float64, Nspecies)
    rmaxsp = zeros(Float64, Nspecies)
    nrsp = zeros(Int64, Nspecies)
    nrspmax = 0
    maxstsp = 40
    nstsp = zeros(Int64, Nspecies)
    nstspmax = 0

    ecvcut = 0.0
    esccut = 0.0

    nsp = zeros(Int64,maxstsp,Nspecies)
    lsp = zeros(Int64,maxstsp,Nspecies)
    ksp = zeros(Int64,maxstsp,Nspecies)
    spcore = zeros(Bool,maxstsp,Nspecies)
    nstcr = 0
    evalsp = zeros(Float64,maxstsp,Nspecies)
    occsp  = zeros(Float64,maxstsp,Nspecies)

    # XXX Will be setup properly later
    rsp = zeros(Float64,1,1)
    rlsp  = zeros(Float64,1,1,1)
    rhosp = zeros(Float64,1,1)
    vrsp  = zeros(Float64,1,1)
    xctsp  = (3,0,0)

    return AtomicSpeciesVars(
        sppath, spfname, spname, spsymb, spzn,
        ptnucl, rnucl, volnucl, nrnucl, nrcnucl,
        vcln, spze, spmass, rminsp, rmaxsp, nrsp,
        nrspmax, maxstsp, nstsp, nstspmax, ecvcut, esccut,
        nsp, lsp, ksp, spcore, nstcr, evalsp, occsp, rsp, rlsp, rhosp, vrsp, xctsp
    )
end

# Setting up vcln
function init_nuclear_pot!( atsp_vars::AtomicSpeciesVars )
    # spherical harmonic for l=m=0
    y00=0.28209479177387814347

    rsp = atsp_vars.rsp
    nrsp = atsp_vars.nrsp
    nspecies = size(nrsp,1)
    ptnucl = atsp_vars.ptnucl
    nrspmax = atsp_vars.nrspmax
    spzn = atsp_vars.spzn

    # determine the nuclear Coulomb potential
    atsp_vars.vcln = zeros(Float64,nrspmax,nspecies)
    vcln = atsp_vars.vcln
    t1 = 1.0/y00
    for is in 1:nspecies
        nr = nrsp[is]
        @views potnucl!(ptnucl, nr, rsp[:,is], spzn[is], vcln[:,is])
        vcln[1:nr,is] = t1*vcln[1:nr,is]
    end
    return
end


text:
using Distributed

nprocs()
addprocs(4)

@distributed for i=1:5
    println("Foo $i")
end

ext: jl

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
x, y = 0, 0

function move_to_border(f, r, side, dif = 0)
    res = 0
    while !isborder(r, side)
        f(r)
        move!(r, side)
        res += dif
    end
    return res
end

function fill_line(r)
    move_to_border(putmarker!, r, Ost)
    putmarker!(r)
    f(x) = nothing
    move_to_border(f, r, West)
end

function fill_all(r)
    global x, y
    f(l) = nothing
    x += move_to_border(f, r, Nord, -1)
    y += move_to_border(f, r, West, -1)
    x += move_to_border(fill_line, r, Sud, 1)
    fill_line(r)
end

fill_all(r)

for i in (1:abs(x))
    move!(r, Nord)
end
for i in (1:abs(y))
    move!(r, Ost)
end
show(r)


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

"""
    NormBall(radius, [metric])

A norm ball with `radius` and `metric`.
Default metric is `Euclidean()`.
"""
struct NormBall{T,M} <: MetricBall
  radius::T
  metric::M
end

function NormBall(radius)
  metric = Euclidean()
  NormBall{typeof(radius),typeof(metric)}(radius, metric)
end

metric(ball::NormBall) = ball.metric

"""
    radius(ball)

Return the radius of the norm `ball`.
"""
radius(ball::NormBall) = ball.radius

function Base.show(io::IO, ball::NormBall{T}) where {T}
  r, d = ball.radius, ball.metric
  print(io, "NormBall{$T}($r, $d)")
end


text:
function iskaprekar(n::Integer)
    n == 1 && return true
    test(a, b) = n == a + b && b ≠ 0
    str = string(n^2)
    any(test(parse(Int, str[1:i]), parse(Int, str[i+1:end])) for i = 1:length(str)-1)
end

@show filter(iskaprekar, 1:10000)
@show count(iskaprekar, 1:10000)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Useful for CDD, LRS and matrix representations

# HyperPlanes and HalfSpaces mixed
abstract type MixedHRep{N, T} <: HRepresentation{N, T} end
# SymPoints and Points mixed and/or Lines and Rays mixed
abstract type MixedVRep{N, T} <: VRepresentation{N, T} end

const MixedRep{N, T} = Union{MixedHRep{N, T}, MixedVRep{N, T}}

# This could be defined for all types Indices but it is not recommended for reps other than Mixed to use this implementation of length as it is inefficient
# a MethodError is more helpful than a hidden inefficiency
function mixedlength(idxs::Indices)
    count = 0
    for idx in idxs
        count += 1
    end
    count
end
Base.length(idxs::Indices{N, T, ElemT, <:MixedRep{N, T}}) where {N, T, ElemT} = mixedlength(idxs)

function mixednext(rep::MixedRep{N, T}, idx::IdxT) where {N, T, ElemT, IdxT<:Index{N, T, ElemT}}
    idx = IdxT(idx.value+1)
    while !done(Indices{N, T, ElemT}(rep), idx) && !isvalid(rep, idx)
        idx = IdxT(idx.value+1)
    end
    idx
end
Base.start(idx::Indices{N, T, ElemT, <:MixedRep{N, T}}) where {N, T, ElemT} = mixednext(idx.rep, Index{N, T, ElemT}(0))
nextindex(rep::MixedRep{N, T}, idx::Index{N, T}) where {N, T} = mixednext(rep, idx)


text:
function get_b1_bot(obj::NeQuickParameters)
    return jcall(obj, "getB1Bot", jdouble, ())
end

function get_b1_top(obj::NeQuickParameters)
    return jcall(obj, "getB1Top", jdouble, ())
end

function get_b2_bot(obj::NeQuickParameters)
    return jcall(obj, "getB2Bot", jdouble, ())
end

function get_be_bot(obj::NeQuickParameters)
    return jcall(obj, "getBEBot", jdouble, ())
end

function get_be_top(obj::NeQuickParameters)
    return jcall(obj, "getBETop", jdouble, ())
end

function get_h0(obj::NeQuickParameters)
    return jcall(obj, "getH0", jdouble, ())
end

function get_hm_e(obj::NeQuickParameters)
    return jcall(obj, "getHmE", jdouble, ())
end

function get_hm_f1(obj::NeQuickParameters)
    return jcall(obj, "getHmF1", jdouble, ())
end

function get_hm_f2(obj::NeQuickParameters)
    return jcall(obj, "getHmF2", jdouble, ())
end

function get_layer_amplitudes(obj::NeQuickParameters)
    return jcall(obj, "getLayerAmplitudes", Vector{jdouble}, ())
end

function get_nm_f2(obj::NeQuickParameters)
    return jcall(obj, "getNmF2", jdouble, ())
end


text:
@nospecialize
"""The view point of showing patterns
"""
function pretty(points_of_view::Dict{Function, Int})
    viewpoint = points_of_view[pretty]

    function and(ps)
        xs = Any[Print.w("(")]
        for p in ps
            push!(xs, p[viewpoint])
            push!(xs, Print.w(" && "))
        end
        pop!(xs)
        if !isempty(ps)
            push!(xs, Print.w(")"))
        end
        Print.seq(xs...)
    end

    function or(ps)
        xs = Any[Print.w("(")]
        for p in ps
            push!(xs, p[viewpoint])
            push!(xs, Print.w(" || "))
        end
        pop!(xs)
        if !isempty(ps)
            push!(xs, Print.w(")"))
        end
        Print.seq(xs...)
    end
    literal(val) = Print.w(string(val))
    wildcard = Print.w("_")

    function decons(comp::PComp, _, ps)
        Print.seq(Print.w(comp.repr), Print.w("("), getindex.(ps, viewpoint)..., Print.w(")"))
    end

    function guard(pred)
        Print.seq(Print.w("when("), Print.w(repr(pred)), Print.w(")"))
    end

    function effect(eff)
        Print.seq(Print.w("do("), Print.w(repr(eff)), Print.w(")"))
    end

    (
        and = and,
        or = or,
        literal = literal,
        wildcard = wildcard,
        decons = decons,
        guard = guard,
        effect = effect
    )
end
@specialize


text:
s = BenchmarkGroup()
suite["core"] = s

function b_edgeiter(g::AGraphOrDiGraph)
    i = 0
    for e in edges(g)
        i += 1
    end
    return i
end

for G in GLIST
    for n in [100]
        g = CompleteGraph(n, G)
        s["edges","edge iter","fullyconn","$g"] = @benchmarkable b_edgeiter($g)
    end
end

for G in DGLIST
    for n in [100]
        g = CompleteDiGraph(n , G)
        s["edges","edge iter","fullyconn","$g"] = @benchmarkable b_edgeiter($g)
    end
end

for G in GLIST
    for n in [1000]
        g = random_regular_graph(n, 5, G)
        s["edges","edge iter","rrg","$g"] = @benchmarkable b_edgeiter($g)
    end
end

for G in DGLIST
    for n in [1000]
        g = random_regular_digraph(n, 5, G)
        s["edges","edge iter","rrg","$g"] = @benchmarkable b_edgeiter($g)
    end
end

function b_rem_vert(g::AGraphOrDiGraph, vs)
    for v in vs
        rem_vertex!(g, v)
    end
end

Random.seed!(17)
n=100
vs100 = [rand(1:n-i+1) for i=1:n÷4]
n=1000
vs1000 = [rand(1:n-i+1) for i=1:n÷4]
for G in GLIST
    for n in [100]
        g = CompleteGraph(n, G)
        s["vertex","rem vert","fullyconn","$g"] = @benchmarkable b_rem_vert($g, $vs100)
    end
end

for G in DGLIST
    for n in [100]
        g = CompleteDiGraph(n , G)
        s["vertex","rem vert","fullyconn","$g"] = @benchmarkable b_rem_vert($g, $vs100)
    end
end

for G in GLIST
    for n in [1000]
        g = random_regular_graph(n, 5, G)
        s["vertex","rem vert","rrg","$g"] = @benchmarkable b_rem_vert($g, $vs1000)
    end
end

for G in DGLIST
    for n in [1000]
        g = random_regular_digraph(n, 5, G)
        s["vertex","rem vert","rrg","$g"] = @benchmarkable b_rem_vert($g, $vs1000)
    end
end


text:
module TestUpdateOrder
    using Base.Test
    using CategoricalArrays
    using CategoricalArrays: DefaultRefType

    pool = CategoricalPool(
        [
            "a",
            "b",
            "c"
        ],
        Dict(
            "a" => convert(DefaultRefType, 1),
            "b" => convert(DefaultRefType, 2),
            "c" => convert(DefaultRefType, 3),
        )
    )

    order = Vector{DefaultRefType}(length(pool.index))

    CategoricalArrays.buildorder!(order, pool.invindex, ["b", "a", "c"])

    @test order[1] == convert(DefaultRefType, 2)
    @test order[2] == convert(DefaultRefType, 1)
    @test order[3] == convert(DefaultRefType, 3)
end


text:
function decompose(_dataset::AbstractVTKStructuredData, target::String = "Faces", decompose_cell_data = false)
    dataset = VTKStructuredData(_dataset)
    _dim = dim(dataset)

    if decompose_cell_data
        if target == "Faces"
            if _dim == 2
                return decompose_to_faces_2d_with_cell_data(dataset)
            elseif _dim == 3
                return decompose_to_faces_3d_with_cell_data(dataset)
            else
                throw("Unsupported dimension.")
            end
        elseif target == "Lines"
            if _dim == 2
                return decompose_to_lines_2d_with_cell_data(dataset)
            elseif _dim == 3
                return decompose_to_lines_3d_with_cell_data(dataset)
            else
                throw("Unsupported dimension.")
            end
        else
            throw("Unsupported target.")
        end
    else
        if target == "Faces"
            if _dim == 2
                return decompose_to_faces_2d_no_cell_data(dataset)
            elseif _dim == 3
                return decompose_to_faces_3d_no_cell_data(dataset)
            else
                throw("Unsupported dimension.")
            end
        elseif target == "Lines"
            if _dim == 2
                return decompose_to_lines_2d_no_cell_data(dataset)
            elseif _dim == 3
                return decompose_to_lines_3d_no_cell_data(dataset)
            else
                throw("Unsupported dimension.")
            end
        else
            throw("Unsupported target.")
        end
    end
end

function decompose_to_faces_2d_with_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _num_of_cells = num_of_cells(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]
    _cell_types = fill(9, _num_of_cells)
    for cind in Iterators.product(1:cextents[1], 1:cextents[2])
        push!(faces, cell_connectivity(dataset, cind))
    end
    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = reshape(dataset.cell_data[m], (_num_of_cells,))
        else
            _cell_data[m] = reshape(dataset.cell_data[m], (_var_dim, _num_of_cells))
        end
    end
    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function decompose_to_faces_3d_with_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]
    for m in keys(dataset.cell_data)
        _cell_data[m] = Float64[]
    end

    for cind in Iterators.product(1:cextents[1], 1:cextents[2], 1:cextents[3])
        hexa_cc = cell_connectivity(dataset, cind)
        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[3], hexa_cc[4]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1:2]...,cind[3]-1]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][:, cind...] +
                        dataset.cell_data[m][:, cind[1:2]...,cind[3]-1]) / 2
                else
                    _cd = dataset.cell_data[m][:, cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[6], hexa_cc[5]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[2] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1],cind[2]-1,cind[3]]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][:, cind...] +
                        dataset.cell_data[m][:, cind[1],cind[2]-1,cind[3]]) / 2
                else
                    _cd = dataset.cell_data[m][:, cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        push!(faces, [hexa_cc[4], hexa_cc[1], hexa_cc[5], hexa_cc[8]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[1] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1]-1,cind[2:3]...]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][:, cind...] +
                        dataset.cell_data[m][:, cind[1]-1,cind[2:3]...]) / 2
                else
                    _cd = dataset.cell_data[m][:, cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        if cind[1] == cextents[1]
            push!(faces, [hexa_cc[2], hexa_cc[3], hexa_cc[7], hexa_cc[6]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:, cind...]
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[2] == cextents[2]
            push!(faces, [hexa_cc[4], hexa_cc[3], hexa_cc[7], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:, cind...]
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[3] == cextents[3]
            push!(faces, [hexa_cc[5], hexa_cc[6], hexa_cc[7], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:, cind...]
                    append!(_cell_data[m], _cd)
                end
            end
        end
    end
    _cell_types = fill(9, length(faces))

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = reshape(_cell_data[m], (length(faces),))
        else
            _cell_data[m] = reshape(_cell_data[m], (_var_dim, length(faces)))
        end
    end

    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function decompose_to_lines_2d_with_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    lines = Vector{Int}[]
    for m in keys(dataset.cell_data)
        _cell_data[m] = Float64[]
    end
    for cind in Iterators.product(1:cextents[1], 1:cextents[2])
        quad_cc = cell_connectivity(dataset, cind)

        push!(lines, [quad_cc[1], quad_cc[2]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[2] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1], cind[2]-1]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[2] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1],cind[2]-1]) / 2
                else
                    _cd = dataset.cell_data[m][:,cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        push!(lines, [quad_cc[1], quad_cc[4]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[1] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1]-1, cind[2]]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[1] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1]-1,cind[2]]) / 2
                else
                    _cd = dataset.cell_data[m][:,cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        if cind[1] == cextents[1]
            push!(lines, [quad_cc[2], quad_cc[3]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:,cind...]
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[2] == cextents[2]
            push!(lines, [quad_cc[3], quad_cc[4]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:,cind...]
                    append!(_cell_data[m], _cd)
                end
            end
        end
        _cell_types = fill(3, length(lines))
    end

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = reshape(dataset.cell_data[m], (length(lines),))
        else
            _cell_data[m] = reshape(dataset.cell_data[m], (_var_dim, length(lines)))
        end
    end

    return VTKPolyData(point_coords, _cell_types, lines, point_data, _cell_data)
end

function decompose_to_lines_3d_with_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    lines = Vector{Int}[]
    for m in keys(dataset.cell_data)
        _cell_data[m] = Float64[]
    end

    for cind in Iterators.product(1:cextents[1], 1:cextents[2], 1:cextents[3])
        hexa_cc = cell_connectivity(dataset, cind)
        push!(lines, [hexa_cc[1], hexa_cc[2]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[2] != 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1:2]...,cind[3]-1] +
                        dataset.cell_data[m][cind[1],cind[2]-1,cind[3]] +
                        dataset.cell_data[m][cind[1],cind[2]-1,cind[3]-1]) / 4
                elseif cind[2] == 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1:2]...,cind[3]-1]) / 2
                elseif cind[2] != 1 && cind[3] == 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1], cind[2]-1,cind[3]]) / 2
                elseif cind[2] == cind[3] == 1
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[2] != 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1:2]...,cind[3]-1] +
                        dataset.cell_data[m][:,cind[1],cind[2]-1,cind[3]] +
                        dataset.cell_data[m][:,cind[1],cind[2]-1,cind[3]-1]) / 4
                elseif cind[2] == 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1:2]...,cind[3]-1]) / 2
                elseif cind[2] != 1 && cind[3] == 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1], cind[2]-1,cind[3]]) / 2
                elseif cind[2] == cind[3] == 1
                    _cd = dataset.cell_data[m][:,cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        push!(lines, [hexa_cc[1], hexa_cc[4]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[1] != 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1:2]...,cind[3]-1] +
                        dataset.cell_data[m][cind[1]-1,cind[2],cind[3]] +
                        dataset.cell_data[m][cind[1]-1,cind[2],cind[3]-1]) / 4
                elseif cind[1] == 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1:2]...,cind[3]-1]) / 2
                elseif cind[1] != 1 && cind[3] == 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1]-1, cind[2],cind[3]]) / 2
                elseif cind[1] == cind[3] == 1
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[1] != 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1:2]...,cind[3]-1] +
                        dataset.cell_data[m][:,cind[1]-1,cind[2],cind[3]] +
                        dataset.cell_data[m][:,cind[1]-1,cind[2],cind[3]-1]) / 4
                elseif cind[1] == 1 && cind[3] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1:2]...,cind[3]-1]) / 2
                elseif cind[1] != 1 && cind[3] == 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1]-1, cind[2],cind[3]]) / 2
                elseif cind[1] == cind[3] == 1
                    _cd = dataset.cell_data[m][:,cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        push!(lines, [hexa_cc[1], hexa_cc[5]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[1] != 1 && cind[2] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1]-1,cind[2:3]...] +
                        dataset.cell_data[m][cind[1],cind[2]-1,cind[3]] +
                        dataset.cell_data[m][cind[1]-1,cind[2]-1,cind[3]]) / 4
                elseif cind[1] == 1 && cind[2] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1],cind[2]-1,cind[3]]) / 2
                elseif cind[1] != 1 && cind[2] == 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1]-1,cind[2],cind[3]]) / 2
                elseif cind[1] == cind[2] == 1
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
            else
                if cind[1] != 1 && cind[2] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1]-1,cind[2:3]...] +
                        dataset.cell_data[m][:,cind[1],cind[2]-1,cind[3]] +
                        dataset.cell_data[m][:,cind[1]-1,cind[2]-1,cind[3]]) / 4
                elseif cind[1] == 1 && cind[2] != 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1],cind[2]-1,cind[3]]) / 2
                elseif cind[1] != 1 && cind[2] == 1
                    _cd = (dataset.cell_data[m][:,cind...] +
                        dataset.cell_data[m][:,cind[1]-1,cind[2],cind[3]]) / 2
                elseif cind[1] == cind[2] == 1
                    _cd = dataset.cell_data[m][:,cind...]
                end
                append!(_cell_data[m], _cd)
            end
        end

        if cind[1] == cextents[1]
            push!(lines, [hexa_cc[2], hexa_cc[3]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[3] != 1
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1:2]..., cind[3]-1]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[3] != 1
                        _cd = (dataset.cell_data[m][:,cind...] +
                            dataset.cell_data[m][:,cind[1:2]...,cind[3]-1]) / 2
                    else
                        _cd = dataset.cell_data[m][:,cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end

            push!(lines, [hexa_cc[3], hexa_cc[7]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[2] != cextents[2]
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1], cind[2]+1, cind[3]]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[2] != cextents[2]
                        _cd = (dataset.cell_data[m][:,cind...] +
                            dataset.cell_data[m][:,cind[1],cind[2]+1,cind[3]]) / 2
                    else
                        _cd = dataset.cell_data[m][:,cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end

            push!(lines, [hexa_cc[6], hexa_cc[7]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[3] != cextents[3]
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1:2]..., cind[3]+1]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[3] != cextents[3]
                        _cd = (dataset.cell_data[m][:,cind...] +
                            dataset.cell_data[m][:,cind[1:2]...,cind[3]+1]) / 2
                    else
                        _cd = dataset.cell_data[m][:,cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end

            push!(lines, [hexa_cc[2], hexa_cc[6]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[2] != 1
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1], cind[2]-1, cind[3]]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[2] != 1
                        _cd = (dataset.cell_data[m][:, cind...] +
                            dataset.cell_data[m][:, cind[1], cind[2]-1, cind[3]]) / 2
                    else
                        _cd = dataset.cell_data[m][:, cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[2] == cextents[2]
            push!(lines, [hexa_cc[3], hexa_cc[4]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[3] != 1
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1:2]..., cind[3]-1]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[3] != 1
                        _cd = (dataset.cell_data[m][:, cind...] +
                            dataset.cell_data[m][:, cind[1:2]..., cind[3]-1]) / 2
                    else
                        _cd = dataset.cell_data[m][:, cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end

            push!(lines, [hexa_cc[7], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[3] != cextents[3]
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1:2]..., cind[3]+1]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[3] != cextents[3]
                        _cd = (dataset.cell_data[m][:, cind...] +
                            dataset.cell_data[m][:, cind[1:2]..., cind[3]+1]) / 2
                    else
                        _cd = dataset.cell_data[m][:, cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end

            push!(lines, [hexa_cc[4], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[1] != 1
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1]-1, cind[2:3]...]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[1] != 1
                        _cd = (dataset.cell_data[m][:, cind...] +
                            dataset.cell_data[m][:, cind[1]-1, cind[2:3]...]) / 2
                    else
                        _cd = dataset.cell_data[m][:, cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[3] == cextents[3]
            push!(lines, [hexa_cc[5], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[1] != 1
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1]-1, cind[2:3]...]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[1] != 1
                        _cd = (dataset.cell_data[m][:, cind...] +
                            dataset.cell_data[m][:, cind[1]-1, cind[2:3]...]) / 2
                    else
                        _cd = dataset.cell_data[m][:, cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end

            push!(lines, [hexa_cc[5], hexa_cc[6]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    if cind[2] != 1
                        _cd = (dataset.cell_data[m][cind...] +
                            dataset.cell_data[m][cind[1], cind[2]-1, cind[3]]) / 2
                    else
                        _cd = dataset.cell_data[m][cind...]
                    end
                    push!(_cell_data[m], _cd)
                else
                    if cind[2] != 1
                        _cd = (dataset.cell_data[m][:, cind...] +
                            dataset.cell_data[m][:, cind[1], cind[2]-1, cind[3]]) / 2
                    else
                        _cd = dataset.cell_data[m][:, cind...]
                    end
                    append!(_cell_data[m], _cd)
                end
            end
        end
    end
    _cell_types = fill(3, length(lines))

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = reshape(dataset.cell_data[m], (length(lines),))
        else
            _cell_data[m] = reshape(dataset.cell_data[m], (_var_dim, length(lines)))
        end
    end

    return VTKPolyData(point_coords, _cell_types, lines, point_data, _cell_data)
end

function decompose_to_faces_2d_no_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _num_of_cells = num_of_cells(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]
    _cell_types = fill(9, _num_of_cells)
    for cind in Iterators.product(1:cextents[1], 1:cextents[2])
        push!(faces, cell_connectivity(dataset, cind))
    end
    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function decompose_to_faces_3d_no_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]

    for cind in Iterators.product(1:cextents[1], 1:cextents[2], 1:cextents[3])
        hexa_cc = cell_connectivity(dataset, cind)
        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[3], hexa_cc[4]])
        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[6], hexa_cc[5]])
        push!(faces, [hexa_cc[4], hexa_cc[1], hexa_cc[5], hexa_cc[8]])
        if cind[1] == cextents[1]
            push!(faces, [hexa_cc[2], hexa_cc[3], hexa_cc[7], hexa_cc[6]])
        end
        if cind[2] == cextents[2]
            push!(faces, [hexa_cc[4], hexa_cc[3], hexa_cc[7], hexa_cc[8]])
        end
        if cind[3] == cextents[3]
            push!(faces, [hexa_cc[5], hexa_cc[6], hexa_cc[7], hexa_cc[8]])
        end
    end
    _cell_types = fill(9, length(faces))

    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function decompose_to_lines_2d_no_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    lines = Vector{Int}[]

    for cind in Iterators.product(1:cextents[1], 1:cextents[2])
        quad_cc = cell_connectivity(dataset, cind)

        push!(lines, [quad_cc[1], quad_cc[2]])
        push!(lines, [quad_cc[1], quad_cc[4]])
        if cind[1] == cextents[1]
            push!(lines, [quad_cc[2], quad_cc[3]])
        end
        if cind[2] == cextents[2]
            push!(lines, [quad_cc[3], quad_cc[4]])
        end
        _cell_types = fill(3, length(lines))
    end

    return VTKPolyData(point_coords, _cell_types, lines, point_data, _cell_data)
end

function decompose_to_lines_3d_no_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    lines = Vector{Int}[]

    for cind in Iterators.product(1:cextents[1], 1:cextents[2], 1:cextents[3])
        hexa_cc = cell_connectivity(dataset, cind)
        push!(lines, [hexa_cc[1], hexa_cc[2]])
        push!(lines, [hexa_cc[1], hexa_cc[4]])
        push!(lines, [hexa_cc[1], hexa_cc[5]])
        if cind[1] == cextents[1]
            push!(lines, [hexa_cc[2], hexa_cc[3]])
            push!(lines, [hexa_cc[3], hexa_cc[7]])
            push!(lines, [hexa_cc[6], hexa_cc[7]])
            push!(lines, [hexa_cc[2], hexa_cc[6]])
        end
        if cind[2] == cextents[2]
            push!(lines, [hexa_cc[3], hexa_cc[4]])
            push!(lines, [hexa_cc[7], hexa_cc[8]])
            push!(lines, [hexa_cc[4], hexa_cc[8]])
        end
        if cind[3] == cextents[3]
            push!(lines, [hexa_cc[5], hexa_cc[8]])
            push!(lines, [hexa_cc[5], hexa_cc[6]])
        end
    end
    _cell_types = fill(3, length(lines))

    return VTKPolyData(point_coords, _cell_types, lines, point_data, _cell_data)
end


text:
#####
##### Advective tracer fluxes
#####

# Calculate the flux of a tracer quantity c through the faces of a cell.
# In this case, the fluxes are given by u*Ax*T̅ˣ, v*Ay*T̅ʸ, and w*Az*T̅ᶻ.
@inline advective_tracer_flux_x(i, j, k, grid, u, c) = Ax_ψᵃᵃᶠ(i, j, k, grid, u) * ℑxᶠᵃᵃ(i, j, k, grid, c)
@inline advective_tracer_flux_y(i, j, k, grid, v, c) = Ay_ψᵃᵃᶠ(i, j, k, grid, v) * ℑyᵃᶠᵃ(i, j, k, grid, c)
@inline advective_tracer_flux_z(i, j, k, grid, w, c) = Az_ψᵃᵃᵃ(i, j, k, grid, w) * ℑzᵃᵃᶠ(i, j, k, grid, c)

#####
##### Tracer advection operator
#####

"""
    div_flux(i, j, k, grid, U, c)

Calculates the divergence of the flux of a tracer quantity c being advected by
a velocity field U = (u, v, w), ∇·(Uc),

    1/V * [δxᶜᵃᵃ(Ax * u * ℑxᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * v * ℑyᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * w * ℑzᵃᵃᶠ(c))]

which will end up at the location `ccc`.
"""
@inline function div_uc(i, j, k, grid, U, c)
    1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, advective_tracer_flux_x, U.u, c) +
                             δyᵃᶜᵃ(i, j, k, grid, advective_tracer_flux_y, U.v, c) +
                             δzᵃᵃᶜ(i, j, k, grid, advective_tracer_flux_z, U.w, c))
end


text:
#---------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Model of an yield reactor
#*----------------------------------------------------------------------
#*
#*   Description:
#*       Modeling of a reactor based on an yield approach.
#*
#*   Assumptions:
#*		* single- and two-phases involved
#*       * steady-state
#*
#*	Specify:
#*		* inlet stream
#*		* component yield or
#*		* reaction yield
#*
#*----------------------------------------------------------------------
#* Author: Rodolfo Rodrigues
#* $Id$
#*--------------------------------------------------------------------
#---------------------------------------------------------------------
#*	only vapour phase
#*--------------------------------------------------------------------
type yield_vap
	yield_vap()=begin
		new(
			tank_vap(),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of reactions",
				:Default=>1
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Key component",
				:Lower=>1,
				:Default=>1
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Overall component rate of reaction"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Fractional conversion of component",
				:Symbol=>"X",
				:Default=>0
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Molar component yield (global)",
				:Symbol=>"Y_G"
			)),(NComp)),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Molar reaction yield (instantaneous)",
				:Symbol=>"Y_I"
			)),(NComp)),
			[
				:(Outlet.F*Outlet.z = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z + rate*_base_1._base_1.Tank.V),
				:(rate*_base_1._base_1.Tank.V = _base_1._base_1.Outletm.F*(yield/(1 + yield(KComp))*_base_1._base_1.Outletm.z(KComp) - _base_1._base_1.Outletm.z)),
				:(rate = yield_*rate(KComp)),
				:(Outlet.P = _base_1._base_1.Outletm.P),
				:(Outlet.F*Outlet.h = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.h),
				:(Outlet.F*Outlet.z([1:_base_1._base_1.NComp]) = _base_1._base_1.Outletm.F*_base_1._base_1.Outletm.z([1:_base_1._base_1.NComp])*(1 - conv([1:_base_1._base_1.NComp]))),
				:(conv([1:_base_1._base_1.NComp]) = 1),
				:(conv([1:_base_1._base_1.NComp]) = 0),
			],
			[
				"Outlet stream","Rate of reaction","Instantaneous yield","Mechanical equilibrium","Energy balance","Molar conversion","Molar conversion","Molar conversion",
			],
			[:NReac,:KComp,],
			[:Outlet,:rate,:conv,:yield,:yield_,]
		)
	end
	_base_1::tank_vap
	NReac::DanaInteger
	KComp::DanaInteger
	Outlet::vapour_stream
	rate::Array{reaction_mol}
	conv::Array{DanaReal}
	yield::Array{DanaReal}
	yield_::Array{DanaReal}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export yield_vap
function setEquationFlow(in::yield_vap)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	if (Outletm.z(i) > 1e-16)
		addEquation(6)
	else
		if (Outlet.z(i) > 0)
			addEquation(7)
			# ?

		else
			addEquation(8)
			# ?

		end
	end
end
function atributes(in::yield_vap,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/cstr"
	fields[:Brief]="Model of a generic vapour-phase yield CSTR"
	fields[:Info]="
== Assumptions ==
* only vapour-phase
* steady-state

== Specify ==
* inlet stream
* component yield or
* reaction yield
"
	drive!(fields,_)
	return fields
end
yield_vap(_::Dict{Symbol,Any})=begin
	newModel=yield_vap()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
@doc doc"""
Approximation de la solution du sous-problème ``q_k(s) = s^{t}g + (1/2)s^{t}Hs``
        avec ``s=-t g_k,t > 0,||s||< \delta_k ``

# Syntaxe
```julia
s1, e1 = Pas_De_Cauchy(gradient,Hessienne,delta)
```

# Entrées
 * **gradfk** : (Array{Float,1}) le gradient de la fonction f appliqué au point ``x_k``
 * **hessfk** : (Array{Float,2}) la Hessienne de la fonction f appliqué au point ``x_k``
 * **delta**  : (Float) le rayon de la région de confiance

# Sorties
 * **s** : (Array{Float,1}) une approximation de la  solution du sous-problème
 * **e** : (Integer) indice indiquant l'état de sortie:
        si g != 0
            si on ne sature pas la boule
              e <- 1
            sinon
              e <- -1
        sinon
            e <- 0

# Exemple d'appel
```julia
g1 = [0; 0]
H1 = [7 0 ; 0 2]
delta1 = 1
s1, e1 = Pas_De_Cauchy(g1,H1,delta1)
```
"""
function Pas_De_Cauchy(g,H,delta)

  a = (1/2)*g'*H*g
  b = norm(g)

  if (b > 0)
    t = (b^2)/(2*a)
    s = -t*g
    e = 1
    if norm(s)>delta
        t = delta/b
        s = -t*g
        e = -1
    end
  else
    t = 0
    s = -t*g
    e = 0
  end

  return s, e
end

max_stars_repo_path: 2A/S7/UE Optimisation et RO/Optimisation/Optinum/src/Pas_De_Cauchy.jl

max_issues_repo_name: Sajid-Badr/ENSEEIHT

max_forks_repo_name: Sajid-Badr/ENSEEIHT


text:
@symbol_func function K_BS(cur_reactor::AbstractReactor)
  cur_K = 4.879

  cur_K *= K_n(cur_reactor)

  cur_K *= 1 + cur_reactor.kappa_95 ^ 2

  cur_K /= 2

  cur_K *= sqrt(cur_reactor.epsilon) ^ 5

  cur_K *= _C_B(cur_reactor)

  cur_K
end

@symbol_func function _C_B(cur_reactor::AbstractReactor)
  cur_nu_n = cur_reactor.nu_n
  cur_nu_T = cur_reactor.nu_T

  cur_func = function (cur_rho)
    cur_value = 1.0 - cur_rho ^ 2

    cur_value ^= cur_nu_n + cur_nu_T - 1

    cur_value *= sqrt(cur_rho) ^ 5

    cur_value /= b_p(cur_reactor.gamma, cur_rho)

    cur_value
  end

  cur_reactor.is_symbolic &&
    ( cur_func = cur_func(rho_sym) )

  cur_C_B = norm_int(cur_func)

  cur_C_B *= 1 + cur_nu_n

  cur_C_B *= 1 + cur_nu_T

  cur_C_B *= cur_nu_n + 0.054 * cur_nu_T

  cur_C_B
end


text:
# Command: Take Screenshot
"""
    screenshot(session::Session)::String
    screenshot(element::Element)::String
    screenshot(source::Union{Session, Element}, sink::Union{AbstractString, IO})

Take Screenshot and optionally saves it.
"""
function screenshot end

function screenshot(session::Session)::String
    @unpack addr, id = session
    response =
        HTTP.get("$addr/session/$id/screenshot", [("Content-Type" => "application/json; charset=utf-8")])
    @assert response.status == 200
    JSON3.read(response.body).value
end
function screenshot(element::Element)::String
    @unpack addr, id = element.session
    response = HTTP.get(
        "$addr/session/$id/element/$(element.id)/screenshot",
        [("Content-Type" => "application/json; charset=utf-8")],
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end
screenshot(source::Union{Session, Element}, sink::Union{AbstractString, IO}) =
    write(sink, base64decode(source))


text:
using MAT

#############
# inference #
#############

function get_data(data_name::String, field::String="data_sq")
    data_file = matopen(data_name)
    data = read(data_file, field)
    close(data_file)

    return data
end

function sample(data::Matrix, n::Integer)
    data_indices = sort!(rand(1:size(data, 1), n))

    return data[data_indices, 1] # 1: x; 2: θ
end

reshape_inferred_data(data::Matrix) = [data[:, i] for i in 1:size(data, 2)]

function get_𝛒_and_𝐰(
    data_name::String;
    n_sample=10, fix_θ=true,
    wf=WignerFunction(LinRange(-3, 3, 101), LinRange(-3, 3, 101), dim=70),
    m=get_model("model")
)
    data = get_data(data_name)

    argv = zeros(Float32, 6)
    for _ in 1:n_sample
        argv += SqState.infer(m, sample(data, 4096))
    end
    argv ./= n_sample

    argv[4:6] ./= sum(argv[4:6])
    r, θ, n̄, c1, c2, c3 = argv
    θ = fix_θ ? zero(typeof(θ)) : θ

    state = SqState.construct_state(r, θ, n̄, c1, c2, c3, wf.m_dim)

    return reshape_inferred_data(real.(𝛒(state))), reshape_inferred_data(wf(state).𝐰_surface)
end

#############
# gen plots #
#############

const color_config = [:colorscale=>"twilight", :cmid=>0, :zmid=>0]

function wigner_surface(data::Vector, width::Integer, height::Integer)
    data = PlotlyJS.surface(z=data; color_config...)
    layout = Layout(title="Wigner Function", width=width, height=height)

    return plot(data, layout)
end

function wigner_heatmap(data::Vector, width::Integer, height::Integer)
    data = PlotlyJS.heatmap(z=data; color_config...)
    layout = Layout(title="Wigner Function", width=width, height=height)

    return plot(data, layout)
end

function density_matrix(data::Vector, width::Integer, height::Integer)
    data = PlotlyJS.heatmap(z=data; color_config...)
    layout = Layout(title="Density Matrix (Real Part)", width=width, height=height)

    return plot(data, layout)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# OffsetArrays (arrays with indexing that doesn't start at 1)

# This test file is designed to exercise support for generic indexing,
# even though offset arrays aren't implemented in Base.

# OffsetArrays v1.3.0
# No compat patch and docstrings
module OffsetArrays

using Base: tail, @propagate_inbounds
using Base: IdentityUnitRange

export OffsetArray, OffsetMatrix, OffsetVector

struct IdOffsetRange{T<:Integer,I<:AbstractUnitRange{T}} <: AbstractUnitRange{T}
    parent::I
    offset::T

    IdOffsetRange{T,I}(r::I, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}} = new{T,I}(r, offset)
end

# Construction/coercion from arbitrary AbstractUnitRanges
function IdOffsetRange{T,I}(r::AbstractUnitRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, o = offset_coerce(I, r)
    return IdOffsetRange{T,I}(rc, convert(T, o+offset))
end
function IdOffsetRange{T}(r::AbstractUnitRange, offset::Integer = 0) where T<:Integer
    rc = convert(AbstractUnitRange{T}, r)::AbstractUnitRange{T}
    return IdOffsetRange{T,typeof(rc)}(rc, convert(T, offset))
end
IdOffsetRange(r::AbstractUnitRange{T}, offset::Integer = 0) where T<:Integer =
    IdOffsetRange{T,typeof(r)}(r, convert(T, offset))

# Coercion from other IdOffsetRanges
IdOffsetRange{T,I}(r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
function IdOffsetRange{T,I}(r::IdOffsetRange) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, offset = offset_coerce(I, r.parent)
    return IdOffsetRange{T,I}(rc, r.offset+offset)
end
function IdOffsetRange{T}(r::IdOffsetRange) where T<:Integer
    return IdOffsetRange(convert(AbstractUnitRange{T}, r.parent), r.offset)
end
IdOffsetRange(r::IdOffsetRange) = r

AbstractUnitRange{T}(r::IdOffsetRange{T}) where {T} = r
AbstractUnitRange{T}(r::IdOffsetRange) where {T} = IdOffsetRange{T}(r)

# TODO: uncomment these when Julia is ready
# # Conversion preserves both the values and the indexes, throwing an InexactError if this
# # is not possible.
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r.parent), r.offset)
# Base.convert(::Type{IdOffsetRange{T,I}}, r::AbstractUnitRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r), 0)

offset_coerce(::Type{Base.OneTo{T}}, r::Base.OneTo) where T<:Integer = convert(Base.OneTo{T}, r), 0
function offset_coerce(::Type{Base.OneTo{T}}, r::AbstractUnitRange) where T<:Integer
    o = first(r) - 1
    return Base.OneTo{T}(last(r) - o), o
end
# function offset_coerce(::Type{Base.OneTo{T}}, r::IdOffsetRange) where T<:Integer
#     rc, o = offset_coerce(Base.OneTo{T}, r.parent)

# Fallback, specialze this method if `convert(I, r)` doesn't do what you need
offset_coerce(::Type{I}, r::AbstractUnitRange) where I<:AbstractUnitRange{T} where T =
    convert(I, r), 0

@inline Base.parent(r::IdOffsetRange) = r.parent
@inline Base.axes(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.axes1(r::IdOffsetRange) = IdOffsetRange(Base.axes1(r.parent), r.offset)
@inline Base.length(r::IdOffsetRange) = length(r.parent)
Base.reduced_index(i::IdOffsetRange) = typeof(i)(first(i):first(i))
# Workaround for #92 on Julia < 1.4
Base.reduced_index(i::IdentityUnitRange{<:IdOffsetRange}) = typeof(i)(first(i):first(i))
for f in [:firstindex, :lastindex]
    @eval Base.$f(r::IdOffsetRange) = $f(r.parent) .+ r.offset
end

@inline function Base.iterate(r::IdOffsetRange)
    ret = iterate(r.parent)
    ret === nothing && return nothing
    return (ret[1] + r.offset, ret[2])
end
@inline function Base.iterate(r::IdOffsetRange, i)
    ret = iterate(r.parent, i)
    ret === nothing && return nothing
    return (ret[1] + r.offset, ret[2])
end

@inline Base.first(r::IdOffsetRange) = first(r.parent) + r.offset
@inline Base.last(r::IdOffsetRange) = last(r.parent) + r.offset

@propagate_inbounds Base.getindex(r::IdOffsetRange, i::Integer) = r.parent[i - r.offset] + r.offset
@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::AbstractUnitRange{<:Integer})
    return r.parent[s .- r.offset] .+ r.offset
end
@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::IdentityUnitRange)
    return IdOffsetRange(r.parent[s .- r.offset], r.offset)
end
@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::IdOffsetRange)
    return IdOffsetRange(r.parent[s.parent .+ (s.offset - r.offset)] .+ (r.offset - s.offset), s.offset)
end

# offset-preserve broadcasting
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(-), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .- x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .+ x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), x::Integer, r::IdOffsetRange{T}) where T =
    IdOffsetRange{T}(x .+ r.parent, r.offset)

Base.show(io::IO, r::IdOffsetRange) = print(io, "OffsetArrays.IdOffsetRange(", first(r), ':', last(r), ")")

# Optimizations
@inline Base.checkindex(::Type{Bool}, inds::IdOffsetRange, i::Real) = Base.checkindex(Bool, inds.parent, i - inds.offset)

struct Origin{T <: Union{Tuple,Int}}
    index::T
end
Origin(I::NTuple{N,Int}) where N = Origin{typeof(I)}(I)
Origin(I::CartesianIndex) = Origin(I.I)
Origin(I1::Int, In::Int...) = Origin((I1, In...))
# Origin(0) != Origin((0, )) but they work the same with broadcasting
Origin(n::Int) = Origin{Int}(n)

(o::Origin)(A::AbstractArray) = o.index .- first.(axes(A))

### Low-level utilities ###

_indexoffset(r::AbstractRange) = first(r) - 1
_indexoffset(i::Integer) = 0
_indexoffset(i::Colon) = 0
_indexlength(r::AbstractRange) = length(r)
_indexlength(i::Integer) = i
_indexlength(i::Colon) = Colon()

_offset(axparent::AbstractUnitRange, ax::AbstractUnitRange) = first(ax) - first(axparent)
_offset(axparent::AbstractUnitRange, ax::Integer) = 1 - first(axparent)

abstract type AxisConversionStyle end
struct SingleRange <: AxisConversionStyle end
struct TupleOfRanges <: AxisConversionStyle end

AxisConversionStyle(::Type) = SingleRange()
AxisConversionStyle(::Type{<:CartesianIndices}) = TupleOfRanges()

_convertTupleAbstractUnitRange(x) = _convertTupleAbstractUnitRange(AxisConversionStyle(typeof(x)), x)
_convertTupleAbstractUnitRange(::SingleRange, x) = (convert(AbstractUnitRange{Int}, x),)
_convertTupleAbstractUnitRange(::TupleOfRanges, x) = convert(Tuple{Vararg{AbstractUnitRange{Int}}}, x)

_toAbstractUnitRanges(t::Tuple) = (_convertTupleAbstractUnitRange(first(t))..., _toAbstractUnitRanges(tail(t))...)
_toAbstractUnitRanges(::Tuple{}) = ()

# ensure that the indices are consistent in the constructor
_checkindices(A::AbstractArray, indices, label) = _checkindices(ndims(A), indices, label)
function _checkindices(N::Integer, indices, label)
    throw_argumenterror(N, indices, label) = throw(ArgumentError(label * " $indices are not compatible with a $(N)D array"))
    N == length(indices) || throw_argumenterror(N, indices, label)
end

# Technically we know the length of CartesianIndices but we need to convert it first, so here we
# don't put it in OffsetAxisKnownLength.
const OffsetAxisKnownLength = Union{Integer,AbstractUnitRange}
const OffsetAxis = Union{OffsetAxisKnownLength,Colon}
const ArrayInitializer = Union{UndefInitializer,Missing,Nothing}

## OffsetArray
struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
    function OffsetArray{T,N,AA}(parent::AA, offsets::NTuple{N,Int}) where {T,N,AA <: AbstractArray}
        @boundscheck overflow_check.(axes(parent), offsets)
        new{T,N,AA}(parent, offsets)
    end
end

const OffsetVector{T,AA <: AbstractArray} = OffsetArray{T,1,AA}

const OffsetMatrix{T,AA <: AbstractArray} = OffsetArray{T,2,AA}

function overflow_check(r, offset::T) where T
    # This gives some performance boost https://github.com/JuliaLang/julia/issues/33273
    throw_upper_overflow_error() = throw(ArgumentError("Boundary overflow detected: offset $offset should be equal or less than $(typemax(T) - last(r))"))
    throw_lower_overflow_error() = throw(ArgumentError("Boundary overflow detected: offset $offset should be equal or greater than $(typemin(T) - first(r))"))

    if offset > 0 && last(r) > typemax(T) - offset
        throw_upper_overflow_error()
    elseif offset < 0 && first(r) < typemin(T) - offset
        throw_lower_overflow_error()
    end
end

# Tuples of integers are treated as offsets
# Empty Tuples are handled here
function OffsetArray(A::AbstractArray, offsets::Tuple{Vararg{Integer}})
    _checkindices(A, offsets, "offsets")
    OffsetArray{eltype(A),ndims(A),typeof(A)}(A, offsets)
end

# These methods are necessary to disallow incompatible dimensions for
# the OffsetVector and the OffsetMatrix constructors
for (FT, ND) in ((:OffsetVector, :1), (:OffsetMatrix, :2))
    @eval function $FT(A::AbstractArray{<:Any,$ND}, offsets::Tuple{Vararg{Integer}})
        _checkindices(A, offsets, "offsets")
        OffsetArray{eltype(A),$ND,typeof(A)}(A, offsets)
    end
    FTstr = string(FT)
    @eval function $FT(A::AbstractArray, offsets::Tuple{Vararg{Integer}})
        throw(ArgumentError($FTstr * " requires a " * string($ND) * "D array"))
    end
end

## OffsetArray constructors
for FT in (:OffsetArray, :OffsetVector, :OffsetMatrix)
    # Nested OffsetArrays may strip off the wrapper and collate the offsets
    @eval function $FT(A::OffsetArray, offsets::Tuple{Vararg{Integer}})
        _checkindices(A, offsets, "offsets")
        $FT(parent(A), map(+, A.offsets, offsets))
    end

    # In general, indices get converted to AbstractUnitRanges.
    # CartesianIndices{N} get converted to N ranges
    @eval function $FT(A::AbstractArray, inds::Tuple{Any,Vararg{Any}})
        $FT(A, _toAbstractUnitRanges(to_indices(A, axes(A), inds)))
    end

    # convert ranges to offsets
    @eval function $FT(A::AbstractArray, inds::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}})
        _checkindices(A, inds, "indices")
        # Performance gain by wrapping the error in a function: see https://github.com/JuliaLang/julia/issues/37558
        throw_dimerr(lA, lI) = throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
        lA = size(A)
        lI = map(length, inds)
        lA == lI || throw_dimerr(lA, lI)
        $FT(A, map(_offset, axes(A), inds))
    end

    @eval $FT(A::AbstractArray, inds::Vararg) = $FT(A, inds)

    @eval $FT(A::AbstractArray, origin::Origin) = $FT(A, origin(A))
end

# array initialization
function OffsetArray{T,N}(init::ArrayInitializer, inds::Tuple{Vararg{OffsetAxisKnownLength}}) where {T,N}
    _checkindices(N, inds, "indices")
    AA = Array{T,N}(init, map(_indexlength, inds))
    OffsetArray{T,N,typeof(AA)}(AA, map(_indexoffset, inds))
end
function OffsetArray{T,N}(init::ArrayInitializer, inds::Tuple) where {T,N}
    OffsetArray{T,N}(init, _toAbstractUnitRanges(inds))
end
OffsetArray{T,N}(init::ArrayInitializer, inds::Vararg) where {T,N} = OffsetArray{T,N}(init, inds)

OffsetArray{T}(init::ArrayInitializer, inds::NTuple{N,OffsetAxisKnownLength}) where {T,N} = OffsetArray{T,N}(init, inds)
function OffsetArray{T}(init::ArrayInitializer, inds::Tuple) where {T}
    OffsetArray{T}(init, _toAbstractUnitRanges(inds))
end
OffsetArray{T}(init::ArrayInitializer, inds::Vararg) where {T} = OffsetArray{T}(init, inds)

Base.IndexStyle(::Type{OA}) where {OA <: OffsetArray} = IndexStyle(parenttype(OA))
parenttype(::Type{OffsetArray{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::OffsetArray) = parenttype(typeof(A))

Base.parent(A::OffsetArray) = A.parent

Base.eachindex(::IndexCartesian, A::OffsetArray) = CartesianIndices(axes(A))
Base.eachindex(::IndexLinear, A::OffsetVector)   = axes(A, 1)

@inline Base.size(A::OffsetArray) = size(parent(A))
@inline Base.size(A::OffsetArray, d) = size(parent(A), d)

@inline Base.axes(A::OffsetArray) = map(IdOffsetRange, axes(parent(A)), A.offsets)
@inline Base.axes(A::OffsetArray, d) = d <= ndims(A) ? IdOffsetRange(axes(parent(A), d), A.offsets[d]) : IdOffsetRange(axes(parent(A), d))
@inline Base.axes1(A::OffsetArray{T,0}) where {T} = IdOffsetRange(axes(parent(A), 1))  # we only need to specialize this one

Base.similar(A::OffsetArray, ::Type{T}, dims::Dims) where T =
    similar(parent(A), T, dims)
function Base.similar(A::AbstractArray, ::Type{T}, inds::Tuple{OffsetAxisKnownLength,Vararg{OffsetAxisKnownLength}}) where T
    B = similar(A, T, map(_indexlength, inds))
    return OffsetArray(B, map(_offset, axes(B), inds))
end

# reshape accepts a single colon
Base.reshape(A::AbstractArray, inds::OffsetAxis...) = reshape(A, inds)
function Base.reshape(A::AbstractArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}})
    AR = reshape(A, map(_indexlength, inds))
    return OffsetArray(AR, map(_offset, axes(AR), inds))
end

# Reshaping OffsetArrays can "pop" the original OffsetArray wrapper and return
# an OffsetArray(reshape(...)) instead of an OffsetArray(reshape(OffsetArray(...)))
Base.reshape(A::OffsetArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) =
    OffsetArray(reshape(parent(A), map(_indexlength, inds)), map(_indexoffset, inds))
# And for non-offset axes, we can just return a reshape of the parent directly
Base.reshape(A::OffsetArray, inds::Tuple{Union{Integer,Base.OneTo},Vararg{Union{Integer,Base.OneTo}}}) = reshape(parent(A), inds)
Base.reshape(A::OffsetArray, inds::Dims) = reshape(parent(A), inds)
Base.reshape(A::OffsetArray, ::Colon) = reshape(parent(A), Colon())
Base.reshape(A::OffsetVector, ::Colon) = A
Base.reshape(A::OffsetVector, ::Tuple{Colon}) = A
Base.reshape(A::OffsetArray, inds::Union{Int,Colon}...) = reshape(parent(A), inds)
Base.reshape(A::OffsetArray, inds::Tuple{Vararg{Union{Int,Colon}}}) = reshape(parent(A), inds)

function Base.similar(::Type{T}, shape::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where {T <: AbstractArray}
    P = T(undef, map(_indexlength, shape))
    OffsetArray(P, map(_offset, axes(P), shape))
end

Base.fill(v, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(Array{typeof(v)}, inds), v)
Base.zeros(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(similar(Array{T}, inds), zero(T))
Base.ones(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(similar(Array{T}, inds), one(T))
Base.trues(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(BitArray, inds), true)
Base.falses(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(BitArray, inds), false)

## Indexing

# Note this gets the index of the parent *array*, not the index of the parent *range*
# Here's how one can think about this:
#   Δi = i - first(r)
#   i′ = first(r.parent) + Δi
# and one obtains the result below.
parentindex(r::IdOffsetRange, i) = i - r.offset

@inline function Base.getindex(A::OffsetArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    J = map(parentindex, axes(A), I)
    @inbounds parent(A)[J...]
end

@inline function Base.getindex(A::OffsetVector, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[parentindex(Base.axes1(A), i)]
end
@propagate_inbounds Base.getindex(A::OffsetArray, i::Int)  = parent(A)[i]

@inline function Base.setindex!(A::OffsetArray{T,N}, val, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    J = @inbounds map(parentindex, axes(A), I)
    @inbounds parent(A)[J...] = val
    A
end

@inline function Base.setindex!(A::OffsetVector, val, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[parentindex(Base.axes1(A), i)] = val
    A
end
@propagate_inbounds function Base.setindex!(A::OffsetArray, val, i::Int)
    parent(A)[i] = val
    A
end

# For fast broadcasting: ref https://discourse.julialang.org/t/why-is-there-a-performance-hit-on-broadcasting-with-offsetarrays/32194
Base.dataids(A::OffsetArray) = Base.dataids(parent(A))
Broadcast.broadcast_unalias(dest::OffsetArray, src::OffsetArray) = parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src)

### Special handling for AbstractRange

const OffsetRange{T} = OffsetArray{T,1,<:AbstractRange{T}}
const IIUR = IdentityUnitRange{S} where S<:AbstractUnitRange{T} where T<:Integer

Base.step(a::OffsetRange) = step(parent(a))

@propagate_inbounds Base.getindex(a::OffsetRange, r::OffsetRange) = OffsetArray(a[parent(r)], r.offsets)
@propagate_inbounds function Base.getindex(a::OffsetRange, r::IdOffsetRange)
    OffsetArray(a.parent[r.parent .+ (r.offset - a.offsets[1])], r.offset)
end
@propagate_inbounds Base.getindex(r::OffsetRange, s::IIUR) =
    OffsetArray(r[s.indices], s)
@propagate_inbounds Base.getindex(a::OffsetRange, r::AbstractRange) = a.parent[r .- a.offsets[1]]
@propagate_inbounds Base.getindex(a::AbstractRange, r::OffsetRange) = OffsetArray(a[parent(r)], r.offsets)

@propagate_inbounds Base.getindex(r::UnitRange, s::IIUR) =
    OffsetArray(r[s.indices], s)

@propagate_inbounds Base.getindex(r::StepRange, s::IIUR) =
    OffsetArray(r[s.indices], s)

# this method is needed for ambiguity resolution
@propagate_inbounds Base.getindex(r::StepRangeLen{T,<:Base.TwicePrecision,<:Base.TwicePrecision}, s::IIUR) where T =
    OffsetArray(r[s.indices], s)

@propagate_inbounds Base.getindex(r::StepRangeLen{T}, s::IIUR) where {T} =
    OffsetArray(r[s.indices], s)

@propagate_inbounds Base.getindex(r::LinRange, s::IIUR) =
    OffsetArray(r[s.indices], s)

function Base.show(io::IO, r::OffsetRange)
    show(io, r.parent)
    o = r.offsets[1]
    print(io, " with indices ", o+1:o+length(r))
end
Base.show(io::IO, ::MIME"text/plain", r::OffsetRange) = show(io, r)

### Some mutating functions defined only for OffsetVector ###

Base.resize!(A::OffsetVector, nl::Integer) = (resize!(A.parent, nl); A)
Base.push!(A::OffsetVector, x...) = (push!(A.parent, x...); A)
Base.pop!(A::OffsetVector) = pop!(A.parent)
Base.append!(A::OffsetVector, items) = (append!(A.parent, items); A)
Base.empty!(A::OffsetVector) = (empty!(A.parent); A)

# These functions keep the summary compact
function Base.inds2string(inds::Tuple{Vararg{Union{IdOffsetRange,IdentityUnitRange{<:IdOffsetRange}}}})
    Base.inds2string(map(UnitRange, inds))
end
Base.showindices(io::IO, ind1::IdOffsetRange, inds::IdOffsetRange...) = Base.showindices(io, map(UnitRange, (ind1, inds...))...)

function Base.showarg(io::IO, a::OffsetArray, toplevel)
    print(io, "OffsetArray(")
    Base.showarg(io, parent(a), false)
    Base.showindices(io, axes(a)...)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(a))
end

function Base.replace_in_print_matrix(A::OffsetArray{<:Any,2}, i::Integer, j::Integer, s::AbstractString)
    J = map(parentindex, axes(A), (i,j))
    Base.replace_in_print_matrix(parent(A), J..., s)
end
function Base.replace_in_print_matrix(A::OffsetArray{<:Any,1}, i::Integer, j::Integer, s::AbstractString)
    ip = parentindex(axes(A,1), i)
    Base.replace_in_print_matrix(parent(A), ip, j, s)
end

function no_offset_view(A::AbstractArray)
    if Base.has_offset_axes(A)
        OffsetArray(A, map(r->1-first(r), axes(A)))
    else
        A
    end
end

no_offset_view(A::OffsetArray) = no_offset_view(parent(A))

end # module

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
function Euler1(N::Int64)

	sum = 0;
	for i in 1:(N-1)
		if (mod(i,3)==0) || (mod(i,5)==0)
			sum+=i
			println(i)
		end
	end

	return sum
end

@time Euler1(1000)
println(Euler1(1000))

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
## --- Dealing with different number representations

    """
    ```julia
    nearest(T, x)
    ```
    Convert `x` to the nearest representable value in type T, rounding if inexact
    """
    function nearest(::Type{T}, x) where T <: Number
        T(max(min(x, typemax(T)), typemin(T)))
    end
    function nearest(::Type{T}, x) where T <: Integer
        round(T, max(min(x, typemax(T)), typemin(T)))
    end
    export nearest

## --- Fast inverse square-root

    """
    ```julia
    fast_inv_sqrt(x)
    ```
    The infamous fast inverse square root of `x`, in 32 and 64 bit versions.
    Can be up to 10x faster than base `1/sqrt(x)`, though with nontrivial loss
    of precision. The implementations here are good to about 4 ppm.
    """
    function fast_inv_sqrt(x::Float64)
        x_2 = 0.5 * x
        result = Base.sub_int(9.603007803048109e153, Base.lshr_int(x,1)) # Floating point magic
        result *= ( 1.5 - (x_2 * result * result )) # Newton's method
        result *= ( 1.5 - (x_2 * result * result )) # Newton's method (again)
        return result
    end
    function fast_inv_sqrt(x::Float32)
        x_2 = 0.5f0 * x
        result = Base.sub_int(1.321202f19, Base.lshr_int(x,1)) # Floating point magic
        result *= ( 1.5f0 - (x_2 * result * result) ) # Newton's method
        result *= ( 1.5f0 - (x_2 * result * result) ) # Newton's method (again)
        return result
    end
    export fast_inv_sqrt

## --- Fast vectorized mean
    function _vmean(A::AbstractArray{T}, ::Colon) where T
        m = zero(T)
        @turbo for i ∈ eachindex(A)
            m += A[i]
        end
        return m / length(A)
    end

## --- Base-10 version of log1p
#
#     log10f(x::Number,from::Number=-1) = log10(abs(x-from)+1)*sign(x-from)
#

## --- Some mathematical constants

    const SQRT2 = sqrt(2)
    const SQRT2PI = sqrt(2*pi)
    const INVSQRT2 = 1/sqrt(2)
    const AN = Union{AbstractArray{<:Number},Number}

## --- Gaussian distribution functions

    """
    ```julia
    normpdf(mu,sigma,x)
    ```
    Probability density function of the Normal (Gaussian) distribution

    ``ℯ^{-(x-μ)^2 / (2σ^2)} / σ√2π``

    with mean `mu` and standard deviation `sigma`, evaluated at `x`
    """
    @inline normpdf(mu,sigma,x) = exp(-(x-mu)*(x-mu) / (2*sigma*sigma)) / (SQRT2PI*sigma)
    @inline normpdf(mu::Number,sigma::Number,x::Number) = exp(-(x-mu)*(x-mu) / (2*sigma*sigma)) / (SQRT2PI*sigma)
    normpdf(mu::AN,sigma::AN,x::AN) = @turbo @. exp(-(x-mu)*(x-mu) / (2*sigma*sigma)) / (SQRT2PI*sigma)
    export normpdf

    """
    ```julia
    normpdf_ll(mu, sigma, x)
    ```
    Fast log likelihood corresponding to the natural logarithm of the probability
    density function of a Normal (Gaussian) distribution with mean `mu` and
    standard deviation `sigma`, evaluated at `x`.

    If `x`, [`mu`, and `sigma`] are given as arrays, the sum of the log likelihood
    over all `x` will be returned.

    See also `normpdf`
    """
    @inline normpdf_ll(mu,sigma,x) = -(x-mu)*(x-mu) / (2*sigma*sigma)
    @inline normpdf_ll(mu::Number,sigma::Number,x::Number) = -(x-mu)*(x-mu) / (2*sigma*sigma)
    function normpdf_ll(mu::Number,sigma::Number,x::AbstractArray)
        inv_s2 = 1/(2*sigma*sigma)
        ll = zero(typeof(inv_s2))
        @turbo for i ∈ eachindex(x)
            ll -= (x[i]-mu)*(x[i]-mu) * inv_s2
        end
        return ll
    end
    function normpdf_ll(mu::AbstractArray,sigma::Number,x::AbstractArray)
        inv_s2 = 1/(2*sigma*sigma)
        ll = zero(typeof(inv_s2))
        @turbo for i ∈ eachindex(x)
            ll -= (x[i]-mu[i])*(x[i]-mu[i]) * inv_s2
        end
        return ll
    end
    function normpdf_ll(mu::Number,sigma::AbstractArray,x::AbstractArray)
        ll = zero(float(eltype(sigma)))
        @turbo for i ∈ eachindex(x)
            ll -= (x[i]-mu)*(x[i]-mu) / (2*sigma[i]*sigma[i])
        end
        return ll
    end
    function normpdf_ll(mu::AbstractArray,sigma::AbstractArray,x::AbstractArray)
        ll = zero(float(eltype(sigma)))
        @turbo for i ∈ eachindex(x)
            ll -= (x[i]-mu[i])*(x[i]-mu[i]) / (2*sigma[i]*sigma[i])
        end
        return ll
    end
    export normpdf_ll

    """
    ```julia
    normcdf(mu,sigma,x)
    ```
    Cumulative distribution function of the Normal (Gaussian) distribution

    ``1/2 + erf(\frac{x-μ}{σ√2})/2``

    with mean `mu` and standard deviation `sigma`, evaluated at `x`.
    """
    @inline normcdf(mu,sigma,x) = 0.5 + 0.5 * erf((x-mu) / (sigma*SQRT2))
    @inline normcdf(mu::Number,sigma::Number,x::Number) = 0.5 + 0.5 * erf((x-mu) / (sigma*SQRT2))
    normcdf(mu::AN,sigma::AN,x::AN) = @turbo @. 0.5 + 0.5 * erf((x-mu) / (sigma*SQRT2))
    export normcdf

    """
    ```julia
    normcdf_ll(mu, sigma, x)
    ```
    Fast log likelihood corresponding to the natural logarithm of the cumulative
    distribution function of a Normal (Gaussian) distribution with mean `mu` and
    standard deviation `sigma`, evaluated at `x`.

    If `x`, [`mu`, and `sigma`] are given as arrays, the sum of the log likelihood
    over all `x` will be returned.

    See also `normcdf`
    """
    @inline function normcdf_ll(xₛ::Number)
        if xₛ < -1.0
            return log(0.5*erfcx(-xₛ * INVSQRT2)) - 0.5*abs2(xₛ)
        else
            return log1p(-0.5*erfc(xₛ * INVSQRT2))
        end
    end
    function normcdf_ll(xₛ::AbstractArray)
        ll = zero(float(eltype(xₛ)))
        @inbounds for i ∈ eachindex(xₛ)
            ll += normcdf_ll(xₛ[i])
        end
        return ll
    end
    @inline function normcdf_ll(mu::Number, sigma::Number, x::Number)
        xₛ = (x - mu) / sigma
        return normcdf_ll(xₛ)
    end
    function normcdf_ll(mu::Number,sigma::Number,x::AbstractArray)
        inv_sigma = 1/sigma
        ll = zero(typeof(inv_sigma))
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu) * inv_sigma
            ll += normcdf_ll(xₛ)
        end
        return ll
    end
    function normcdf_ll(mu::AbstractArray,sigma::Number,x::AbstractArray)
        inv_sigma = 1/sigma
        ll = zero(typeof(inv_sigma))
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu[i]) * inv_sigma
            ll += normcdf_ll(xₛ)
        end
        return ll
    end
    function normcdf_ll(mu::Number,sigma::AbstractArray,x::AbstractArray)
        ll = zero(float(eltype(sigma)))
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu) / sigma[i]
            ll += normcdf_ll(xₛ)
        end
        return ll
    end
    function normcdf_ll(mu::AbstractArray,sigma::AbstractArray,x::AbstractArray)
        ll = zero(float(eltype(sigma)))
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu[i]) / sigma[i]
            ll += normcdf_ll(xₛ)
        end
        return ll
    end
    export normcdf_ll

    """
    ```julia
    normcdf_ll!(mu, sigma, x)
    ```
    Fast log likelihood corresponding to the natural logarithm of the cumulative
    distribution function of a Normal (Gaussian) distribution with mean `mu` and
    standard deviation `sigma`, evaluated at `x`.

    As `normcdf_ll`, but in-place (using `x` as a buffer).
    """
    function normcdf_ll!(xₛ::AbstractArray)
        @inbounds for i ∈ eachindex(xₛ)
            xₛ[i] = normcdf_ll(xₛ[i])
        end
        ll = zero(float(eltype(xₛ)))
        @turbo for i ∈ eachindex(xₛ)
            ll += xₛ[i]
        end
        return ll
    end
    function normcdf_ll!(mu::Number,sigma::Number,x::AbstractArray)
        inv_sigma = 1/sigma
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu) * inv_sigma
            x[i] = normcdf_ll(xₛ)
        end
        ll = zero(typeof(inv_sigma))
        @turbo for i ∈ eachindex(x)
            ll += x[i]
        end
        return ll
    end
    function normcdf_ll!(mu::AbstractArray,sigma::Number,x::AbstractArray)
        inv_sigma = 1/sigma
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu[i]) * inv_sigma
            x[i] = normcdf_ll(xₛ)
        end
        ll = zero(typeof(inv_sigma))
        @turbo for i ∈ eachindex(x)
            ll += x[i]
        end
        return ll
    end
    function normcdf_ll!(mu::Number,sigma::AbstractArray,x::AbstractArray)
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu) / sigma[i]
            x[i] = normcdf_ll(xₛ)
        end
        ll = zero(float(eltype(sigma)))
        @turbo for i ∈ eachindex(x)
            ll += x[i]
        end
        return ll
    end
    function normcdf_ll!(mu::AbstractArray,sigma::AbstractArray,x::AbstractArray)
        @inbounds for i ∈ eachindex(x)
            xₛ = (x[i] - mu[i]) / sigma[i]
            x[i] = normcdf_ll(xₛ)
        end
        ll = zero(float(eltype(sigma)))
        @turbo for i ∈ eachindex(x)
            ll += x[i]
        end
        return ll
    end
    export normcdf_ll!

    """
    ```julia
    normcdf!(result,mu,sigma,x)
    ```
    In-place version of `normcdf`
    """
    function normcdf!(result::Array, mu::Number, sigma::Number, x::AbstractArray)
        T = eltype(result)
        inv_sigma_sqrt2 = one(T)/(sigma*T(SQRT2))
        @turbo for i ∈ 1:length(x)
            result[i] = T(0.5) + T(0.5) * erf((x[i]-mu) * inv_sigma_sqrt2)
        end
        return result
    end
    export normcdf!

    """
    ```julia
    norm_quantile(F::Number)
    ```
    How far away from the mean (in units of sigma) should we expect proportion
    F of the samples to fall in a standard Gaussian (Normal[0,1]) distribution
    """
    @inline norm_quantile(F) = SQRT2*erfinv(2*F-1)
    export norm_quantile

    """
    ```julia
    norm_width(N::Number)
    ```
    How dispersed (in units of sigma) should we expect a sample of N numbers
    drawn from a standard Gaussian (Normal[0,1]) distribution to be?
    """
    @inline norm_width(N) = 2*norm_quantile(1 - 1/(2N))
    export norm_width

    """
    ```julia
    normproduct(μ1, σ1, μ2, σ2)
    ```
    The integral of the product of two normal distributions N[μ1,σ1] * N[μ2,σ2].
    This is itself just another Normal distribution! Specifically, one with
    variance σ1^2 + σ2^2, evaluated at distance |μ1-μ2| from the mean
    """
    normproduct(μ1, σ1, μ2, σ2) = normpdf(μ1, sqrt.(σ1.*σ1 + σ2.*σ2), μ2)
    export normproduct

    """
    ```julia
    normproduct_ll(μ1, σ1, μ2, σ2)
    ```
    Log likelihood corresponding to the integral of N[μ1,σ1] * N[μ2,σ2]
    As `normproduct`, but using the fast log likelihood of a Normal distribution
    """
    normproduct_ll(μ1, σ1, μ2, σ2) = normpdf_ll(μ1, sqrt.(σ1.*σ1 + σ2.*σ2), μ2)
    export normproduct_ll

## --- Geometry

    """
    ```julia
    inpolygon(x,y,point)
    ```
    Check if a 2D polygon defined by the arrays `x`, `y` contains a given `point`.
    Returns boolean (true or false)
    """
    function inpolygon(x,y,point)
        # Check that we have the right kind of input data
        if length(x) != length(y)
            error("polygon must have equal number of x and y points\n")
        end
        if length(x) < 3
            error("polygon must have at least 3 points\n")
        end
        if length(point) != 2
            error("point must be an ordered pair (x,y)\n")
        end

        # Extract x and y data of point
        point_x = point[1]
        point_y = point[2]
        # For first point, previous point is last
        x_here = x[end]
        y_here = y[end]
        # Ensure we are not sitting parallel to a vertex by infinitessimally moving the point
        if y_here == point_y
            point_y = nextfloat(float(point_y))
        end
        if x_here == point_x
            point_x = nextfloat(float(point_x))
        end

        # Check how many times a line projected right along x-axis from point intersects the polygon
        intersections = 0
        @inbounds for i ∈ eachindex(x)
            # Recycle our vertex
            x_last = copy(x_here)
            y_last = copy(y_here)

            # Get a new vertex
            x_here = x[i]
            y_here = y[i]

            # Ensure we are not sitting parallel to a vertex by infinitessimally moving the point
            if y_here == point_y
                point_y = nextfloat(float(point_y))
            end
            if x_here == point_x
                point_x = nextfloat(float(point_x))
            end

            if y_last > point_y && y_here > point_y
                # If both ys above point, no intersection
                continue
            elseif y_last < point_y && y_here < point_y
                # If both ys below point, no intersection
                continue
            elseif x_last < point_x && x_here < point_x
                # If both x's left of point, no intersection
                continue
            elseif x_last > point_x && x_here > point_x
                # By elimination
                # We have one y above and y below our point
                # If both y's are right of line, then definite intersection
                intersections += 1
                continue
            else
                # By elimination
                # One y above and one y below
                # One x to the right and one x to the left
                # We must project
                dy = y_here - y_last
                if abs(dy) > 0
                    dx = x_here - x_last
                    inv_slope = dx / dy
                    x_proj = x_last + (point_y - y_last) * inv_slope
                    if x_proj > point_x
                        intersections += 1
                    end
                end
            end
        end

        # If number of intersections is odd, point is in the polygon
        return Bool(mod(intersections,2))
    end
    export inpolygon

    """
    ```julia
    (rows, columns) = find_grid_inpolygon(grid_x, grid_y, poly_x, poly_y)
    ```
    Find the indexes of grid points that fall within a polygon for a grid with
    cell centers given by grid_x (j-columns of grid) and grid_y (i-rows of grid).
    Returns a list of rows and columns in the polygon
    """
    function find_grid_inpolygon(grid_x, grid_y, poly_x, poly_y)
        # Check that we have the right kind of input data
        if length(poly_x) != length(poly_y)
            error("polygon must have equal number of x and y points\n")
        end
        if length(poly_x) < 3
            error("polygon must have at least 3 points\n")
        end

        # Find maximum x and y range of polygon
        (xmin, xmax) = extrema(poly_x)
        (ymin, ymax) = extrema(poly_y)

        # Find the matrix indices within the range of the polygon (if any)
        column_inrange = findall((grid_x .>= xmin) .& (grid_x .<= xmax))
        row_inrange = findall((grid_y .>= ymin) .& (grid_y .<= ymax))

        # Keep a list of matrix indexes in the polygon
        row = Array{Int}(undef,length(column_inrange) * length(row_inrange))
        column = Array{Int}(undef,length(column_inrange) * length(row_inrange))
        n = 0
        for j = 1:length(column_inrange)
            for i = 1:length(row_inrange)
                point = [grid_x[column_inrange[j]], grid_y[row_inrange[i]]]
                if inpolygon(poly_x, poly_y, point)
                    n += 1
                    row[n] = row_inrange[i]
                    column[n] = column_inrange[j]
                end
            end
        end

        return (row[1:n], column[1:n])
    end
    export find_grid_inpolygon

    """
    ```julia
    arcdistance(latᵢ,lonᵢ,lat,lon)
    ```
    Calculate the distance on a sphere between the point (`latᵢ`,`lonᵢ`) and any
    number of points in (`lat`,`lon`).
    Latitude and Longitude should be specified in decimal degrees
    """
    function arcdistance(latᵢ,lonᵢ,lat,lon)
        # Argument for acos()
        arg = sin.(latᵢ .* pi/180) .* sin.(lat .* pi/180) .+ cos.(latᵢ*pi/180) .* cos.(lat .* pi/180).*cos.((lonᵢ .- lon) .* pi/180)

        # Avoid domain errors from imprecise sine and cosine math
        arg[arg .> 1.0] .= 1.0
        arg[arg .< -1.0] .= -1.0

        # Calculate angular distance
        theta = 180/pi .* acos.(arg)
        return theta
    end
    export arcdistance

## --- Linear regression

    """
    ```julia
    (a,b) = linreg(x::AbstractVector, y::AbstractVector)
    ```

    Returns the coefficients for a simple linear least-squares regression of
    the form `y = a + bx`
    """
    linreg(x::AbstractVector, y::AbstractVector) = hcat(fill!(similar(x), 1), x) \ y
    export linreg

## --- End of File

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
## demo functions:

"""
    demo1(gridChoice::String)

Demonstrate basic functionalities (load grid, arithmetic, exchange, gradient,
etc.). Call sequence:

```
!isdir("GRID_LLC90") ? error("missing files") : nothing

(D,Dexch,Darr,DD)=MeshArrays.demo1("LLC90");
```
"""
function demo1(gridChoice::String)

    mygrid=GridSpec(gridChoice)

    D=mygrid.read(mygrid.path*"Depth.data",MeshArray(mygrid,mygrid.ioPrec))

    1000 .+ D
    D .+ 1000
    D+D
    D .- 1000
    1000 .- D
    D-D
    1000*D
    D*1000
    D*D
    D/1000
    1000 ./ D
    D/D

    Dexch=exchange(D,4)
    Darr=mygrid.write(D)
    DD=mygrid.read(Darr,D)
    DD .== D

    GridVariables=GridLoad(mygrid)

    (dFLDdx, dFLDdy)=gradient(GridVariables["YC"],GridVariables)
    (dFLDdxEx,dFLDdyEx)=exchange(dFLDdx,dFLDdy,4)

    H=GridVariables["hFacC"]
    Ha=mygrid.write(H)
    Hb=mygrid.read(Ha,H)

    println(typeof(H))
    println(size(H))
    println(H.fSize[1])
    println(size(H.f[1]))
    if ndims(H)==3&&size(H,3)>40
        println(typeof(view(H,:,:,40)))
    elseif ndims(H.f)==2&&size(H.f,2)>40
        println(typeof(view(H,:,40)))
    end

    return (D,Dexch,Darr,DD)

end

##

"""
    demo2()

Demonstrate higher level functions using `smooth`. Call sequence:

```
(Rini,Rend,DXCsm,DYCsm)=MeshArrays.demo2();

using Plots; plotlyjs()
include(joinpath(dirname(pathof(MeshArrays)),"Plots.jl"))
heatmap(Rend,title="smoothed noise",clims=(-0.5,0.5))
heatmap(Rini,title="raw noise",clims=(-0.5,0.5))

#gr(); contour(Rend,clims=(-0.5,0.5),levels=-0.5:0.1:0.5,fill=true)
```

"""
function demo2()

    #Pre-requisite: either load predefined grid using `demo1` or call `GridOfOnes`
    isdir("GRID_LLC90") ? GridVariables=GridLoad(GridSpec("LLC90")) : GridVariables=GridOfOnes("cs",6,100)

    (Rini,Rend,DXCsm,DYCsm)=demo2(GridVariables)
end

function demo2(GridVariables::Dict)

    mygrid=GridVariables["XC"].grid

    #initialize 2D field of random numbers
    tmp1=randn(Float32,Tuple(mygrid.ioSize))
    Rini=mygrid.read(tmp1,MeshArray(mygrid,Float32))

    #apply land mask
    if ndims(GridVariables["hFacC"].f[1])>2
        tmp1=mask(view(GridVariables["hFacC"],:,:,1),NaN,0)
    elseif ndims(GridVariables["hFacC"].f)>1
        #tmp1=mask(view(GridVariables["hFacC"],:,1),NaN,0)
        tmp1=similar(Rini)
        for i=1:length(tmp1.fIndex); tmp1[i]=GridVariables["hFacC"][i,1]; end;
        tmp1=mask(tmp1,NaN,0)
    else
        tmp1=mask(GridVariables["hFacC"],NaN,0)
    end
    msk=fill(1.,tmp1) + 0. *tmp1;
    Rini=msk*Rini;

    #specify smoothing length scales in x, y directions
    DXCsm=3*GridVariables["DXC"]; DYCsm=3*GridVariables["DYC"];

    #apply smoother
    Rend=smooth(Rini,DXCsm,DYCsm,GridVariables);

    return (Rini,Rend,DXCsm,DYCsm)

end

"""
    demo3()

Demonstrate ocean transport computations. Call sequence:

```
(UV,LC,Tr)=MeshArrays.demo3();
using Plots; plot(Tr/1e6,title="meridional transport")

using Plots; plotlyjs()
include(joinpath(dirname(pathof(MeshArrays)),"Plots.jl"))
heatmap(UV["U"],title="U comp.",clims=(-1e7,1e7))
heatmap(UV["V"],title="V comp.",clims=(-1e7,1e7))
```
"""
function demo3()

    mygrid=GridSpec("LLC90")
    GridVariables=GridLoad(mygrid)

    TrspX=mygrid.read(mygrid.path*"TrspX.bin",MeshArray(mygrid,Float32))
    TrspY=mygrid.read(mygrid.path*"TrspY.bin",MeshArray(mygrid,Float32))
    TauX=mygrid.read(mygrid.path*"TauX.bin",MeshArray(mygrid,Float32))
    TauY=mygrid.read(mygrid.path*"TauY.bin",MeshArray(mygrid,Float32))
    SSH=mygrid.read(mygrid.path*"SSH.bin",MeshArray(mygrid,Float32))

    (UV, LC, Tr)=demo3(TrspX,TrspY,GridVariables)

end

function demo3(U::MeshArray,V::MeshArray,GridVariables::Dict)

    LC=LatitudeCircles(-89.0:89.0,GridVariables)

    #UV=Dict("U"=>U,"V"=>V,"dimensions"=>["x","y","z","t"],"factors"=>["dxory","dz"])
    UV=Dict("U"=>U,"V"=>V,"dimensions"=>["x","y"])

    Tr=Array{Float64,1}(undef,length(LC));
    for i=1:length(LC)
       Tr[i]=ThroughFlow(UV,LC[i],GridVariables)
    end

    return UV, LC, Tr

end