text:
module DiracDeltasTests

using Test
using Gridap.TensorValues
using Gridap.Helpers
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.CellData

domain = (0,1,0,1)
cells = (30,30)
model = CartesianDiscreteModel(domain,cells)

Ω = Triangulation(model)

vfun(x) = x[1] + 3*x[2]
v = CellField(vfun,Ω)

degree = 2
δ = DiracDelta{1}(model,degree,tags=5)
@test sum(δ(v)) ≈ 0.5

#using Gridap.Visualization
#writevtk(get_triangulation(δ),"trian_d")
#writevtk(Ω,"trian")

δ = DiracDelta{0}(model,tags=2)
@test sum(δ(v)) ≈ 1

δ = DiracDelta{0}(model,tags=4)
@test sum(δ(v)) ≈ 4

@test sum(δ(3.0)) ≈ 3.0
@test sum(δ(x->2*x)) ≈ VectorValue(2,2)

δ = DiracDelta{0}(model,tags=[2,4])
@test sum(δ(v)) ≈ 5

#using Gridap
#
#order = 2
#reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
#V = FESpace(model,reffe,dirichlet_tags=[1,2,5])
#
#dΩ = Measure(Ω,2*order)
#
#δ = DiracDelta{0}(model,tags=4)
#
#const E = 1
#const ν = 0.33
#const λ = (E*ν)/((1+ν)*(1-2*ν))
#const μ = E/(2*(1+ν))
#σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
#
#a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
#l(v) = δ( VectorValue(-1,0)⋅v )
#
#op = AffineFEOperator(a,l,V,V)
#uh = solve(op)
#
#writevtk(Ω,"results",cellfields=["uh"=>uh])
#writevtk(δ.Γ,"d",cellfields=["uh"=>uh])

end # module


text:
"""
** Connectance of a network**

Returns the connectance of a square matrix, defined as ``S/L^2``.

"""
function connectance(S::Int64, L::Int64)
    @assert S > 1
    C = L/S^2
    @assert C < 1.0
    @assert C > 0.0
    return C
end

function connectance(A::Array{Int64, 2})
    return connectance(size(A, 1), sum(A))
end

"""
**Niche model of food webs**

Takes a number of species `S` and a number of interactions `L`, and returns
a food web with predators in rows, and preys in columns. This function is
used internally by `nichemodel` called with a connectance.

"""
function nichemodel(S::Int64, L::Int64)

    C = connectance(S, L)

    # Beta distribution parameter
    β = 1/(2*C)-1.0

    # Pre-allocate the network
    A = zeros(Int64, (S, S))

    # Generate body size
    n = sort(rand(Uniform(0.0, 1.0), S))

    # Pre-allocate centroids
    c = zeros(Float64, S)

    # Generate random ranges
    r = n .* rand(Beta(1.0, β), S)

    # Generate random centroids
    for s in 1:S
        c[s] = rand(Uniform(r[s]/2, n[s]))
    end

    # The smallest species has a body size and range of 0
    n[n.==minimum(n)] .= 0.0
    r[n.==minimum(n)] .= 0.0

    for consumer in 1:S
        for resource in 1:S
            if n[resource] < c[consumer] + r[consumer]
                if n[resource] > c[consumer] - r[consumer]
                    A[consumer, resource] = 1
                end
            end
        end
    end

    return A

end

"""
**Niche model of food webs**

Takes a number of species `S` and a connectance `C`, and returns a food web
with predators in rows, and preys in columns. Note that the connectance is
first transformed into an integer number of interactions.

This function has two keyword arguments:

1. `tolerance` is the allowed error on tolerance (see below)

2. `toltype` is the type or error, and can be `:abs` (absolute) and `:rel`
(relative). Relative tolerance is the amount of error allowed, relative to the
desired connectance value. If the simulated network has a tolerance x, the
target connectance is c, then the relative error is |1-x/c|.

"""
function nichemodel(S::Int64, C::Float64; tolerance::Float64=0.05, toltype::Symbol=:abs)
    @assert C < 1.0
    @assert C > 0.0
    @assert tolerance > 0.0
    @assert toltype ∈ [:abs, :rel]
    L = round(Int64, C * S^2)
    A = nichemodel(S, L)
    if toltype == :abs
      tolfunc = (x) -> abs(x-C) < tolerance
    else
      tolfunc = (x) -> abs(1-x/C) < tolerance
    end
    while !(tolfunc(BioEnergeticFoodWebs.connectance(A)))
      A = nichemodel(S, L)
    end
    return A
end

ext: jl

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


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

module RangesTest

using Test
using Dates

using InteractiveUtils: subtypes

let
    for T in (Dates.Date, Dates.DateTime)
        f1 = T(2014); l1 = T(2013, 12, 31)
        f2 = T(2014); l2 = T(2014)
        f3 = T(-2000); l3 = T(2000)
        f4 = typemin(T); l4 = typemax(T)

        for P in subtypes(Dates.DatePeriod)
            for pos_step in (P(1), P(2), P(50), P(2048), P(10000))
                # empty range
                dr = f1:pos_step:l1
                len = length(dr)
                @test len == 0
                @test isa(len, Int64)
                @test isempty(dr)
                @test first(dr) == f1
                @test last(dr) < f1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr), dr) == Int64[]
                @test [dr;] == T[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) < f1
                @test last(reverse(dr)) >= f1
                @test issorted(dr)
                @test sortperm(dr) == 1:1:0
                @test !(f1 in dr)
                @test !(l1 in dr)
                @test !(f1 - pos_step in dr)
                @test !(l1 + pos_step in dr)
                @test dr == []
                @test hash(dr) == hash([])

                for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                    dr = f:pos_step:l
                    len = length(dr)
                    @test len > 0
                    @test isa(len, Int64)
                    @test !isempty(dr)
                    @test first(dr) == f
                    @test last(dr) <= l
                    @test minimum(dr) == first(dr)
                    @test maximum(dr) == last(dr)
                    @test dr[1] == f
                    @test dr[end] <= l
                    @test iterate(dr) == (first(dr), (length(dr), 1))

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr), dr) == [1:len;]
                        @test length([dr;]) == len
                        @test dr == dr1
                        @test hash(dr) == hash(dr1)
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test last(reverse(dr)) == f
                    @test issorted(dr)
                    @test f in dr

                end
            end
            for neg_step in (P(-1), P(-2), P(-50), P(-2048), P(-10000))
                # empty range
                dr = l1:neg_step:f1
                len = length(dr)
                @test len == 0
                @test isa(len, Int64)
                @test isempty(dr)
                @test first(dr) == l1
                @test last(dr) > l1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr), dr) == Int64[]
                @test [dr;] == T[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) > l1
                @test last(reverse(dr)) <= l1
                @test issorted(dr)
                @test sortperm(dr) == 1:1:0
                @test !(l1 in dr)
                @test !(l1 in dr)
                @test !(l1 - neg_step in dr)
                @test !(l1 + neg_step in dr)
                @test dr == []
                @test hash(dr) == hash([])

                for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                    dr = l:neg_step:f
                    len = length(dr)
                    @test len > 0
                    @test isa(len, Int64)
                    @test !isempty(dr)
                    @test first(dr) == l
                    @test last(dr) >= f
                    @test minimum(dr) == last(dr)
                    @test maximum(dr) == first(dr)
                    @test dr[1] == l
                    @test dr[end] >= f
                    @test iterate(dr) == (first(dr), (length(dr), 1))

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr), dr) == [1:len;]
                        @test length([dr;]) == len
                        @test dr == dr1
                        @test hash(dr) == hash(dr1)
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test issorted(dr) == (len <= 1)
                    @test l in dr
                end
            end
        end
        if T == Dates.DateTime
            for P in subtypes(Dates.TimePeriod)
                P in (Dates.Microsecond, Dates.Nanosecond) && continue
                for pos_step in (P(1), P(2), P(50), P(2048), P(10000))
                    # empty range
                    dr = f1:pos_step:l1
                    len = length(dr)
                    @test len == 0
                    @test isa(len, Int64)
                    @test isempty(dr)
                    @test first(dr) == f1
                    @test last(dr) < f1
                    @test length([i for i in dr]) == 0
                    @test_throws ArgumentError minimum(dr)
                    @test_throws ArgumentError maximum(dr)
                    @test_throws BoundsError dr[1]
                    @test findall(in(dr), dr) == Int64[]
                    @test [dr;] == T[]
                    @test isempty(reverse(dr))
                    @test length(reverse(dr)) == 0
                    @test first(reverse(dr)) < f1
                    @test last(reverse(dr)) >= f1
                    @test issorted(dr)
                    @test sortperm(dr) == 1:1:0
                    @test !(f1 in dr)
                    @test !(l1 in dr)
                    @test !(f1 - pos_step in dr)
                    @test !(l1 + pos_step in dr)
                    @test dr == []
                    @test hash(dr) == hash([])

                    for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                        dr = f:pos_step:l
                        len = length(dr)
                        @test len > 0
                        @test isa(len, Int64)
                        @test !isempty(dr)
                        @test first(dr) == f
                        @test last(dr) <= l
                        @test minimum(dr) == first(dr)
                        @test maximum(dr) == last(dr)
                        @test dr[1] == f
                        @test dr[end] <= l
                        @test iterate(dr) == (first(dr), (length(dr), 1))

                        if len < 10000
                            dr1 = [i for i in dr]
                            @test length(dr1) == len
                            @test findall(in(dr), dr) == [1:len;]
                            @test length([dr;]) == len
                            @test dr == dr1
                            @test hash(dr) == hash(dr1)
                        end
                        @test !isempty(reverse(dr))
                        @test length(reverse(dr)) == len
                        @test last(reverse(dr)) == f
                        @test issorted(dr)
                        @test f in dr

                    end
                end
                for neg_step in (P(-1), P(-2), P(-50), P(-2048), P(-10000))
                    # empty range
                    dr = l1:neg_step:f1
                    len = length(dr)
                    @test len == 0
                    @test isa(len, Int64)
                    @test isempty(dr)
                    @test first(dr) == l1
                    @test last(dr) > l1
                    @test length([i for i in dr]) == 0
                    @test_throws ArgumentError minimum(dr)
                    @test_throws ArgumentError maximum(dr)
                    @test_throws BoundsError dr[1]
                    @test findall(in(dr), dr) == Int64[]
                    @test [dr;] == T[]
                    @test isempty(reverse(dr))
                    @test length(reverse(dr)) == 0
                    @test first(reverse(dr)) > l1
                    @test last(reverse(dr)) <= l1
                    @test issorted(dr)
                    @test sortperm(dr) == 1:1:0
                    @test !(l1 in dr)
                    @test !(l1 in dr)
                    @test !(l1 - neg_step in dr)
                    @test !(l1 + neg_step in dr)
                    @test dr == []
                    @test hash(dr) == hash([])

                    for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                        dr = l:neg_step:f
                        len = length(dr)
                        @test len > 0
                        @test isa(len, Int64)
                        @test !isempty(dr)
                        @test first(dr) == l
                        @test last(dr) >= f
                        @test minimum(dr) == last(dr)
                        @test maximum(dr) == first(dr)
                        @test dr[1] == l
                        @test dr[end] >= f
                        @test iterate(dr) == (first(dr), (length(dr), 1))

                        if len < 10000
                            dr1 = [i for i in dr]
                            @test length(dr1) == len
                            @test findall(in(dr), dr) == [1:len;]
                            @test length([dr;]) == len
                            @test dr == dr1
                            @test hash(dr) == hash(dr1)
                        end
                        @test !isempty(reverse(dr))
                        @test length(reverse(dr)) == len
                        @test last(reverse(dr)) >= l
                        @test issorted(dr) == (len <= 1)
                        @test l in dr

                    end
                end
            end
        end
    end
end

# Dates are physical units, and ranges should require an explicit step.
# See #19896 and https://discourse.julialang.org/t/type-restriction-on-unitrange/6557/12
if VERSION >= v"1.0.0"
    @test_throws MethodError Dates.DateTime(2013, 1, 1):Dates.DateTime(2013, 2, 1)
end

# All the range representations we want to test
# Date ranges
dr  = Dates.DateTime(2013, 1, 1):Dates.Day(1):Dates.DateTime(2013, 2, 1)
dr1 = Dates.DateTime(2013, 1, 1):Dates.Day(1):Dates.DateTime(2013, 1, 1)
dr2 = Dates.DateTime(2013, 1, 1):Dates.Day(1):Dates.DateTime(2012, 2, 1) # empty range
dr3 = Dates.DateTime(2013, 1, 1):Dates.Day(-1):Dates.DateTime(2012) # negative step
# Big ranges
dr4 = Dates.DateTime(0):Dates.Day(1):Dates.DateTime(20000, 1, 1)
dr5 = Dates.DateTime(0):Dates.Day(1):Dates.DateTime(200000, 1, 1)
dr6 = Dates.DateTime(0):Dates.Day(1):Dates.DateTime(2000000, 1, 1)
dr7 = Dates.DateTime(0):Dates.Day(1):Dates.DateTime(20000000, 1, 1)
dr8 = Dates.DateTime(0):Dates.Day(1):Dates.DateTime(200000000, 1, 1)
dr9 = typemin(Dates.DateTime):Dates.Day(1):typemax(Dates.DateTime)
# Other steps
dr10 = typemax(Dates.DateTime):Dates.Day(-1):typemin(Dates.DateTime)
dr11 = typemin(Dates.DateTime):Dates.Week(1):typemax(Dates.DateTime)

dr12 = typemin(Dates.DateTime):Dates.Month(1):typemax(Dates.DateTime)
dr13 = typemin(Dates.DateTime):Dates.Year(1):typemax(Dates.DateTime)

dr14 = typemin(Dates.DateTime):Dates.Week(10):typemax(Dates.DateTime)
dr15 = typemin(Dates.DateTime):Dates.Month(100):typemax(Dates.DateTime)
dr16 = typemin(Dates.DateTime):Dates.Year(1000):typemax(Dates.DateTime)
dr17 = typemax(Dates.DateTime):Dates.Week(-10000):typemin(Dates.DateTime)
dr18 = typemax(Dates.DateTime):Dates.Month(-100000):typemin(Dates.DateTime)
dr19 = typemax(Dates.DateTime):Dates.Year(-1000000):typemin(Dates.DateTime)
dr20 = typemin(Dates.DateTime):Dates.Day(2):typemax(Dates.DateTime)

drs = Any[dr, dr1, dr2, dr3, dr4, dr5, dr6, dr7, dr8, dr9, dr10,
          dr11, dr12, dr13, dr14, dr15, dr16, dr17, dr18, dr19, dr20]
drs2 = map(x->Dates.Date(first(x)):step(x):Dates.Date(last(x)), drs)

@test map(length, drs) == map(x->size(x)[1], drs)
@test map(length, drs) == map(x->length(Dates.Date(first(x)):step(x):Dates.Date(last(x))), drs)
@test map(length, drs) == map(x->length(reverse(x)), drs)
@test all(x->findall(in(x), x)==[1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == range(last(x), step=-step(x), length=length(x)), drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test all(drs[1:3]) do dd
    for (i, d) in enumerate(dd)
        @test d == (first(dd) + Dates.Day(i - 1))
    end
    true
end
@test_throws MethodError dr .+ 1
a = Dates.DateTime(2013, 1, 1)
b = Dates.DateTime(2013, 2, 1)
@test map!(x->x + Dates.Day(1), Vector{Dates.DateTime}(undef, 32), dr) == [(a + Dates.Day(1)):Dates.Day(1):(b + Dates.Day(1));]
@test map(x->x + Dates.Day(1), dr) == [(a + Dates.Day(1)):Dates.Day(1):(b + Dates.Day(1));]

@test map(x->a in x, drs[1:4]) == [true, true, false, true]
@test a in dr
@test b in dr
@test Dates.DateTime(2013, 1, 3) in dr
@test Dates.DateTime(2013, 1, 15) in dr
@test Dates.DateTime(2013, 1, 26) in dr
@test !(Dates.DateTime(2012, 1, 1) in dr)

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test length(b:Dates.Day(-1):a) == 32
@test length(b:Dates.Day(1):a) == 0
@test length(b:Dates.Day(1):a) == 0
@test length(a:Dates.Day(2):b) == 16
@test last(a:Dates.Day(2):b) == Dates.DateTime(2013, 1, 31)
@test length(a:Dates.Day(7):b) == 5
@test last(a:Dates.Day(7):b) == Dates.DateTime(2013, 1, 29)
@test length(a:Dates.Day(32):b) == 1
@test last(a:Dates.Day(32):b) == Dates.DateTime(2013, 1, 1)
@test (a:Dates.Day(1):b)[1] == Dates.DateTime(2013, 1, 1)
@test (a:Dates.Day(1):b)[2] == Dates.DateTime(2013, 1, 2)
@test (a:Dates.Day(1):b)[7] == Dates.DateTime(2013, 1, 7)
@test (a:Dates.Day(1):b)[end] == b
@test first(a:Dates.Day(1):Dates.DateTime(20000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.DateTime(200000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.DateTime(2000000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.DateTime(20000000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.DateTime(200000000, 1, 1)) == a
@test first(a:Dates.Day(1):typemax(Dates.DateTime)) == a
@test first(typemin(Dates.DateTime):Dates.Day(1):typemax(Dates.DateTime)) == typemin(Dates.DateTime)

# Date ranges
dr  = Dates.Date(2013, 1, 1):Dates.Day(1):Dates.Date(2013, 2, 1)
dr1 = Dates.Date(2013, 1, 1):Dates.Day(1):Dates.Date(2013, 1, 1)
dr2 = Dates.Date(2013, 1, 1):Dates.Day(1):Dates.Date(2012, 2, 1) # empty range
dr3 = Dates.Date(2013, 1, 1):Dates.Day(-1):Dates.Date(2012, 1, 1) # negative step
# Big ranges
dr4 = Dates.Date(0):Dates.Day(1):Dates.Date(20000, 1, 1)
dr5 = Dates.Date(0):Dates.Day(1):Dates.Date(200000, 1, 1)
dr6 = Dates.Date(0):Dates.Day(1):Dates.Date(2000000, 1, 1)
dr7 = Dates.Date(0):Dates.Day(1):Dates.Date(20000000, 1, 1)
dr8 = Dates.Date(0):Dates.Day(1):Dates.Date(200000000, 1, 1)
dr9 = typemin(Dates.Date):Dates.Day(1):typemax(Dates.Date)
# Other steps
dr10 = typemax(Dates.Date):Dates.Day(-1):typemin(Dates.Date)
dr11 = typemin(Dates.Date):Dates.Week(1):typemax(Dates.Date)
dr12 = typemin(Dates.Date):Dates.Month(1):typemax(Dates.Date)
dr13 = typemin(Dates.Date):Dates.Year(1):typemax(Dates.Date)
dr14 = typemin(Dates.Date):Dates.Week(10):typemax(Dates.Date)
dr15 = typemin(Dates.Date):Dates.Month(100):typemax(Dates.Date)
dr16 = typemin(Dates.Date):Dates.Year(1000):typemax(Dates.Date)
dr17 = typemax(Dates.Date):Dates.Week(-10000):typemin(Dates.Date)
dr18 = typemax(Dates.Date):Dates.Month(-100000):typemin(Dates.Date)
dr19 = typemax(Dates.Date):Dates.Year(-1000000):typemin(Dates.Date)
dr20 = typemin(Dates.Date):Dates.Day(2):typemax(Dates.Date)

drs = Any[dr, dr1, dr2, dr3, dr4, dr5, dr6, dr7, dr8, dr9, dr10,
          dr11, dr12, dr13, dr14, dr15, dr16, dr17, dr18, dr19, dr20]

@test map(length, drs) == map(x->size(x)[1], drs)
@test all(x->findall(in(x), x) == [1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == last(x): - step(x):first(x), drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test all(drs[1:3]) do dd
    for (i, d) in enumerate(dd)
        @test d == (first(dd) + Dates.Day(i - 1))
    end
    true
end
@test_throws MethodError dr .+ 1
a = Dates.Date(2013, 1, 1)
b = Dates.Date(2013, 2, 1)
@test map!(x->x + Dates.Day(1), Vector{Dates.Date}(undef, 32), dr) == [(a + Dates.Day(1)):Dates.Day(1):(b + Dates.Day(1));]
@test map(x->x + Dates.Day(1), dr) == [(a + Dates.Day(1)):Dates.Day(1):(b + Dates.Day(1));]

@test map(x->a in x, drs[1:4]) == [true, true, false, true]
@test a in dr
@test b in dr
@test Dates.Date(2013, 1, 3) in dr
@test Dates.Date(2013, 1, 15) in dr
@test Dates.Date(2013, 1, 26) in dr
@test !(Dates.Date(2012, 1, 1) in dr)

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test length(b:Dates.Day(-1):a) == 32
@test length(b:Dates.Day(1):a) == 0
@test length(a:Dates.Day(2):b) == 16
@test last(a:Dates.Day(2):b) == Dates.Date(2013, 1, 31)
@test length(a:Dates.Day(7):b) == 5
@test last(a:Dates.Day(7):b) == Dates.Date(2013, 1, 29)
@test length(a:Dates.Day(32):b) == 1
@test last(a:Dates.Day(32):b) == Dates.Date(2013, 1, 1)
@test (a:Dates.Day(1):b)[1] == Dates.Date(2013, 1, 1)
@test (a:Dates.Day(1):b)[2] == Dates.Date(2013, 1, 2)
@test (a:Dates.Day(1):b)[7] == Dates.Date(2013, 1, 7)
@test (a:Dates.Day(1):b)[end] == b
@test first(a:Dates.Day(1):Dates.Date(20000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.Date(200000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.Date(2000000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.Date(20000000, 1, 1)) == a
@test first(a:Dates.Day(1):Dates.Date(200000000, 1, 1)) == a
@test first(a:Dates.Day(1):typemax(Dates.Date)) == a
@test first(typemin(Dates.Date):Dates.Day(1):typemax(Dates.Date)) == typemin(Dates.Date)

@test length(typemin(Dates.Date):Dates.Week(1):typemax(Dates.Date)) == 26351950414948059
# Big Month/Year ranges
@test length(typemin(Dates.Date):Dates.Month(1):typemax(Dates.Date)) == 6060531933867600
@test length(typemin(Dates.Date):Dates.Year(1):typemax(Dates.Date)) == 505044327822300
@test length(typemin(Dates.DateTime):Dates.Month(1):typemax(Dates.DateTime)) == 3507324288
@test length(typemin(Dates.DateTime):Dates.Year(1):typemax(Dates.DateTime)) == 292277024

@test length(typemin(Dates.DateTime):Dates.Week(1):typemax(Dates.DateTime)) == 15250284420
@test length(typemin(Dates.DateTime):Dates.Day(1):typemax(Dates.DateTime)) == 106751990938
@test length(typemin(Dates.DateTime):Dates.Hour(1):typemax(Dates.DateTime)) == 2562047782512
@test length(typemin(Dates.DateTime):Dates.Minute(1):typemax(Dates.DateTime)) == 153722866950720
@test length(typemin(Dates.DateTime):Dates.Second(1):typemax(Dates.DateTime)) == 9223372017043200
@test length(typemin(DateTime):Dates.Millisecond(1):typemax(DateTime)) == 9223372017043199001

c = Dates.Date(2013, 6, 1)
@test length(a:Dates.Month(1):c) == 6
@test [a:Dates.Month(1):c;] == [a + Dates.Month(1)*i for i in 0:5]
@test [a:Dates.Month(2):Dates.Date(2013, 1, 2);] == [a]
@test [c:Dates.Month(-1):a;] == reverse([a:Dates.Month(1):c;])

@test length(range(Date(2000), step=Dates.Day(1), length=366)) == 366
let n = 100000
    local a = Dates.Date(2000)
    for i = 1:n
        @test length(range(a, step=Dates.Day(1), length=i)) == i
    end
    return a + Dates.Day(n)
end

@test typeof(step(Dates.DateTime(2000):Dates.Day(1):Dates.DateTime(2001))) == Dates.Day

a = Dates.Date(2013, 1, 1)
b = Dates.Date(2013, 2, 1)
d = Dates.Date(2020, 1, 1)
@test length(a:Dates.Year(1):d) == 8
@test first(a:Dates.Year(1):d) == a
@test last(a:Dates.Year(1):d) == d
@test length(a:Dates.Month(12):d) == 8
@test first(a:Dates.Month(12):d) == a
@test last(a:Dates.Month(12):d) == d
@test length(a:Dates.Week(52):d) == 8
@test first(a:Dates.Week(52):d) == a
@test last(a:Dates.Week(52):d) == Dates.Date(2019, 12, 24)
@test length(a:Dates.Day(365):d) == 8
@test first(a:Dates.Day(365):d) == a
@test last(a:Dates.Day(365):d) == Dates.Date(2019, 12, 31)

a = Dates.Date(2013, 1, 1)
b = Dates.Date(2013, 2, 1)
@test length(a:Dates.Year(1):Dates.Date(2020, 2, 1)) == 8
@test length(a:Dates.Year(1):Dates.Date(2020, 6, 1)) == 8
@test length(a:Dates.Year(1):Dates.Date(2020, 11, 1)) == 8
@test length(a:Dates.Year(1):Dates.Date(2020, 12, 31)) == 8
@test length(a:Dates.Year(1):Dates.Date(2021, 1, 1)) == 9
@test length(Dates.Date(2000):Dates.Year(-10):Dates.Date(1900)) == 11
@test length(Dates.Date(2000, 6, 23):Dates.Year(-10):Dates.Date(1900, 2, 28)) == 11
@test length(Dates.Date(2000, 1, 1):Dates.Year(1):Dates.Date(2000, 2, 1)) == 1

let n = 100000
    local a, b
    a= b = Dates.Date(0)
    for i = 1:n
        @test length(a:Dates.Year(1):b) == i
        b += Dates.Year(1)
    end
end

let n = 10000,
    a = Dates.Date(1985, 12, 5),
    b = Dates.Date(1986, 12, 27),
    c = Dates.DateTime(1985, 12, 5),
    d = Dates.DateTime(1986, 12, 27)
    for i = 1:n
        @test length(a:Dates.Month(1):b) == 13
        @test length(a:Dates.Year(1):b) == 2
        @test length(c:Dates.Month(1):d) == 13
        @test length(c:Dates.Year(1):d) == 2
        a += Dates.Day(1)
        b += Dates.Day(1)
    end
end

let n = 100000
    local a, b
    a = b = Dates.Date(2000)
    for i = 1:n
        @test length(a:Dates.Month(1):b) == i
        b += Dates.Month(1)
    end
end

@test length(Dates.Year(1):Dates.Year(1):Dates.Year(10)) == 10
@test length(Dates.Year(10):Dates.Year(-1):Dates.Year(1)) == 10
@test length(Dates.Year(10):Dates.Year(-2):Dates.Year(1)) == 5
@test_throws OverflowError length(typemin(Dates.Year):Dates.Year(1):typemax(Dates.Year))
@test_throws MethodError Dates.Date(0):Dates.DateTime(2000)
@test_throws MethodError Dates.Date(0):Dates.Year(10)
@test length(range(Dates.Date(2000), step=Dates.Day(1), length=366)) == 366
@test last(range(Dates.Date(2000), step=Dates.Day(1), length=366)) == Dates.Date(2000, 12, 31)
@test last(range(Dates.Date(2001), step=Dates.Day(1), length=365)) == Dates.Date(2001, 12, 31)
@test last(range(Dates.Date(2000), step=Dates.Day(1), length=367)) == last(range(Dates.Date(2000), step=Dates.Month(12), length=2)) == last(range(Dates.Date(2000), step=Dates.Year(1), length=2))
@test last(range(Dates.DateTime(2000), step=Dates.Day(366), length=2)) == last(range(Dates.DateTime(2000), step=Dates.Hour(8784), length=2))

# Issue 5
lastdaysofmonth = [Dates.Date(2014, i, Dates.daysinmonth(2014, i)) for i=1:12]
@test [Date(2014, 1, 31):Dates.Month(1):Date(2015);] == lastdaysofmonth

# Range addition/subtraction:
let d = Dates.Day(1)
    @test (Dates.Date(2000):d:Dates.Date(2001)) + d == (Dates.Date(2000) + d:d:Dates.Date(2001) + d)
    @test (Dates.Date(2000):d:Dates.Date(2001)) - d == (Dates.Date(2000) - d:d:Dates.Date(2001) - d)
end

# Time ranges
dr  = Dates.Time(23, 1, 1):Dates.Second(1):Dates.Time(23, 2, 1)
dr1 = Dates.Time(23, 1, 1):Dates.Second(1):Dates.Time(23, 1, 1)
dr2 = Dates.Time(23, 1, 1):Dates.Second(1):Dates.Time(22, 2, 1) # empty range
dr3 = Dates.Time(23, 1, 1):Dates.Minute(-1):Dates.Time(22, 1, 1) # negative step
# Big ranges
dr8 = typemin(Dates.Time):Dates.Second(1):typemax(Dates.Time)
dr9 = typemin(Dates.Time):Dates.Nanosecond(1):typemax(Dates.Time)
# Other steps
dr10 = typemax(Dates.Time):Dates.Microsecond(-1):typemin(Dates.Time)
dr11 = typemin(Dates.Time):Dates.Millisecond(1):typemax(Dates.Time)
dr12 = typemin(Dates.Time):Dates.Minute(1):typemax(Dates.Time)
dr13 = typemin(Dates.Time):Dates.Hour(1):typemax(Dates.Time)
dr14 = typemin(Dates.Time):Dates.Millisecond(10):typemax(Dates.Time)
dr15 = typemin(Dates.Time):Dates.Minute(100):typemax(Dates.Time)
dr16 = typemin(Dates.Time):Dates.Hour(1000):typemax(Dates.Time)
dr17 = typemax(Dates.Time):Dates.Millisecond(-10000):typemin(Dates.Time)
dr18 = typemax(Dates.Time):Dates.Minute(-100):typemin(Dates.Time)
dr19 = typemax(Dates.Time):Dates.Hour(-10):typemin(Dates.Time)
dr20 = typemin(Dates.Time):Dates.Microsecond(2):typemax(Dates.Time)

drs = Any[dr, dr1, dr2, dr3, dr8, dr9, dr10,
          dr11, dr12, dr13, dr14, dr15, dr16, dr17, dr18, dr19, dr20]

@test map(length, drs) == map(x->size(x)[1], drs)
@test all(x->findall(in(x), x) == [1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == last(x): - step(x):first(x), drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test_throws MethodError dr .+ 1

a = Dates.Time(23, 1, 1)
@test map(x->a in x, drs[1:4]) == [true, true, false, true]
@test a in dr

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test !(1 in Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05))
@test !(Complex(1, 0) in Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05))
@test !(π in Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05))
@test !("a" in Date(2017, 01, 01):Dates.Day(1):Date(2017, 01, 05))

@test hash(Any[Date("2018-1-03"), Date("2018-1-04"), Date("2018-1-05")]) ==
      hash([Date("2018-1-03"), Date("2018-1-04"), Date("2018-1-05")]) ==
      hash(Date("2018-1-03"):Day(1):Date("2018-1-05"))

end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
# Tests for entity-based learning
function t_entity_networklearner()

#########################################
# Test the entity-based NetworkLearner  #
#########################################
N= 100							# Number of entitites
inferences = [:ic, :rl, :gs, :unknown]			# Collective inferences
rlearners = [:rn, :wrn, :bayesrn, :unknown]		# Relational learners
nAdj = 2						# Number of adjacencies to generate
X = rand(1,N); 						# Training data

nlmodel=[]

#####################
# Column-major case #
#####################

# Initializations
ft=x->vec(x)
Xo = rand(1,N)
update = trues(N)

# Train and test methods for relational model
fr_train=(x)->sum(x[1], dims=2);
fr_exec=(m,x)->sum(x.-m, dims=1)

amv = sparse.(Symmetric.([sprand(Float64, N,N, 0.5) for i in 1:nAdj]));
adv = adjacency.(amv);

for infopt in inferences
	for rlopt in rlearners
		Test.@test try
			# Train NetworkLearner
			nlmodel=fit(NetworkLearnerEnt, Xo, update,
			       adv, fr_train,fr_exec;
			       learner=rlopt,
			       inference=infopt,
			       f_targets=ft,
			       normalize=false, maxiter = 5
			)

			# Make modifications of adjacencies, estimates
			# ...

			# Re-run inference
			infer!(nlmodel)
			true
		catch
			false
		end
	end
end

##################
# Row-major case #
##################

# Initializations
ft=x->vec(x)
Xo = rand(N,1)
update = trues(N)

# Train and test methods for relational model
fr_train=(x)->sum(x[1], dims=1);
fr_exec=(m,x)->sum(x.-m, dims=2)

amv = sparse.(Symmetric.([sprand(Float64, N,N, 0.5) for i in 1:nAdj]));
adv = adjacency.(amv);

for infopt in inferences
	for rlopt in rlearners
		Test.@test try
			# Train NetworkLearner
			nlmodel=fit(NetworkLearnerEnt, Xo, update,
			       adv, fr_train,fr_exec;
			       learner=rlopt,
			       inference=infopt,
			       f_targets=ft,
			       normalize=false, maxiter = 5,
			       obsdim=1
			)

			# Make modifications of adjacencies, estimates
			# ...

			# Re-run inference
			infer!(nlmodel)
			true
		catch
			false
		end
	end
end

# Test show methods
buf = IOBuffer()
Test.@test try
	show(buf,nlmodel)
	true
catch
	false
end

Test.@test try
	show(buf,NetworkLearning.NetworkLearnerState(rand(2,2), trues(2), ObsDim.Constant{2}()))
	true
catch
	false
end

end


text:
using Test
using SymDiff

@test differentiate(:(x / 1), :x) == :((1 * 1 + x * 0) / 1 ^ 2)


text:
function LVGNA_alignments(testing=true)

    files = readdir(sparse_matrices_path)

    if testing
        files = files[6:7]
        files =  ["worm_PHY1.smat","worm_Y2H1.smat"]
        result_path = all_results_path*"LVGNA_alignments/test/"
        alphas = [1.0]
        betas = [0.0]
        samples = 300000
    else
        result_path = all_results_path*"LVGNA_alignments/"
        alphas = [.5,1.0]
        betas = [0.0,1.0,1e1,1e2]
        samples = 3000000
    end

    orders = [3]
        # generates the tensors using TuranShadow as smat
        # is needed for post processing. Use ssten files for
        # previously generated triangle adjacency tensors, but
        # post processing isn't supported.

        # _MultiMotif_M() flag is used to generate tensors,
        # replace with _M() for non-generating routines.

    if !isdir(result_path)
        mkpath(result_path)
    end

    #
    #  --  LowRankTAME  --  #
    #
    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=KlauAlgo(),method=LowRankTAME_MultiMotif_M(),iter=15,
                                    alphas,betas,profile=true,no_matching=false,motif=Clique(),orders,samples=samples)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LowRankTAME_alphas:$(alphas)_betas:$(betas)_iter:15_order:3_postProcessing:KlauAlgo_profile:true_samples:$(samples)_tol:1e-6_results.json")
    save_results(result_path*filename,X)

    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=LocalSearch(),method=LowRankTAME_MultiMotif_M(),iter=15,
                                    alphas,betas,profile=true,no_matching=false,motif=Clique(),orders,samples=samples)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LowRankTAME_alphas:$(alphas)_betas:$(betas)_iter:15_order:$(orders)_postProcessing:LocalSearch_profile:true_tol:1e-6_results.json")
    save_results(result_path*filename,X)

    #=   Low Rank Matching Routines
    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=KlauAlgo(),method=LowRankTAME_MultiMotif_M(),iter=15,
                                    alphas,betas,profile=true,no_matching=false,motif=Clique(),orders,samples=samples,low_rank_matching = true)
    filename =  filter(x-> !isspace(x),LVGNA_pairwiseAlignment_LowRankTAME_lrm_alphas:$(alphas)_betas:$(betas)_iter:15_order:3_postProcessing:KlauAlgo_profile:true_samples:$(samples)_tol:1e-6_results.json")
    save_results(result_path*filename,X)

    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=LocalSearch(),method=LowRankTAME_MultiMotif_M(),iter=15,
                                    alphas,betas,profile=true,no_matching=false,motif=Clique(),orders,samples=samples,low_rank_matching = true)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LowRankTAME_lrm_alphas:$(alphas)_betas:$(betas)_iter:15_order:3_postProcessing:LocalSearch_profile:true_samples:$(samples)_tol:1e-6_results.json")
    save_results(result_path*filename,X)
    =#

    #
    #  --  Lambda - TAME  --  #
    #
    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=KlauAlgo(),method=ΛTAME_MultiMotif_M(),iter=15,matchingMethod=ΛTAME_GramMatching(),
                                    alphas,betas,profile=true,no_matching=false,orders,samples=samples)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LambdaTAME_alphas:$(alphas)_betas:$(betas)_iter:15_matchingMethod:GramMatching_orders:$(orders)_postProcessing:KlauAlgo_profile:true_tol:1e-6_results.json")
    save_results(result_path*filename,X)

    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=LocalSearch(),method=ΛTAME_MultiMotif_M(),iter=15,matchingMethod=ΛTAME_GramMatching(),
                                    alphas,betas,profile=true,no_matching=false,motif=Clique(),orders,samples=samples)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LambdaTAME_alphas:$(alphas)_betas:$(betas)_iter:15_matchingMethod:GramMatching_orders:$(orders)_postProcessing:LocalSearch_profile:true_samples:$(samples)_tol:1e-6_results.json")
    save_results(result_path*filename,X)

    #= Low Rank Matching Routines
    #02/11/22 TODO: must fix samples and orders being passed in for calling rankOneMatching

    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=LocalSearch(),method=ΛTAME_M(),iter=15,matchingMethod=ΛTAME_rankOneMatching(),
                                    alphas,betas,profile=true,no_matching=false,motif=Clique(),orders=[3],samples=samples)
    filename =  filter(x-> !isspace(x)"LVGNA_pairwiseAlignment_LambdaTAME_alphas:$(alphas)_betas:$(betas)_iter:15_matchingMethod:rankOneMatching_order:3_postProcessing:LocalSearch_profile:true_samples:$(samples)_tol:1e-6_results.json")
    save_results(result_path*filename,X)
    =#

    #
    #  -- LowRankEigenAlign --  #
    #
    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=KlauAlgo(),method=LowRankEigenAlign_M(),profile=true)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LREA_postProcessing:KlauAlgo_profile:true_results.json")
    save_results(result_path*filename,X)

    X = distributed_pairwise_alignment(files,sparse_matrices_path;postProcessing=LocalSearch(),method=LowRankEigenAlign_M(),profile=true)
    filename =  filter(x-> !isspace(x),"LVGNA_pairwiseAlignment_LREA_postProcessing:LocalSearch_profile:true_results.json")
    save_results(result_path*filename,X)

    #=
    #
    #   (Julia) TAME implementation
    #
    X = distributed_pairwise_alignment(files,sparse_matrices_path;method=TAME_M(),iter=15,
                                       alphas = [1.0,.5],betas = [0.0,1e1,1e2,1e3],
                                       profile=true,no_matching=false)
    =#

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
module CartesianFDM

using Base.Iterators

using LinearAlgebra
using SparseArrays

using Reexport
@reexport using Symbolics

const subscripts = ("\u2081", "\u2082", "\u2083")
const TupleN{T,N} = NTuple{N,T}

export scalar, vector

export Periodic, periodic
export NonPeriodic, nonperiodic

export Dirichlet, dir
export Neumann, neu
export Mixed

export star

export CartesianFDMContext, cartesianfdmcontext

export mask, gradient, divergence, strainrate, divergence2, dissipation, permanent

export linearize

export spacing, coordinate

export potentialflow

include("symbolics.jl")
include("topology.jl")
include("boundary.jl")
include("stencil.jl")
include("operators.jl")
include("calculus.jl")
include("linearization.jl")
include("mesh.jl")
include("flow.jl")

end


text:
# Autogenerated wrapper script for ERFA_jll for i686-w64-mingw32
export liberfa

JLLWrappers.@generate_wrapper_header("ERFA")
JLLWrappers.@declare_library_product(liberfa, "liberfa-1.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liberfa,
        "bin\\liberfa-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ASTInterpreter2, Test

function read_and_parse(filename)
    src = read(filename, String)
    ex = Base.parse_input_line(src; filename=filename)
end

module Toplevel end

@testset "toplevel" begin
    stack = ASTInterpreter2.JuliaStackFrame[]
    ASTInterpreter2.interpret!(stack, Toplevel, read_and_parse("toplevel_script.jl"))

    @test isconst(Toplevel, :StructParent)
    @test isconst(Toplevel, :Struct)
    @test isconst(Toplevel, :MyInt8)

    s = Toplevel.Struct([2.0])

    @test Toplevel.f1(0) == 1
    @test Toplevel.f1(0.0) == 2
    @test Toplevel.f1(0.0f0) == 3
    @test Toplevel.f2("hi") == -1
    @test Toplevel.f2(UInt16(1)) == UInt16
    @test Toplevel.f2(3.2) == 0
    @test Toplevel.f2(view([1,2], 1:1)) == 2
    @test Toplevel.f2([1,2]) == 3
    @test Toplevel.f2(reshape(view([1,2], 1:2), 2, 1)) == 4
    if VERSION >= v"1.1.0"
        @test Toplevel.f3(1, 1) == 1
    end
    @test Toplevel.f3(1, :hi) == 2
    @test Toplevel.f3(UInt16(1), :hi) == Symbol
    @test Toplevel.f3(rand(2, 2), :hi, :there) == 2
    @test_throws MethodError Toplevel.f3([1.0], :hi, :there)
    @test Toplevel.f4(1, 1.0) == 1
    @test Toplevel.f4(1, 1) == Toplevel.f4(1) == 2
    @test Toplevel.f4(UInt(1), "hey", 2) == 3
    @test Toplevel.f4(rand(2,2)) == 2
    @test Toplevel.f5(Int8(1); y=22) == 22
    @test Toplevel.f5(Int16(1)) == 2
    @test Toplevel.f5(Int32(1)) == 3
    @test Toplevel.f5(Int64(1)) == 4
    @test Toplevel.f5(rand(2,2); y=7) == 2
    @test Toplevel.f6(1, "hi"; z=8) == 1
    @test Toplevel.f7(1, (1, :hi)) == 1
    @test Toplevel.f8(0) == 1
    @test Toplevel.f9(3) == 9
    @test Toplevel.f9(3.0) == 3.0
    @test s("hello") == [2.0]
    @test Toplevel.Struct{Float32}(Dict(1=>"two")) == 4
    @test Toplevel.first_two_funcs == (Toplevel.f1, Toplevel.f2)
    if VERSION >= v"1.2.0-DEV.239"
        @test isconst(Toplevel, :first_two_funcs)
    else
        @test_broken isconst(Toplevel, :first_two_funcs)
    end
    @test Toplevel.myint isa Toplevel.MyInt8
    @test_throws UndefVarError Toplevel.ffalse(1)
    @test Toplevel.ftrue(1) == 3
    @test Toplevel.fctrue(0) == 1
    @test_throws UndefVarError Toplevel.fcfalse(0)
    @test !Toplevel.Consts.b2
    @test Toplevel.fb1true(0) == 1
    @test_throws UndefVarError Toplevel.fb1false(0)
    @test Toplevel.fb2false(0) == 1
    @test_throws UndefVarError Toplevel.fb2true(0)
    @test Toplevel.fstrue(0) == 1
    @test Toplevel.fouter(1) === 2
    @test Toplevel.feval1(1.0) === 1
    @test Toplevel.feval1(1.0f0) === 1
    @test_throws MethodError Toplevel.feval1(1)
    @test Toplevel.feval2(1.0, Int8(1)) == 2
    @test length(s) === nothing
    @test size(s) === nothing
    @test Toplevel.nbytes(Float32) == 4
    @test Toplevel.typestring(1.0) == "Float64"
    @test Toplevel._feval3(0) == 3
    @test Toplevel.feval_add!(0) == 1
    @test Toplevel.feval_min!(0) == 1
    @test Toplevel.paramtype(Vector{Int8}) == Int8
    @test Toplevel.paramtype(Vector) == Toplevel.NoParam
    @test Toplevel.Inner.g() == 5
    @test Toplevel.Inner.InnerInner.g() == 6

    @test @interpret(Toplevel.f1(0)) == 1
    @test @interpret(Toplevel.f1(0.0)) == 2
    @test @interpret(Toplevel.f1(0.0f0)) == 3
    @test @interpret(Toplevel.f2("hi")) == -1
    @test @interpret(Toplevel.f2(UInt16(1))) == UInt16
    @test @interpret(Toplevel.f2(3.2)) == 0
    @test @interpret(Toplevel.f2(view([1,2], 1:1))) == 2
    @test @interpret(Toplevel.f2([1,2])) == 3
    @test @interpret(Toplevel.f2(reshape(view([1,2], 1:2), 2, 1))) == 4
    if VERSION >= v"1.1.0"
        @test @interpret(Toplevel.f3(1, 1)) == 1
    end
    @test @interpret(Toplevel.f3(1, :hi)) == 2
    @test @interpret(Toplevel.f3(UInt16(1), :hi)) == Symbol
    @test @interpret(Toplevel.f3(rand(2, 2), :hi, :there)) == 2
    @test_throws ErrorException("no unique matching method found for the specified argument types") @interpret(Toplevel.f3([1.0], :hi, :there))
    @test @interpret(Toplevel.f4(1, 1.0)) == 1
    @test @interpret(Toplevel.f4(1, 1)) == @interpret(Toplevel.f4(1)) == 2
    @test @interpret(Toplevel.f4(UInt(1), "hey", 2)) == 3
    @test @interpret(Toplevel.f4(rand(2,2))) == 2
    @test @interpret(Toplevel.f5(Int8(1); y=22)) == 22
    @test @interpret(Toplevel.f5(Int16(1))) == 2
    @test @interpret(Toplevel.f5(Int32(1))) == 3
    @test @interpret(Toplevel.f5(Int64(1))) == 4
    @test @interpret(Toplevel.f5(rand(2,2); y=7)) == 2
    @test @interpret(Toplevel.f6(1, "hi"; z=8)) == 1
    @test @interpret(Toplevel.f7(1, (1, :hi))) == 1
    @test @interpret(Toplevel.f8(0)) == 1
    @test @interpret(Toplevel.f9(3)) == 9
    @test @interpret(Toplevel.f9(3.0)) == 3.0
    @test @interpret(s("hello")) == [2.0]
    @test @interpret(Toplevel.Struct{Float32}(Dict(1=>"two"))) == 4
    @test_throws UndefVarError @interpret(Toplevel.ffalse(1))
    @test @interpret(Toplevel.ftrue(1)) == 3
    @test @interpret(Toplevel.fctrue(0)) == 1
    @test_throws UndefVarError @interpret(Toplevel.fcfalse(0))
    @test @interpret(Toplevel.fb1true(0)) == 1
    @test_throws UndefVarError @interpret(Toplevel.fb1false(0))
    @test @interpret(Toplevel.fb2false(0)) == 1
    @test_throws UndefVarError @interpret(Toplevel.fb2true(0))
    @test @interpret(Toplevel.fstrue(0)) == 1
    @test @interpret(Toplevel.fouter(1)) === 2
    @test @interpret(Toplevel.feval1(1.0)) === 1
    @test @interpret(Toplevel.feval1(1.0f0)) === 1
    @test_throws ErrorException("no unique matching method found for the specified argument types") @interpret(Toplevel.feval1(1))
    @test @interpret(Toplevel.feval2(1.0, Int8(1))) == 2
    @test @interpret(length(s)) === nothing
    @test @interpret(size(s)) === nothing
    @test @interpret(Toplevel.nbytes(Float32)) == 4
    @test @interpret(Toplevel.typestring(1.0)) == "Float64"
    @test @interpret(Toplevel._feval3(0)) == 3
    @test @interpret(Toplevel.feval_add!(0)) == 1
    @test @interpret(Toplevel.feval_min!(0)) == 1
    @test @interpret(Toplevel.paramtype(Vector{Int8})) == Int8
    @test @interpret(Toplevel.paramtype(Vector)) == Toplevel.NoParam
    @test @interpret(Toplevel.Inner.g()) == 5
    @test @interpret(Toplevel.Inner.InnerInner.g()) == 6

    # Check that nested expressions are handled appropriately (module-in-block, internal `using`)
    ex = quote
       module Testing
       if true
           using ASTInterpreter2
       end
       end
   end
   ASTInterpreter2.interpret!(stack, Toplevel, ex)
   @test Toplevel.Testing.JuliaStackFrame === JuliaStackFrame
end


text:
#########################
### Forest fire model ###
#########################
using Agents
using Random

mutable struct Tree{T<:Integer} <: AbstractAgent
  id::T
  pos::Tuple{T, T}
  status::T  # 1 is green, 2 is burning, 3 is burned
end

mutable struct Forest{T<:AbstractSpace, Y<:AbstractVector, Z<:AbstractFloat} <: AbstractModel
  space::T
  agents::Y
  scheduler::Function
  d::Z  # forest density
end

mutable struct MyGrid{T<:Integer, Y<:AbstractVector} <: AbstractSpace
  dimensions::Tuple{T, T}
  space::SimpleGraph
  agent_positions::Y
end

function model_initiation(;d, griddims, seed)
  Random.seed!(seed)
  # initialize the model
  # we start the model without creating the agents first
  agent_positions = [Int64[] for i in 1:gridsize(griddims)]
  mygrid = MyGrid(griddims, grid(griddims, false, false), agent_positions)
  forest = Forest(mygrid, Array{Tree}(undef, 0), random_activation, d)

  # create and add trees to each node with probability d, which determines the density of the forest. If a tree is on one edge of the grid (x=0), set it on fire.
  idcounter = 0
  for node in 1:gridsize(forest.space.dimensions)
    pp = rand()
    if pp <= forest.d
      idcounter += 1
      if vertex_to_coord(node, forest)[1] == 1
        tree = Tree(idcounter, (1,1), 2)
      else
        tree = Tree(idcounter, (1,1), 1)
      end
      add_agent!(tree, node, forest)
    end
  end
  return forest
end

function forest_step!(forest)
  shuffled_nodes = Random.shuffle(1:length(forest.agents))
  still_burning = false
  for node in shuffled_nodes  # randomly go through the trees
    tree = forest.agents[node]
    if tree.status == 2  # if it is has been burning
      still_burning = true
      tree.status = 3
      # set its neighbor trees on fire
      neighbor_cells = node_neighbors(tree, forest)
      for cell in neighbor_cells
        treeid = get_node_contents(cell, forest)
        if length(treeid) != 0 # the cell is not empty
          treen = forest.agents[treeid[1]]
          if treen.status == 1
            treen.status = 2
          end
        end
      end
    end
  end
  if still_burning == false
    return "All trees burned down."
  end
end


text:
@testset "Structured Matrices" begin
    @testset "/ and \\ on Square Matrixes" begin
        @testset "//, $T on the RHS" for T in (Diagonal, UpperTriangular, LowerTriangular)
            RHS = T(randn(T == Diagonal ? 10 : (10, 10)))
            Y = randn(5, 10)
            Ȳ = randn(size(/(Y, RHS))...)
            rrule_test(/, Ȳ, (Y, randn(size(Y))), (RHS, randn(size(RHS))))
        end

        @testset "\\ $T on LHS" for T in (Diagonal, UpperTriangular, LowerTriangular)
            LHS = T(randn(T == Diagonal ? 10 : (10, 10)))
            y = randn(10)
            ȳ = randn(size(\(LHS, y))...)
            rrule_test(\, ȳ, (LHS, randn(size(LHS))), (y, randn(10)))
            Y = randn(10, 10)
            Ȳ = randn(10, 10)
            rrule_test(\, Ȳ, (LHS, randn(size(LHS))), (Y, randn(size(Y))))
        end
    end

    @testset "Diagonal" begin
        N = 3
        rrule_test(Diagonal, randn(N, N), (randn(N), randn(N)))
        D = Diagonal(randn( N))
        rrule_test(Diagonal, D, (randn(N), randn(N)))
        # Concrete type instead of UnionAll
        rrule_test(typeof(D), D, (randn(N), randn(N)))

        # TODO: replace this with a `rrule_test` once we have that working
        # see https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/24
        res, pb = rrule(Diagonal, [1, 4])
        @test pb(10*res) == (NO_FIELDS, [10, 40])
        comp = Composite{typeof(res)}(; diag=10*res.diag)  # this is the structure of Diagonal
        @test pb(comp) == (NO_FIELDS, [10, 40])
    end

    @testset "::Diagonal * ::AbstractVector" begin
        N = 3
        rrule_test(
            *,
            randn(N),
            (Diagonal(randn(N)), Diagonal(randn(N))),
            (randn(N), randn(N)),
        )
    end
    @testset "diag" begin
        N = 7
        rrule_test(diag, randn(N), (randn(N, N), randn(N, N)))
        rrule_test(diag, randn(N), (Diagonal(randn(N)), randn(N, N)))
        rrule_test(diag, randn(N), (randn(N, N), Diagonal(randn(N))))
        rrule_test(diag, randn(N), (Diagonal(randn(N)), Diagonal(randn(N))))
        VERSION ≥ v"1.3" && @testset "k=$k" for k in (-1, 0, 2)
            M = N - abs(k)
            rrule_test(diag, randn(M), (randn(N, N), randn(N, N)), (k, nothing))
        end
    end
    @testset "diagm" begin
        @testset "without size" begin
            M, N = 7, 9
            s = (8, 8)
            a, ā = randn(M), randn(M)
            b, b̄ = randn(M), randn(M)
            c, c̄ = randn(M - 1), randn(M - 1)
            ȳ = randn(s)
            ps = (0 => a, 1 => b, 0 => c)
            y, back = rrule(diagm, ps...)
            @test y == diagm(ps...)
            ∂self, ∂pa, ∂pb, ∂pc = back(ȳ)
            @test ∂self === NO_FIELDS
            ∂a_fd, ∂b_fd, ∂c_fd = j′vp(_fdm, (a, b, c) -> diagm(0 => a, 1 => b, 0 => c), ȳ, a, b, c)
            for (p, ∂px, ∂x_fd) in zip(ps, (∂pa, ∂pb, ∂pc), (∂a_fd, ∂b_fd, ∂c_fd))
                ∂px = unthunk(∂px)
                @test ∂px isa Composite{typeof(p)}
                @test ∂px.first isa AbstractZero
                @test ∂px.second ≈ ∂x_fd
            end
        end
        VERSION ≥ v"1.3" && @testset "with size" begin
            M, N = 7, 9
            a, ā = randn(M), randn(M)
            b, b̄ = randn(M), randn(M)
            c, c̄ = randn(M - 1), randn(M - 1)
            ȳ = randn(M, N)
            ps = (0 => a, 1 => b, 0 => c)
            y, back = rrule(diagm, M, N, ps...)
            @test y == diagm(M, N, ps...)
            ∂self, ∂M, ∂N, ∂pa, ∂pb, ∂pc = back(ȳ)
            @test ∂self === NO_FIELDS
            @test ∂M === DoesNotExist()
            @test ∂N === DoesNotExist()
            ∂a_fd, ∂b_fd, ∂c_fd = j′vp(_fdm, (a, b, c) -> diagm(M, N, 0 => a, 1 => b, 0 => c), ȳ, a, b, c)
            for (p, ∂px, ∂x_fd) in zip(ps, (∂pa, ∂pb, ∂pc), (∂a_fd, ∂b_fd, ∂c_fd))
                ∂px = unthunk(∂px)
                @test ∂px isa Composite{typeof(p)}
                @test ∂px.first isa AbstractZero
                @test ∂px.second ≈ ∂x_fd
            end
        end
    end
    @testset "$(SymHerm)(::AbstractMatrix{$T}, :$(uplo))" for
        SymHerm in (Symmetric, Hermitian),
        T in (Float64, ComplexF64),
        uplo in (:U, :L)

        N = 3
        @testset "frule" begin
            x = randn(T, N, N)
            Δx = randn(T, N, N)
            # can't use frule_test here because it doesn't yet ignore nothing tangents
            Ω = SymHerm(x, uplo)
            Ω_ad, ∂Ω_ad = frule((Zero(), Δx, Zero()), SymHerm, x, uplo)
            @test Ω_ad == Ω
            ∂Ω_fd = jvp(_fdm, z -> SymHerm(z, uplo), (x, Δx))
            @test ∂Ω_ad ≈ ∂Ω_fd
        end
        @testset "rrule" begin
            x = randn(T, N, N)
            ∂x = randn(T, N, N)
            ΔΩ = randn(T, N, N)
            @testset "back(::$MT)" for MT in (Matrix, LowerTriangular, UpperTriangular)
                rrule_test(SymHerm, MT(ΔΩ), (x, ∂x), (uplo, nothing))
            end
            @testset "back(::Diagonal)" begin
                rrule_test(SymHerm, Diagonal(ΔΩ), (x, Diagonal(∂x)), (uplo, nothing))
            end
        end
    end
    @testset "$(f)(::$(SymHerm){$T}) with uplo=:$uplo" for f in (Matrix, Array),
        SymHerm in (Symmetric, Hermitian),
        T in (Float64, ComplexF64),
        uplo in (:U, :L)

        N = 3
        x = SymHerm(randn(T, N, N), uplo)
        Δx = randn(T, N, N)
        ∂x = SymHerm(randn(T, N, N), uplo)
        ΔΩ = f(SymHerm(randn(T, N, N), uplo))
        frule_test(f, (x, Δx))
        frule_test(f, (x, SymHerm(Δx, uplo)))
        rrule_test(f, ΔΩ, (x, ∂x))
    end
    @testset "$f" for f in (Adjoint, adjoint, Transpose, transpose)
        n = 5
        m = 3
        rrule_test(f, randn(m, n), (randn(n, m), randn(n, m)))
        rrule_test(f, randn(1, n), (randn(n), randn(n)))
    end
    @testset "$T" for T in (UpperTriangular, LowerTriangular)
        n = 5
        rrule_test(T, T(randn(n, n)), (randn(n, n), randn(n, n)))
    end
    @testset "$Op" for Op in (triu, tril)
        n = 7
        rrule_test(Op, randn(n, n), (randn(n, n), randn(n, n)))
        @testset "k=$k" for k in -2:2
            rrule_test(Op, randn(n, n), (randn(n, n), randn(n, n)), (k, nothing))
        end
    end
end


text:
function init_test_budget(lastWorkAge)
    ageEfficiencyV = collect(LinRange(1.0, 3.0, lastWorkAge));
    ageEfficiencyV[20 : lastWorkAge] .= 2.0;
    wage = 2.0;
    intRate = 1.05;
    taxRate = 0.4;
    retireTransfer = 3.9;
    cMin = 1e-3;
    b = Budget(ageEfficiencyV, wage, intRate, taxRate, retireTransfer, cMin);
    @assert validate_budget(b);
    return b
end

function validate_budget(b :: Budget)
    isValid = validate_age_efficiencies(b.ageEfficiencyV)  &&
        (wage(b) > 0.0)  &&  (b.intRate > 0.8)  &&  (b.intRate < 1.2)  &&
        (b.taxRate >= 0.0)  &&  (b.taxRate < 0.9);
    return isValid
end

function validate_age_efficiencies(ageEfficiencyV :: Vector{Float64})
    isValid = all(ageEfficiencyV .> 0.0);
    return isValid
end

## -------  Helpers

c_min(b :: Budget) = b.cMin;
# betar(b :: Budget) = b.discFactor * b.intRate;
wage(b :: Budget) = b.wage;
workspan(b) = length(b.ageEfficiencyV);
isretired(b :: Budget, t :: Integer) = t > workspan(b);

## ---------  Work phase

function earnings(b :: Budget, t :: Integer, efficiency)
    if isretired(b, t)
        earn = 0.0;
    else
        earn = wage(b) .* (1.0 - b.taxRate) .* efficiency .* b.ageEfficiencyV[t];
    end
    return earn
end

function retire_transfer(b :: Budget, t :: Integer)
    if isretired(b, t)
        transfer = b.retireTransfer;
    else
        transfer = 0.0;
    end
    return transfer
end

capital_income(b :: Budget, k) = b.intRate .* k;

income(b :: Budget, t :: Integer, k, efficiency) =
    earnings(b, t, efficiency) .+
    retire_transfer(b, t) .+
    capital_income(b, k);

kprime(b :: Budget, t :: Integer, k, c, efficiency) =
    income(b, t, k, efficiency) .- c;

consumption(b :: Budget, t :: Integer, k, kPrime, efficiency) =
    income(b, t, k, efficiency) .- kPrime;

# Max kPrime consistent with c > c_min
kprime_max(b :: Budget, t :: Integer, k, efficiency) =
    kprime(b, t, k, c_min(b), efficiency);

# ## ----------  Retirement
# # No efficiency state

# retired_income(b :: Budget, k) = b.retireTransfer .+ b.intRate .* k;
# retired_kprime(b :: Budget, k, c) = retired_income(b, k) .- c;
# retired_consumption(b :: Budget, k, kPrime) = retired_income(b, k) .- kPrime;

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


text:
function scalar_rhs!(du::AbstractMatrix, rhs1, f_face, f_interaction, dgl, dgr)
    ncell, nsp = size(du)
    @threads for ppp1 = 1:nsp
        for i = 2:ncell-1
            @inbounds du[i, ppp1] =
                -(
                    rhs1[i, ppp1] +
                    (f_interaction[i] - f_face[i, 1]) * dgl[ppp1] +
                    (f_interaction[i+1] - f_face[i, 2]) * dgr[ppp1]
                )
        end
    end

    return nothing
end

function scalar_rhs!(du::CuDeviceMatrix, rhs1, f_face, f_interaction, dgl, dgr)
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    idy = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    strx = blockDim().x * gridDim().x
    stry = blockDim().y * gridDim().y

    ncell, nsp = size(du)
    for ppp1 = idy:stry:nsp
        for i = idx+1:strx:ncell-1
            @inbounds du[i, ppp1] =
                -(
                    rhs1[i, ppp1] +
                    (f_interaction[i] - f_face[i, 1]) * dgl[ppp1] +
                    (f_interaction[i+1] - f_face[i, 2]) * dgr[ppp1]
                )
        end
    end

    return nothing
end


text:
abstract type HeuristicsCoStateInitializer{HOT <: Heuristics.Optimizer} end

function initialize!(hcsi::HeuristicsCoStateInitializer)
    optimize!(hcsi.ho, hcsi.hoOpts)
    hcsi.λh .= hcsi.ho.results.xbest
    return hcsi.ho.results
end

function GetInitializedCostates(hcsi::HeuristicsCoStateInitializer)
    return hcsi.λh
end


text:
function gencodeelemface(filename, npm, ncase)

    foldername = "app";
    opufile = "opu" * filename;
    cpufile = "cpu" * filename;
    gpufile = "gpu" * filename;

    tmp = "";
    for k = 1:npm
        tmp = tmp * "#include " * "\"" * opufile * string(k) *  ".cpp" * "\"" * "\n";
    end
    tmp = tmp * "\n";

    if ncase==1
        tmp = tmp * "template <typename T> void " * opufile;
        tmp = tmp * "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uinf, T *param, T time, int modelnumber, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw)\n";
        tmp = tmp * "{\n";
        for k = 1:npm
            if k == 1
                tmp = tmp * "\tif (modelnumber == " * string(k) * ")\n";
            else
                tmp = tmp * "\telse if (modelnumber == " * string(k) * ")\n";
            end
            tmp = tmp * "\t\t" * opufile * string(k) * "(f, xdg, udg, odg, wdg, uinf, param, time, modelnumber, ng, nc, ncu, nd, ncx, nco, ncw);\n";
        end
        tmp = tmp * "}\n\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int);\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int);\n";
    elseif ncase==2
        tmp = tmp * "template <typename T> void " * opufile;
        tmp = tmp * "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uinf, T *param, T time, int modelnumber, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw, int nce, int npe, int ne)\n";
        tmp = tmp * "{\n";
        for k = 1:npm
            if k == 1
                tmp = tmp * "\tif (modelnumber == " * string(k) * ")\n";
            else
                tmp = tmp * "\telse if (modelnumber == " * string(k) * ")\n";
            end
            tmp = tmp * "\t\t" * opufile * string(k) * "(f, xdg, udg, odg, wdg, uinf, param, time, modelnumber, ng, nc, ncu, nd, ncx, nco, ncw, nce, npe, ne);\n";
        end
        tmp = tmp * "}\n\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int, int, int, int);\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int, int, int, int);\n";
    elseif ncase==3
        tmp = tmp * "template <typename T> void " * opufile;
        tmp = tmp * "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uhg, T *nlg, T *tau, T *uinf, T *param, T time, int modelnumber, int ib, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw)\n";
        tmp = tmp * "{\n";
        for k = 1:npm
            if k == 1
                tmp = tmp * "\tif (modelnumber == " * string(k) * ")\n";
            else
                tmp = tmp * "\telse if (modelnumber == " * string(k) * ")\n";
            end
            tmp = tmp * "\t\t" * opufile * string(k) * "(f, xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time, modelnumber, ib, ng, nc, ncu, nd, ncx, nco, ncw);\n";
        end
        tmp = tmp * "}\n\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int, int);\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int, int);\n";
    elseif ncase==4
        tmp = tmp * "template <typename T> void " * opufile;
        tmp = tmp * "(T *f, T *xdg, T *udg1, T *udg2,  T *odg1, T *odg2,  T *wdg1, T *wdg2,  T *uhg, T *nlg, T *tau, T *uinf, T *param, T time, int modelnumber, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw)\n";
        tmp = tmp * "{\n";
        for k = 1:npm
            if k == 1
                tmp = tmp * "\tif (modelnumber == " * string(k) * ")\n";
            else
                tmp = tmp * "\telse if (modelnumber == " * string(k) * ")\n";
            end
            tmp = tmp * "\t\t" * opufile * string(k) * "(f, xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time, modelnumber, ng, nc, ncu, nd, ncx, nco, ncw);\n";
        end
        tmp = tmp * "}\n\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double time, int, int, int, int, int, int, int, int);\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float time, int, int, int, int, int, int, int, int);\n";
    elseif ncase==5
        tmp = tmp * "template <typename T> void " * opufile;
        tmp = tmp * "(T *f, T *xdg, T *uinf, T *param, int modelnumber, int ng, int ncx, int nce, int npe, int ne)\n";
        tmp = tmp * "{\n";
        for k = 1:npm
            if k == 1
                tmp = tmp * "\tif (modelnumber == " * string(k) * ")\n";
            else
                tmp = tmp * "\telse if (modelnumber == " * string(k) * ")\n";
            end
            tmp = tmp * "\t\t" * opufile * string(k) * "(f, xdg, uinf, param, modelnumber, ng, ncx, nce, npe, ne);\n";
        end
        tmp = tmp * "}\n\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(double *, double *, double *, double *, int, int, int, int, int, int);\n";
        tmp = tmp * "template void " * opufile;
        tmp = tmp * "(float *, float *, float *, float *, int, int, int, int, int, int);\n";
    end

    stropu = tmp;
    ioopu = open(foldername * "/" * opufile * ".cpp", "w");
    write(ioopu, stropu);
    close(ioopu);

    strcpu = replace(stropu, "opu" => "cpu");
    iocpu = open(foldername * "/" * cpufile * ".cpp", "w");
    write(iocpu, strcpu);
    close(iocpu);

    strgpu = replace(stropu, "opu" => "gpu");
    strgpu = replace(strgpu, "cpp" => "cu");
    iogpu = open(foldername * "/" * gpufile * ".cu", "w");
    write(iogpu, strgpu);
    close(iogpu);
end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
using SimJulia
using Distributions
using HypothesisTests
using Printf
using ResumableFunctions
using Random

#set inputs
runs =100
max_customer =150
seed = 123
theta_arrival = 7
theta_service = 5.0
servers = 1
simulation_time = 500
#initialize arrays
arrival_time_Array=Float64[]
departure_time_Array=Float64[]
service_time_Array=Float64[]
time_in_system_Array=Float64[]
time_in_queue_Array=Float64[]
aggregate_time_in_service_Array=Float64[]
aggregate_time_in_queue_Array=Float64[]
aggregate_time_in_system_Array=Float64[]
customers_served=Int[]

#simulation logic
function get_arrival_time()
    inter_arrival = Exponential(theta_arrival)
    push!(arrival_time_Array,rand(inter_arrival))
    return sum(arrival_time_Array)
end

@resumable function queue(env::Environment,name::Int, server::Resource,theta_service::Float64)
    @yield timeout(env, get_arrival_time())
    arrive = now(env)
    #println(name, " arrive at ", arrive)
    @yield request(server)
    served = now(env)
    #println(name, " served at ",served)
    service_time = rand(Exponential(theta_service))
    push!(service_time_Array,service_time)
    @yield timeout(env, service_time)
    depart=now(env)
    push!(departure_time_Array,depart)
    #println(name, " leaving system at ",depart)
    time_in_system=depart-arrive
    push!(time_in_system_Array,time_in_system)
    time_in_queue= time_in_system-service_time
    push!(time_in_queue_Array,time_in_queue)
    @yield release(server)
end
#
for i in 1:runs
    Random.seed!(seed)
    arrival_time_Array=Float64[]
    departure_time_Array=Float64[]
    service_time_Array=Float64[]
    time_in_system_Array=Float64[]
    time_in_queue_Array=Float64[]
    queue
    sim= Simulation()
    server = Resource(sim,servers)
    for i in 1:max_customer
        @process queue(sim, i, server, theta_service)
    end
    run(sim,simulation_time)
    seed=seed+1
#create aggregates

    mean_time_in_service=mean(service_time_Array)
    mean_time_in_system=mean(time_in_system_Array)
    mean_time_in_queue=mean(time_in_queue_Array)
    push!(customers_served,length(departure_time_Array))
    #println("Mean time in service ",mean_time_in_service)
    #println("Mean time in system ",mean_time_in_system)
    #println("Mean time in queue ", mean_time_in_queue)
    #println("Customers served ", customers_served)
    push!(aggregate_time_in_queue_Array,mean_time_in_queue)
    push!(aggregate_time_in_service_Array,mean_time_in_service)
    push!(aggregate_time_in_system_Array,mean_time_in_system)
end
println("Runs: ", runs, "  Sim time: ", simulation_time, "  Seed: ", seed)
println("CAUTION max customers is: ", max_customer, "  ensure it is well over customers served\n")

println("Mean time in service ", mean(aggregate_time_in_service_Array))
println("95% confidecn interval",ci(OneSampleTTest(aggregate_time_in_service_Array)),"\n")

println("Mean time in system ", mean(aggregate_time_in_system_Array))
println("95% confidecn interval",ci(OneSampleTTest(aggregate_time_in_system_Array)),"\n")

println("Mean time in queue ", mean(aggregate_time_in_queue_Array))
println("95% confidecn interval",ci(OneSampleTTest(aggregate_time_in_queue_Array)),"\n")

println("Mean Customers served ", mean(customers_served))
println("95% confidecn interval",ci(OneSampleTTest(customers_served)),"\n")

max_stars_repo_name: jackson-perry/julia-projects

max_issues_repo_name: jackson-perry/julia-projects

max_forks_repo_name: jackson-perry/julia-projects


text:
# replstr(x) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80)), MIME("text/plain"), x), x)
# showstr(x) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80)), x), x)

@testset "Bias correction" begin

d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
Random.seed!(42)
data = randn(2, 2, 10957)
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
obs = ClimateTools.ClimGrid(axisdata, variable = "tasmax", dimension_dict=dimension_dict)
ref = obs * 1.05
fut = obs * 1.05

D = qqmap(obs, ref, fut, method = "Additive", detrend=false)
@test round(D[1][1, 1, 1], digits=5) == -0.55603#68761463861
D = qqmap(obs, ref, fut, method = "Additive", detrend=true)
@test round(D[1][1, 1, 1], digits=2) == -0.56#5603#74620422795

# ===================================
# Ensure we have similar statistics
# ===================================
d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
Random.seed!(42)
data = randn(2, 2, 10957)
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
obs = ClimateTools.ClimGrid(axisdata, variable = "tasmax", dimension_dict=dimension_dict)
ref = obs + 2.0
fut = obs + 2.0

D = qqmap(obs, ref, fut, method="Additive", detrend=false)
@test mean(obs[1][1,1,:]) - mean(fut[1][1,1,:]) ≈ -1.99999999999
@test mean(obs[1][1,1,:]) - mean(D[1][1,1,:]) ≈ -0.00029395108

ref = obs * 2.0
fut = obs * 2.0
D = qqmap(obs, ref, fut, method="Multiplicative", detrend=false)
@test std(obs[1][1,1,:]) - std(fut[1][1,1,:]) ≈ -1.012096874
@test std(obs[1][1,1,:]) - std(D[1][1,1,:]) ≈ -5.3838558178e-5

# ============================
Random.seed!(42)
data = randn(2, 2, 10957)
axisdata = AxisArray(data .- minimum(data), Axis{:lon}(1:2), Axis{:lat}(1:2),Axis{:time}(d))
obs = ClimateTools.ClimGrid(axisdata, variable = "tasmax", dimension_dict=dimension_dict)
ref = obs * 1.05
fut = obs * 1.05

D = qqmap(obs, ref, fut, method = "Multiplicative", detrend=false)
@test round(D[1][1, 1, 1], digits=5)== 3.78144#13533272675

# Create a ClimGrid with a clear trend
x = 1:10957
x = x
d = DateTime(1961,1,1):Day(1):DateTime(1990,12,31)
data = zeros(2, 2, 10957)
data[1,1,:] .= 0.0 .+ x
data[1,2,:] .= 0.0 .+ 0.2.*x .+ x.^2
data[2,1,:] .= 0.0 .+ 0.04.*x - 0.01 .* x.^3
data[2,2,:] .= 0.0 .+ 0.14.*x.^2 .- 0.032 .* x .^4
axisdata = AxisArray(data, Axis{:lon}(1:2), Axis{:lat}(1:2), Axis{:time}(d))
C = ClimGrid(axisdata)
poly = polyfit(C)
val = polyval(C, poly)
D = C - val
@test D[1] == (C - val)[1]

end


text:
@testset "init" begin
    @testset "_print_num_threads_warning" begin
        withenv("SUPPRESS_GAIUS_WARNING" => "false") do
            @test_logs (:warn, "The system has 16 cores. However, Julia was started with 1 thread. We recommend starting Julia with at least 16 threads to take advantage of Gaius's multithreading algorithms. To suppress this warning, set the environment variable SUPPRESS_GAIUS_WARNING=true") match_mode=:any Gaius._print_num_threads_warning(16, 1)
        end
    end

    @testset "_string_to_bool" begin
        @test !Gaius._string_to_bool("")
        @test !Gaius._string_to_bool("foobarbaz")
        @test !Gaius._string_to_bool("false")
        @test !Gaius._string_to_bool("0")
        @test Gaius._string_to_bool("true")
        @test Gaius._string_to_bool("1")
    end

    @testset "_is_suppress_warning" begin
        withenv("SUPPRESS_GAIUS_WARNING" => nothing) do
            @test !Gaius._is_suppress_warning()
        end
        withenv("SUPPRESS_GAIUS_WARNING" => "") do
            @test !Gaius._is_suppress_warning()
        end
        withenv("SUPPRESS_GAIUS_WARNING" => "false") do
            @test !Gaius._is_suppress_warning()
        end
        withenv("SUPPRESS_GAIUS_WARNING" => "0") do
            @test !Gaius._is_suppress_warning()
        end
        withenv("SUPPRESS_GAIUS_WARNING" => "true") do
            @test Gaius._is_suppress_warning()
        end
        withenv("SUPPRESS_GAIUS_WARNING" => "1") do
            @test Gaius._is_suppress_warning()
        end
    end

    @testset "_pluralize" begin
        @test Gaius._pluralize("foo", "bar", 0) == "bar"
        @test Gaius._pluralize("foo", "bar", 1) == "foo"
        @test Gaius._pluralize("foo", "bar", 2) == "bar"
        @test Gaius._pluralize("foo", "bar", 3) == "bar"
    end

    @testset "_pluralize_threads" begin
        @test Gaius._pluralize_threads(0) == "0 threads"
        @test Gaius._pluralize_threads(1) == "1 thread"
        @test Gaius._pluralize_threads(2) == "2 threads"
        @test Gaius._pluralize_threads(3) == "3 threads"
    end

    @testset "_pluralize_cores" begin
        @test Gaius._pluralize_cores(0) == "0 cores"
        @test Gaius._pluralize_cores(1) == "1 core"
        @test Gaius._pluralize_cores(2) == "2 cores"
        @test Gaius._pluralize_cores(3) == "3 cores"
    end
end


text:
# Test of Dijkstra's algorithm for shorest paths

using Graphs
using Base.Test

# g1: the example in CLRS (2nd Ed.)
g1 = simple_inclist(5)

g1_wedges = [
    (1, 2, 10.),
    (1, 3, 5.),
    (2, 3, 2.),
    (3, 2, 3.),
    (2, 4, 1.),
    (3, 5, 2.),
    (4, 5, 4.),
    (5, 4, 6.),
    (5, 1, 7.),
    (3, 4, 9.) ]

ne = length(g1_wedges)
eweights1 = zeros(ne)
for i = 1 : ne
    we = g1_wedges[i]
    add_edge!(g1, we[1], we[2])
    eweights1[i] = we[3]
end

@assert num_vertices(g1) == 5
@assert num_edges(g1) == 10

s1 = dijkstra_shortest_paths(g1, eweights1, 1)

@test s1.parents == [1, 3, 1, 2, 3]
@test s1.dists == [0., 8., 5., 9., 7.]
@test s1.colormap == [2, 2, 2, 2, 2]

s1 = dijkstra_shortest_paths(g1, eweights1, [1])

@test s1.parents == [1, 3, 1, 2, 3]
@test s1.dists == [0., 8., 5., 9., 7.]
@test s1.colormap == [2, 2, 2, 2, 2]

g1ex = graph([i for i in 1:5], ExEdge{Int}[], is_directed=true)

#todo -- this should be doable in a comprehension, I can't get the types to work
for (i,v) in enumerate(g1_wedges)
    ed = ExEdge(i, v[1], v[2])
    ed.attributes["length"] = v[3]
    add_edge!(g1ex, ed)
end

edgel = AttributeEdgePropertyInspector{Float64}("length")
s1ex = dijkstra_shortest_paths(g1ex, edgel, [1])

@test s1ex.parents == [1, 3, 1, 2, 3]
@test s1ex.dists == [0., 8., 5., 9., 7.]
@test s1ex.colormap == [2, 2, 2, 2, 2]

# Check early termination

type EndWhenNode <: AbstractDijkstraVisitor
  n::Int
end

function Graphs.include_vertex!(visitor::EndWhenNode, u, v, d)
  v != visitor.n
end

s1b = dijkstra_shortest_paths(g1, eweights1, [1], visitor=EndWhenNode(5))

@test s1b.parents == [1, 3, 1, 3, 3]
@test s1b.dists == [0.0, 8.0, 5.0, 14.0, 7.0]
@test s1b.colormap == [2, 1, 2, 1, 2]

# g2: the example in Wikipedia
g2 = simple_inclist(6, is_directed=false)

g2_wedges = [
    (5, 6, 9.),
    (5, 4, 6.),
    (6, 3, 2.),
    (4, 3, 11.),
    (6, 1, 14.),
    (3, 1, 9.),
    (3, 2, 10.),
    (4, 2, 15.),
    (1, 2, 7.) ]

ne = length(g2_wedges)
eweights2 = zeros(ne)
for i = 1 : ne
    we = g2_wedges[i]
    add_edge!(g2, we[1], we[2])
    eweights2[i] = we[3]
end

@assert num_vertices(g2) == 6
@assert num_edges(g2) == 9

s2 = dijkstra_shortest_paths(g2, eweights2, 1)

@test s2.parents == [1, 1, 1, 3, 6, 3]
@test s2.dists == [0., 7., 9., 20., 20., 11.]
@test s2.colormap == [2, 2, 2, 2, 2, 2]

s2 = dijkstra_shortest_paths(g2, eweights2, [1])

@test s2.parents == [1, 1, 1, 3, 6, 3]
@test s2.dists == [0., 7., 9., 20., 20., 11.]
@test s2.colormap == [2, 2, 2, 2, 2, 2]

g3 = simple_graph(4)
add_edge!(g3,1,2); add_edge!(g3,1,3); add_edge!(g3,2,3); add_edge!(g3,3,4); add_edge!(g3,4,3); add_edge!(g3,3,1)

s3 = dijkstra_shortest_paths(g3,2)
sps = enumerate_paths(vertices(g3), s3.parent_indices)
@test length(sps) == 4
@test sps[1] == [2,3,1]
@test sps[2] == [2]
@test sps[3] == [2, 3]
@test sps[4] == [2, 3, 4]

sps = enumerate_paths(vertices(g3), s3.parent_indices, [2,4])
@test length(sps) == 2
@test sps[1] == [2]
@test sps[2] == [2, 3, 4]

sps = enumerate_paths(vertices(g3), s3.parent_indices, 4)
@test sps == [2, 3, 4]

g4 = Graphs.inclist([4,5,6,7],is_directed=true)
add_edge!(g4,4,5); add_edge!(g4,4,6); add_edge!(g4,5,6); add_edge!(g4,6,7)

s4 = dijkstra_shortest_paths(g4,5)
sps = enumerate_indices(s4.parent_indices)
@test length(sps) == 4
@test sps[1] == []
@test sps[2] == [2]
@test sps[3] == [2, 3]
@test sps[4] == [2, 3, 4]

sps = enumerate_indices(s4.parent_indices, [2,4])
@test length(sps) == 2
@test sps[1] == [2]
@test sps[2] == [2, 3, 4]

sps = enumerate_indices(s4.parent_indices, 4)
@test sps == [2, 3, 4]

sps = enumerate_paths(vertices(g4), s4.parent_indices)
@test length(sps) == 4
@test sps[1] == []
@test sps[2] == [5]
@test sps[3] == [5, 6]
@test sps[4] == [5, 6, 7]

sps = enumerate_paths(vertices(g4), s4.parent_indices, [2,4])
@test length(sps) == 2
@test sps[1] == [5]
@test sps[2] == [5, 6, 7]

sps = enumerate_paths(vertices(g4), s4.parent_indices, 4)
@test sps == [5, 6, 7]


text: abstract type CameraBase end


text:
"""
    Functions in this file read the subfind output.

"""

struct SubfindHeader
    nhalos::Int32                       # number of halos in the output file
    nsubhalos::Int32                    # number of subhalos in the output file
    nfof::Int32                         # number of particles in the FoF
    ngroups::Int32                      # number of large groups in the output file
    time::Float64                       # time / scale factor of the simulation
    z::Float64                          # redshift of the simulation
    tothalos::UInt32                    # total number of halos over all output files
    totsubhalos::UInt32                 # total number of subhalos over all output files
    totfof::UInt32                      # total number of particles in the FoF
    totgroups::UInt32                   # total number of large groups over all output files
    num_colors::Int32                   # number of colors
    boxsize::Float64                    # total size of the simulation box
    omega_0::Float64                    # Omega matter
    omega_l::Float64                    # Omega dark enery
    h0::Float64                         # little h
    flag_doubleprecision::Int32         # 1 if snapshot is in double precision, else 0
    flag_ic_info::Int32
end

function read_subfind_header(filename::String)

    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize == 8
        swap = 0
        snap_format = 2
    elseif blocksize == 256
        swap = 0
        snap_format = 1
    else
        blocksize = bswap(blocksize)
        if blocksize == 8
            swap = 1
            snap_format = 2
        elseif blocksize == 256
            swap = 1
            snap_format = 1
        else
            println("incorrect file format encountered when reading header of", filename)
        end
    end

    #println("Reading snapshot format: ", snap_format)

    if snap_format == 2
        seek(f, 16)
        skip_line = read(f, Int32)
    end

    nhalos = read(f, Int32)
    nsubhalos = read(f, Int32)
    nfof = read(f, Int32)
    ngroups = read(f, Int32)

    # rest of numbers and massarr are obsolete
    DUMMY = read!(f, Array{Int32,1}(undef,2))
    DUMMY = read!(f, Array{Float64,1}(undef,6))

    time = read(f, Float64)
    z = read(f, Float64)

    # flags not needed
    DUMMY = read(f, Int32)
    DUMMY = read(f, Int32)

    tothalos = read(f, UInt32)
    totsubhalos = read(f, UInt32)
    totfof = read(f, UInt32)
    totgroups = read(f, UInt32)

    DUMMY = read!(f, Array{Int32,1}(undef,2))
    DUMMY = read(f, Int32)

    num_files = read(f, Int32)
    boxsize = read(f, Float64)
    omega_0 = read(f, Float64)
    omega_l = read(f, Float64)
    h0 = read(f, Float64)

    DUMMY = read(f, Int32)
    DUMMY = read(f, Int32)
    DUMMY = read!(f, Array{UInt32,1}(undef,6))
    DUMMY = read(f, Int32)

    flag_doubleprecision = read(f, Int32)
    flag_ic_info = read(f, Int32)

    close(f)

    return SubfindHeader(nhalos, nsubhalos, nfof, ngroups,
                         time, z,
                         tothalos, totsubhalos, totfof, totgroups,
                         num_files,
                         boxsize, omega_0, omega_l, h0,
                         flag_doubleprecision,
                         flag_ic_info)
end

function read_subfind(filename::String, blockname::String)

    info = read_info(filename)

    info = info[getfield.(info, :block_name) .== blockname][1]

    parttype = findall(info.is_present .== 1)[1] - 1

    return read_snap(filename, blockname, parttype)
end

using ProgressMeter
using Base.Threads

struct SubfindID
    file::Int64
    id::Int64
end

@inline function select_file(filebase::String, i::Int, nfiles::Int)

    if nfiles > 1
        sub_input = filebase * ".$i"
        if !isfile(sub_input)
            error("File $sub_input not present!")
        end
    else  # nfiles > 1
        sub_input = filebase

        if !isfile(sub_input)
            sub_input = filebase * ".0"
        end

        if !isfile(sub_input)
            error("Subfind file not present!")
        end

    end  # nfiles > 1

    return sub_input
end

"""
    find_most_massive_halo(filebase::String [, nfiles::Int=1])

Reads the selected file and its subfiles and returns position, virial radius
and a SubfindID object that contains the subfile which contains the most massive
halo and the position in the block.

"""
function find_most_massive_halo(filebase::String, nfiles::Int=1)

    Mmax     = 0.0
    POS      = zeros(Float32, 3)
    RVIR     = 0.0
    max_file = 0
    max_id   = 0

    @showprogress "Reading files..." for i = 0:nfiles-1

        sub_input = select_file(filebase, i, nfiles)

        M = read_subfind(sub_input, "MVIR")
        max_test = findmax(M)

        if max_test[1] > Mmax

            max_file = i
            max_id   = max_test[2][1]

            # store position and virial radius of most massive halo
            POS  = read_subfind(sub_input, "GPOS")[max_id,:]
            RVIR = read_subfind(sub_input, "RVIR")[max_id]

            # store new maximum mass
            Mmax = max_test[1]
        end

    end # for

    return POS, RVIR, SubfindID(max_file, max_id)
end

struct SubfindFilter
    file::Int64
    index::Array{Int64}
end

"""
    filter_subfind(filebase::String, blockname::String, filter_function [, nfiles::Integer=1])

Selects entries in subfind block that fulfill the 'filter_funcion' requirements and
returns a 'SubfindFilter' object.

# Examples
```jldoctest
julia> find_mass_gt_1e15(M) = ( (M > 1.e15) ? true : false )
find_mass_gt_1e15 (generic function with 1 method)
julia> filtered_subfind = filter_subfind(filebase, "MVIR", find_mass_gt_1e15)
[...]
```

"""
function filter_subfind(filebase::String, blockname::String, filter_function, nfiles::Integer=1)

    # allocate empty array for SubfindFilter objects
    A = Vector{SubfindFilter}(undef, 0)

    # loop over all files in parallel
    @threads for i = 0:nfiles-1

        sub_input = select_file(filebase, i, nfiles)

        block = read_subfind(sub_input, blockname)

        selection = filter_function.(block)
        correct_selection = findall(selection)

        if length(selection[selection]) > 0

            # create array of integers for easy storing
            id_array = collect(1:length(selection))[correct_selection]

            # create SubfindFilter object and push it to the storage array
            push!(A, SubfindFilter(i, id_array))

        end # if entries > 0

    end # for loop

    return A
end

function read_particles_in_halo(sub_file::String, snap_filebase::String,
                                halo_id::Int64, blocks::Vector{String};
                                parttype::Integer=0, Rhalf_factor::AbstractFloat=5.0)

    error("Not finished!")

    # read subfind data

    pos      = read_subfind(sub_input, "SPOS")[halo_id,:]
    rhalf    = read_subfind(sub_input, "RHMS")[halo_id]

    sub_info = read_info(sub_file)
    pos = get_block_positions(sub_file)
    ids = read_block_with_offset(sub_info, "KEY",
                                      info = sub_info[getfield.(sub_info, :block_name) .== "PID"][1],
                                      parttype = 0)
end

"""
            Sandbox
"""

"""
using GadJet
import GadJet.Info_Line
import GadJet.read_info_line
import GadJet.read_info
import GadJet.read_particles_in_volume
import GadJet.read_subfind

snap_filebase = "/HydroSims/Magneticum/Box2/hr_bao/snapdir_140/snap_140"
sub_input  = "/HydroSims/Magneticum/Box2/hr_bao/groups_140/sub_140.50"
halo_id    = UInt64(1)
pos        = read_subfind(sub_input, "SPOS")[halo_id,:]
rhalf      = read_subfind(sub_input, "RHMS")[halo_id]
ids_offset = read_subfind(sub_input, "SOFF")[halo_id]
n_ids      = read_subfind(sub_input, "SLEN")[halo_id]

Int64(n_ids)

tot_ids      = read_subfind(sub_input, "PID")

int_ids    = Int64.(tot_ids)

sub_info = read_info(sub_input)

ids_info = sub_info[getfield.(sub_info, :block_name) .== "PID"][1]

ids_info.data_type = UInt64

ids_pos = get_block_positions(sub_input)["PID"]
ids = Array{ids_info.data_type, 1}(undef, 0)

tot_offset = ids_pos# + ids_offset
ids = read_block_with_offset(sub_input, ids, ids_pos, ids_info, ids_offset, [0], n_ids, [n_ids])

typeof(pos) <: DenseArray{Float32,1}

supertype(supertype(supertype(typeof(pos))))

read_particles_in_volume(snap_filebase, "ID", pos, 5rhalf)

2^63 - 1

a = 1:100

b = [3, 1, 17, 22, 12]

selected = [a .== b[i] for i = 1:length(b)]

c = reduce(+, selected)
"""

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
using UnsteadyFlowSolvers

alpha_init = 10. *pi/180
alphadot_init = 0.
h_init = 0.
hdot_init = 0.
u = 0.467
udot = 0
kinem = KinemPar2DOF(alpha_init, h_init, alphadot_init, hdot_init, u)

x_alpha = 0.05
r_alpha = 0.5
kappa = 0.05
w_alpha = 1.
w_h = 1.
w_alphadot = 0.
w_hdot = 0.
cubic_h_1 = 1.
cubic_h_3 = 0.
cubic_alpha_1 = 1.
cubic_alpha_3 = 0.
strpar = TwoDOFPar(x_alpha, r_alpha, kappa, w_alpha, w_h, w_alphadot, w_hdot, cubic_h_1, cubic_h_3, cubic_alpha_1, cubic_alpha_3)

#Kinematic definitions to initialise surface (initial conditions)
alphadef = ConstDef(alpha_init)
hdef = ConstDef(h_init)
udef = ConstDef(1.) #This is relative to uref
startkinem = KinemDef(alphadef, hdef, udef)

lespcrit = [0.11;]

pvt = 0.35

c = 1.

surf = TwoDSurf("FlatPlate", pvt, startkinem, lespcrit, c=c, uref=u)

curfield = TwoDFlowField()

dtstar = 0.015

nsteps = 500

t_tot = nsteps * dtstar / u

startflag = 0

writeflag = 0
writeInterval = dtstar * nsteps/20.

writeInterval = t_tot/20.

delvort = delSpalart(500, 12, 1e-5)

mat, surf, curfield = ldvm2DOF(surf, curfield, strpar, kinem, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

makeForcePlots2D()

#cleanWrite()

max_stars_repo_name: MaximeUofG/UNSflowTutorials

max_issues_repo_name: MaximeUofG/UNSflowTutorials

max_forks_repo_name: MaximeUofG/UNSflowTutorials


text:
# General API

inds(A::AbstractArray, d) = 1:size(A, d)
inds(A::AbstractArray{T,N}) where {T,N} = ntuple(d->inds(A,d), Val(N))

eachindex(x::ArrayIndexingWrapper) = each(index(x))

function show(io::IO, W::ArrayIndexingWrapper)
    print(io, "iteration hint over ", hint_string(W), " of a ", summary(W.data), " over the region ", W.indexes)
end

function show(io::IO, iter::ContigCartIterator)
    print(io, "Cartesian iterator with:\n  domain ", iter.arrayrange, "\n  start  ", first(iter.columnrange), "\n  stop   ", last(iter.columnrange))
end

parent(W::ArrayIndexingWrapper) = W.data
parent(F::Follower) = F.value

stripwrapper(A::AbstractArray) = A
stripwrapper(A::ArrayIndexingWrapper) = parent(A)

"""
`index(A)`
`index(A, indexes...)`

`index` creates an "iteration hint" that records the region of `A`
that you wish to iterate over. The iterator will return the indexes,
rather than values, of `A`. "iteration hints" are not iterables; to
create an iterator from a hint, call `each` on the resulting object.

In contrast to `eachindex` iteration over a subarray of `A`, the
indexes are for `A` itself.

See also: `value`, `stored`, `each`.
"""
index(w::ArrayIndexingWrapper{A,I,isindex,isstored}) where {A,I,isindex,isstored} = ArrayIndexingWrapper{A,I,true,isstored}(w.data, w.indexes)

"""
`value(A)`
`value(A, indexes...)`

`value` creates an "iteration hint" that records the region of `A`
that you wish to iterate over. The iterator will return the values,
rather than the indexes, of `A`. "iteration hints" are not iterables; to
create an iterator from a hint, call `each` on the resulting object.

See also: `index`, `stored`, `each`.
"""
value(w::ArrayIndexingWrapper{A,I,isindex,isstored}) where {A,I,isindex,isstored} = ArrayIndexingWrapper{A,I,true,isstored}(w.data, w.indexes)

"""
`stored(A)`
`stored(A, indexes...)`

`stored` creates an "iteration hint" that records the region of `A`
that you wish to iterate over. The iterator will return just the
stored values of `A`. "iteration hints" are not iterables; to create
an iterator from a hint, call `each` on the resulting object.

See also: `index`, `value`, `each`.
"""
stored(w::ArrayIndexingWrapper{A,I,isindex,isstored}) where {A,I,isindex,isstored} = ArrayIndexingWrapper{A,I,isindex,true}(w.data, w.indexes)

allindexes(A::AbstractArray{T,N}) where {T,N} = ntuple(d->Colon(),Val(N))

index(A::AbstractArray) = index(A, allindexes(A))
index(A::AbstractArray, I::IterIndex...) = index(A, I)
index(A::AbstractArray{T,N}, indexes::NTuple{N,IterIndex}) where {T,N} = ArrayIndexingWrapper{typeof(A),typeof(indexes),true,false}(A, indexes)

value(A::AbstractArray) = value(A, allindexes(A))
value(A::AbstractArray, I::IterIndex...) = value(A, I)
value(A::AbstractArray{T,N}, indexes::NTuple{N,IterIndex}) where {T,N} = ArrayIndexingWrapper{typeof(A),typeof(indexes),false,false}(A, indexes)

stored(A::AbstractArray) = stored(A, allindexes(A))
stored(A::AbstractArray, I::IterIndex...) = stored(A, I)
stored(A::AbstractArray{T,N}, indexes::NTuple{N,IterIndex}) where {T,N} = ArrayIndexingWrapper{typeof(A),typeof(indexes),false,true}(A, indexes)

"""
`each(iterhint)`
`each(iterhint, indexes...)`

`each` instantiates the iterator associated with `iterhint`. In
conjunction with `index` and `stored`, you may choose to iterate over
either indexes or values, as well as choosing whether to iterate over
all elements or just the stored elements.
"""
each(A::AbstractArray) = each(A, allindexes(A))
each(A::AbstractArray, indexes::IterIndex...) = each(A, indexes)
each(A::AbstractArray{T,N}, indexes::NTuple{N,IterIndex}) where {T,N} = each(ArrayIndexingWrapper{typeof(A),typeof(indexes),false,false}(A, indexes))

# Fallback definitions for each
each(W::ArrayIndexingWrapper{A,I,false,isstored}) where {A,I,isstored} = (itr = each(index(W)); ValueIterator{A,typeof(itr)}(W.data, itr))
each(W::ArrayIndexingWrapper{A,NTuple{N,<:CSlice},true,isstored}) where {A,N,isstored} = eachindex(W.data)
each(W::ArrayIndexingWrapper{A,I,true,isstored}) where {A,I,isstored} = _each(contiguous_index(W.indexes), W)

_each(::Contiguous, W) = contiguous_iterator(W)
_each(::Any, W) = CartesianIndices(ranges(W))

function iterate(vi::ValueIterator, s...)
    v = iterate(vi.iter, s...)
    v === nothing && return v
    idx, s = v
    return vi.data[idx], s
end

function iterate(iter::SyncedIterator, state...)
    v = iterate(iter.iter, state...)
    v === nothing && return v
    item, newstate = v
    mapf(iter.itemfuns, iter.items, item), newstate
end

function iterate(itr::FirstToLastIterator, i=(itr.itr,))
    v = iterate(i...)
    v === nothing && return v
    idx, s = v
    itr.parent[idx], (i[1], s)
end

# SyncedIterator(iter, funcs) = SyncedIterator{typeof(iter), Base.typesof(funcs)}(iter, funcs)

function sync(A::AllElements, B::AllElements)
    checksame_inds(A, B)
    _sync(checksame_storageorder(A, B), A, B)
end

function sync(A::AllElements, B::AllElements...)
    checksame_inds(A, B...)
    _sync(checksame_storageorder(A, B...), A, B...)
end

_sync(::Type{Val{true}}, A, B) = zip(each(A), each(B))
_sync(::Type{Val{false}}, A, B) = zip(columnmajoriterator(A), columnmajoriterator(B))
_sync(::Type{Val{true}}, As...) = zip(map(each, As)...)
_sync(::Type{Val{false}}, As...) = zip(map(columnmajoriterator, As)...)

# For stored, see sync_stored.jl

### Utility methods

"""
`mapf(fs, objs, x)` is similar to `map(f, a, b)`, except instead of mapping one
function over many objects, it maps many function/object pairs over one
`x`. `fs` should be a tuple-of-functions, and `objs` a tuple-of-containers.
"""
mapf(fs::Tuple{Vararg{<:Any,N}}, objs::Tuple{Vararg{<:Any,N}}, x) where N = _mapf((), fs, objs, x)
_mapf(out, ::Tuple{}, ::Tuple{}, x) = out
@inline function _mapf(out, fs, objs, x)
    f, obj = fs[1], objs[1]
    ret = _mapf((out..., f(obj, x)), Base.tail(fs), Base.tail(objs), x)
end

storageorder(::Array) = FirstToLast()
storageorder(::PermutedDimsArray{T,N,AA,perm}) where {T,N,AA,perm} = OtherOrder{perm}()
storageorder(A::ReshapedArray) = _so(storageorder(parent(A)))
storageorder(A::AbstractArray) = storageorder(parent(A)) # parent required!

storageorder(W::ArrayIndexingWrapper) = storageorder(parent(W))

_so(o::FirstToLast) = o
_so(::Any) = NoOrder() # reshape + permutedims => undefined

hint_string(::ArrayIndexingWrapper{A,I,false,false}) where {A,I} = "values"
hint_string(::ArrayIndexingWrapper{A,I,true,false}) where {A,I} = "indexes"
hint_string(::ArrayIndexingWrapper{A,I,false,true}) where {A,I} = "stored values"
hint_string(::ArrayIndexingWrapper{A,I,true,true}) where {A,I} = "indexes of stored values"

ranges(W) = ranges((), W.data, 1, W.indexes...)
ranges(out, A, d) = out
@inline ranges(out, A, d, i, I...) = ranges((out..., i), A, d+1, I...)
@inline ranges(out, A, d, i::CSlice, I...) = ranges((out..., inds(A, d)), A, d+1, I...)

checksame_inds(::Type{Bool}, A::ArrayOrWrapper) = true
checksame_inds(::Type{Bool}, A::ArrayOrWrapper, B::ArrayOrWrapper) = extent_inds(A) == extent_inds(B)
checksame_inds(::Type{Bool}, A, B, C...) = checksame_inds(Bool, A, B) && checksame_inds(Bool, B, C...)
checksame_inds(A) = checksame_inds(Bool, A)
checksame_inds(A, B) = checksame_inds(Bool, A, B) || throw(DimensionMismatch("extent inds $(extent_inds(A)) and $(extent_inds(B)) do not match"))
checksame_inds(A, B, C...) = checksame_inds(A, B) && checksame_inds(B, C...)

# extent_inds drops sliced dimensions
extent_inds(A::AbstractArray) = inds(A)
extent_inds(W::ArrayIndexingWrapper) = _extent_inds((), W.data, 1, W.indexes...)
_extent_inds(out, A, d) = out
@inline _extent_inds(out, A, d, ::Int, indexes...) = _extent_inds(out, A, d+1, indexes...)
@inline _extent_inds(out, A, d, i, indexes...) = _extent_inds((out..., inds(A, d)), A, d+1, indexes...)

# extent_dims indicates which dimensions have extended size
extent_dims(A::AbstractArray{T,N}) where {T,N} = ntuple(identity,Val(N))
extent_dims(W::ArrayIndexingWrapper) = _extent_dims((), 1, W.indexes...)
_extent_dims(out, d::Integer) = out
@inline _extent_dims(out, d, i1::Union{UnitRange{Int},CSlice}, indexes...) = _extent_dims((out..., d), d+1, indexes...)
@inline _extent_dims(out, d, i1, indexes...) = _extent_dims(out, d+1, indexes...)

columnmajoriterator(A::AbstractArray) = columnmajoriterator(IndexStyle(A), A)
columnmajoriterator(::IndexLinear, A) = A
columnmajoriterator(::IndexCartesian, A) = FirstToLastIterator(A, CartesianIndices(size(A)))

columnmajoriterator(W::ArrayIndexingWrapper) = CartesianIndices(ranges(W))

checksame_storageorder(A) = Val{true}
checksame_storageorder(A, B) = _sso(storageorder(A), storageorder(B))
checksame_storageorder(A, B, C...) = checksame_storageorder(_sso(storageorder(A), storageorder(B)), B, C...)
checksame_storageorder(::Type{Val{true}}, A, B...) = checksame_storageorder(A, B...)
checksame_storageorder(::Type{Val{false}}, A, B...) = Val{false}
_sso(::FirstToLast, ::FirstToLast) = Val{true}
_sso(::OtherOrder{p}, ::OtherOrder{p}) where {p} = Val{true}
_sso(::StorageOrder, ::StorageOrder) = Val{false}

# indexes is contiguous if it's one of:
#    Colon...
#    Colon..., Union{UnitRange,Int}, Int...
@inline contiguous_index(I) = contiguous_index(Contiguous(), I...)
@inline contiguous_index(c::Contiguous, ::CSlice, I...) = contiguous_index(c, I...)
@inline contiguous_index(::Contiguous, ::Any, I...) = contiguous_index(MaybeContiguous(), I...)
@inline contiguous_index(c::MaybeContiguous, ::Int, I...) = contiguous_index(c, I...)
@inline contiguous_index(::MaybeContiguous, ::Any, I...) = NonContiguous()
@inline contiguous_index(::Contiguity) = Contiguous()  # won't get here for NonContiguous

contiguous_iterator(W) = _contiguous_iterator(W, IndexStyle(parent(W)))
function _contiguous_iterator(W, ::IndexLinear)
    f, l = firstlast(W)
    f:l
end
_contiguous_iterator(W, ::IndexCartesian) = CartesianIndices(ranges(W))

# Return the "corners" of an iteration range
function firstlast(W)
    A = parent(W)
    f = firstlast(first, A, W.indexes)
    l = firstlast(last, A, W.indexes)
    f, l
end

# This effectively implements 2-argument map, but without allocating
# intermediate tuples
sub2ind = Base._sub2ind
@inline firstlast(f, A, indexes) = sub2ind(size(A), _firstlast((), f, A, indexes...)...)
@inline _firstlast(out, f, A) = out
@inline function _firstlast(out, f, A, i1, indexes...)
    d = length(out)+1
    _firstlast((out..., f(inds(A, d)[i1])), f, A, indexes...)
end


text:
# Laplacian Eigenmaps
# -------------------
# Laplacian Eigenmaps for Dimensionality Reduction and Data Representation,
# M. Belkin, P. Niyogi, Neural Computation, June 2003; 15 (6):1373-1396

"""
    LEM{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction

The `LEM` type represents a Laplacian eigenmaps model constructed for `T` type data with a help of the `NN` nearest neighbor algorithm.
"""
struct LEM{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction
    λ::AbstractVector{T}
    ɛ::T
    proj::Projection{T}
    nearestneighbors::NN
    component::AbstractVector{Int}
end

## properties
outdim(R::LEM) = size(R.proj, 1)
eigvals(R::LEM) = R.λ
neighbors(R::LEM) = R.nearestneighbors.k
vertices(R::LEM) = R.component

## show
summary(io::IO, R::LEM) = print(io, "Laplacian Eigenmaps(outdim = $(outdim(R)), neighbors = $(neighbors(R)), ɛ = $(R.ɛ))")

## interface functions
"""
    fit(LEM, data; k=12, maxoutdim=2, ɛ=1.0, nntype=BruteForce)

Fit a Laplacian eigenmaps model to `data`.

# Arguments
* `data`: a matrix of observations. Each column of `data` is an observation.

# Keyword arguments
* `k`: a number of nearest neighbors for construction of local subspace representation
* `maxoutdim`: a dimension of the reduced space.
* `nntype`: a nearest neighbor construction class (derived from `AbstractNearestNeighbors`)
* `ɛ`: a Gaussian kernel variance (the scale parameter)

# Examples
```julia
M = fit(LEM, rand(3,100)) # construct Laplacian eigenmaps model
R = transform(M)          # perform dimensionality reduction
```
"""
function fit(::Type{LEM}, X::AbstractMatrix{T};
        k::Int=12, maxoutdim::Int=2, ɛ::Real=1.0, nntype=BruteForce) where {T<:Real}
    # Construct NN graph
    NN = fit(nntype, X, k)
    D, E = knn(NN, X)
    G, C = largest_component(SimpleWeightedGraph(adjmat(D,E)))

    # Compute weights
    W = weights(G)
    W .^= 2
    W ./= maximum(W)

    W[W .> eps(T)] = exp.(-W[W .> eps(T)] ./ convert(T,ɛ))
    D = spdiagm(0=>vec(sum(W, dims=2)))
    L = D - W

    λ, V = decompose(L, D, maxoutdim)
    return LEM{nntype, T}(λ, ɛ, transpose(V), NN, C)
end

"""
    transform(R::LEM)

Transforms the data fitted to the Laplacian eigenmaps model `R` into a reduced space representation.
"""
transform(R::LEM) = R.proj