text:
"""
    HalfIntegerArrays.OneTo(n)

Define an AbstractUnitRange that behaves like `1:n`, with the added distinction that the lower limit is guaranteed (by the type
system) to be `1`. The elements are integers, but of a `HalfInteger` type.
"""
struct OneTo{T<:HalfInteger} <: AbstractUnitRange{T}
    stop :: T

    OneTo{T}(stop) where {T<:HalfInteger} = new(max(zero(T), stop))
    function OneTo{T}(r::AbstractRange{Q}) where {T<:HalfInteger, Q}
        throwstart(r) = throw(ArgumentError("first element must be 1, got $(first(r))"))
        throwstep(r)  = throw(ArgumentError("step must be 1, got $(step(r))"))
        first(r) == one(Q) || throwstart(r)
        step(r)  == one(Q) || throwstep(r)
        return new(max(zero(T), last(r)))
    end
end
OneTo(stop::T) where {T<:Integer} = OneTo{Half{T}}(stop)
OneTo(stop::T) where {T<:HalfInteger} = OneTo{T}(stop)
OneTo(r::AbstractRange{T}) where {T<:Integer} = OneTo{Half{T}}(r)
OneTo(r::AbstractRange{T}) where {T<:HalfInteger} = OneTo{T}(r)
OneTo(r::OneTo{T}) where {T<:HalfInteger} = r
OneTo{T}(r::OneTo{T}) where {T<:HalfInteger} = r
OneTo{T}(r::Union{OneTo,Base.OneTo}) where {T<:HalfInteger} = OneTo{T}(r.stop)

Base.first(::OneTo{Q}) where {Q} = oneunit(Q)
Base.step(::OneTo{Q}) where {Q} = oneunit(Q)
Base.last(b::OneTo) = b.stop

Base.unsafe_length(r::OneTo) = unsafeInt(r.stop)
Base.length(r::OneTo) = Base.unsafe_length(r)

@inline function Base.getindex(v::OneTo{T}, i::Integer) where T
    @boundscheck ((i > 0) & (i <= v.stop)) || throw_boundserror(v, i)
    convert(T, i)
end
@inline function Base.getindex(r::OneTo{T}, s::Union{OneTo,Base.OneTo}) where T
    @boundscheck checkbounds(r, s)
    OneTo(T(s.stop))
end

# conversions
Base.convert(::Type{AbstractUnitRange{Q}}, b::Base.OneTo) where {Q<:HalfInteger} = OneTo{Q}(b.stop)
Base.convert(::Type{AbstractUnitRange{Q}}, b::OneTo) where {Q<:HalfInteger} = OneTo{Q}(b.stop)
Base.convert(::Type{AbstractUnitRange{Q}}, b::OneTo) where {Q<:Integer} = Base.OneTo{Q}(b)

Base.OneTo{Q}(b::OneTo) where {Q<:Integer} = Base.OneTo{Q}(unsafeInt(b.stop))
Base.OneTo(b::OneTo) = Base.OneTo(unsafeInt(b.stop))

"""
    ro = IdOffsetRange(r::AbstractUnitRange, offset::HalfInteger)

Construct an "identity offset range".
Numerically, `collect(ro) == collect(r) .+ offset`, with the additional property that
`axes(ro,1) = axes(r, 1) .+ offset`.
When `r` starts at `1`, then `ro[i] == i` and even `ro[ro] == ro`, i.e., it's the "identity," which is the
origin of the "Id" in `IdOffsetRange`. The element type is a `HalfInteger`.
"""
struct IdOffsetRange{T<:HalfInteger,I<:AbstractUnitRange{<:Integer}} <: AbstractUnitRange{T}
    parent::I
    offset::T

    function IdOffsetRange{T,I}(r::I, o::T) where {T<:HalfInteger,I<:AbstractUnitRange{<:Integer}}
        new{T,I}(r,o)
    end
end

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

IdOffsetRange(r::AbstractUnitRange{T}) where {T<:Integer} = IdOffsetRange(r, zero(Half{T}))

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

function offset_coerce(::Type{R}, r::AbstractUnitRange{T}) where {T,R<:Union{Base.OneTo,OneTo}}
    o = first(r) - one(eltype(r))
    return R(floor(Int,last(r) - o)), T(o)
end

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

Base.eltype(::IdOffsetRange{T,<:AbstractUnitRange{Q}}) where {T,Q} = promote_type(T,Q)

@inline offset(iter::IdOffsetRange) = iter.offset
@inline Base.parent(r::IdOffsetRange) = r.parent
@inline Base.axes(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.axes1(r::IdOffsetRange) = IdOffsetRange(Base.axes1(parent(r)), offset(r))
@inline Base.unsafe_indices(r::IdOffsetRange) = (r,)
@inline Base.unsafe_indices(r::Base.Slice{<:IdOffsetRange}) = (r.indices,)
@inline Base.length(r::IdOffsetRange) = length(parent(r))
@inline Base.compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{IdOffsetRange}, I::Tuple) =
    Base.compute_linindex(parent, I) - stride1*first(inds[1])
Base.reduced_index(i::IdOffsetRange) = typeof(i)(first(parent(i)):first(parent(i)),offset(i))
Base.reduced_index(i::IdentityUnitRange{<:IdOffsetRange}) = Base.reduced_index(i.indices)

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

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

function unsafeUnitRangeInt(r::AbstractUnitRange{HalfInt})
    unsafeInt(first(r)):unsafeInt(last(r))
end

@propagate_inbounds function Base.getindex(r::IdOffsetRange{HalfInt}, i::HalfInt)
    r.parent[unsafeInt(i - r.offset)] + r.offset
end
@propagate_inbounds function Base.getindex(r::IdOffsetRange{HalfInt}, s::AbstractUnitRange{HalfInt})
    return r.parent[unsafeUnitRangeInt(s .- r.offset)] .+ r.offset
end
@propagate_inbounds function Base.getindex(r::IdOffsetRange{HalfInt}, s::IdOffsetRange{HalfInt})
    p = r.parent[unsafeUnitRangeInt(s .- r.offset)]
    return IdOffsetRange(p, r.offset)
end
for DT in [:Integer,:HalfInteger]
    @eval @propagate_inbounds function Base.getindex(r::IdOffsetRange, i::$DT)
    	r.parent[Integer(i - r.offset)] + r.offset
    end
    @eval @propagate_inbounds function Base.getindex(r::IdOffsetRange, s::AbstractUnitRange{<:$DT})
        return r.parent[map(Integer, s .- r.offset)] .+ r.offset
    end
end

@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::IdOffsetRange)
    p = r.parent[map(Integer,s .- r.offset)]
    return IdOffsetRange(p, r.offset)
end

# Optimizations
function Base.checkindex(::Type{Bool}, r::IdOffsetRange, i::Real)
    Base.checkindex(Bool, r.parent, i - r.offset) && isinteger(i - r.offset)
end
function Base.checkindex(::Type{Bool}, r::Base.Slice{<:IdOffsetRange}, i::Real)
    Base.checkindex(Bool, r.indices, i - r.offset)
end

function Base.isassigned(r::IdOffsetRange, I::Integer...)
    (checkbounds(Bool, r, I...) && isinteger(r.offset)) || return false
    J = to_parentindices(axes(r),I)
    isassigned(parent(r), J...)
end
function Base.isassigned(r::IdOffsetRange, I::Real...)
    (checkbounds(Bool, r, I...) && indicescompatible((r.offset,), I)) || return false
    J = to_parentindices(axes(r),I)
    isassigned(parent(r), J...)
end

function Base.hash(r::IdOffsetRange, h::UInt)
    h = hash(parent(r), h)
    h = hash(r.offset, h)
    h = hash(:IdOffsetRange, h)
    return h
end


text:
using BlackBoxOptim

# Another problem from Hans W. Borchers, constrained 15-points:
function points15(x)
    d = 2.0
    for i = 1:14, j = (i+1):15
        s = (x[i]-x[j])^2 + (x[i+15]-x[j+15])^2 + (x[i+30]-x[j+30])^2
        s = sqrt(s)
        if s < d
          d = s
        end
    end
    return -d
end

Dims = 15 * 3

# The best known minimum is around -0.62 but we do not clearly know if there are better
# ones so we aim somewhat lower...
prob = BlackBoxOptim.minimization_problem(points15, "Constrained points 15",
        (0.0, 1.0), 45, -0.64)

# This is a hard problem so run for some time...
MaxMinutes = 0.50

# To run XNES et al on it we need to add a penalty for going outside the (0,1) box.
function penalty(x, range)
  lower = x .< range[1]
  higher = x .> range[2]
  100 * (norm(range[1] - x[lower]) + norm(x[higher] - range[2]))
end

penalized_points15(x) = points15(x) + penalty(x, (0.0, 1.0))

# Run XNES
result = bboptimize(penalized_points15; SearchRange = (0.0, 1.0), NumDimensions = Dims,
  MaxTime = MaxMinutes * 60, Method = :xnes)
best_xnes = best_candidate(result)

# Run DXNES
result = bboptimize(penalized_points15; SearchRange = (0.0, 1.0), NumDimensions = Dims,
  MaxTime = MaxMinutes * 60, Method = :dxnes)
best_dxnes = best_candidate(result)

# Run the default (Adaptive DE) method for same amount of time on penalized function
result = bboptimize(penalized_points15; SearchRange = (0.0, 1.0), NumDimensions = Dims,
  MaxTime = MaxMinutes * 60)
best_de_pen = best_candidate(result)

# But DE need no penalization since it respects search range bounds:
result = bboptimize(points15; SearchRange = (0.0, 1.0), NumDimensions = Dims,
  MaxTime = MaxMinutes * 60)
best_de = best_candidate(result)

# Run several methods several times:
BlackBoxOptim.repeated_bboptimize(2, prob, Dims, [
  :generating_set_search,
  :adaptive_de_rand_1_bin_radiuslimited,
  :random_search,
  ],
  MaxMinutes * 60)

println("xNES points15 fitness = ", points15(best_xnes))
println("dxNES points15 fitness = ", points15(best_dxnes))
println("Adaptive DE (default, penalized) points15 fitness = ", points15(best_de_pen))
println("Adaptive DE (default) points15 fitness = ", points15(best_de))

ext: jl


text:
using Stheno: BlockMean, BlockCrossKernel, BlockKernel, ConstantMean, ZeroMean
using Stheno: FiniteMean, FiniteKernel, FiniteCrossKernel, map, pairwise
using FillArrays

@testset "block" begin

    # Test BlockMean.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x1, x2 = randn(rng, N), randn(rng, N′)
        X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        μ1, μ2 = ConstantMean(1.0), ZeroMean{Float64}()
        μ, Dx, DX = BlockMean([μ1, μ2]), BlockData([x1, x2]), BlockData([X1, X2])

        @test μ == BlockMean(μ1, μ2)
        @test length(μ) == Inf
        @test map(μ, DX) == vcat(map(μ1, X1), map(μ2, X2))
        @test map(μ, Dx) == vcat(map(μ1, x1), map(μ2, x2))

        # mean_function_tests(μ, BlockData([X1, X2]))
        # mean_function_tests(μ, BlockData([x1, x2]))
        # unary_map_tests(μ, BlockData([X1, X2]))

        # Tests for finite case.
        μ1f, μ2f = FiniteMean(μ1, X1), FiniteMean(μ2, X2)
        cat_μ = BlockMean([μ1f, μ2f])
        @test eachindex(cat_μ) == BlockData([eachindex(μ1f), eachindex(μ2f)])
        @test AbstractVector(cat_μ) ==
            BlockVector([AbstractVector(μ1f), AbstractVector(μ2f)])
    end

    # Test BlockCrossKernel.
    let
        rng, N1, N2, N1′, N2′, D = MersenneTwister(123456), 5, 6, 2, 7, 8
        X0, X0′ = ColsAreObs(randn(rng, D, N1)), ColsAreObs(randn(rng, D, N2))
        X1, X1′ = ColsAreObs(randn(rng, D, N1′)), ColsAreObs(randn(rng, D, N2′))
        X2, X2′ = ColsAreObs(randn(rng, D, N1)), ColsAreObs(randn(rng, D, N2))
        k11, k12, k21, k22 =  EQ(), ZeroKernel{Float64}(), ZeroKernel{Float64}(), EQ()
        k = BlockCrossKernel([k11 k12; k21 k22])

        @test k == k

        @test size(k) == (Inf, Inf)
        @test size(k, 1) == Inf
        @test size(k, 2) == Inf

        row1 = hcat(pairwise(k11, X0, X0′), pairwise(k12, X0, X1′))
        row2 = hcat(pairwise(k21, X1, X0′), pairwise(k22, X1, X1′))
        @test pairwise(k, BlockData([X0, X1]), BlockData([X0′, X1′])) == vcat(row1, row2)

        # cross_kernel_tests(k, [X0, X0′], [X2, X2′], [X1, X1′])

        # Tests for finite case.
        k11f, k12f = FiniteCrossKernel(k11, X1, X1), FiniteCrossKernel(k12, X1, X0)
        k21f, k22f = FiniteCrossKernel(k21, X0, X1), FiniteCrossKernel(k22, X0, X0)
        kf = BlockCrossKernel([k11f k12f; k21f k22f])
        @test eachindex(kf, 1) == BlockData([eachindex(k11f, 1), eachindex(k22f, 1)])
        @test eachindex(kf, 2) == BlockData([eachindex(k11f, 2), eachindex(k22f, 2)])
        @test AbstractMatrix(kf) isa BlockMatrix
        @test AbstractMatrix(kf) ==
            BlockMatrix(reshape(AbstractMatrix.([k11f, k21f, k12f, k22f]), 2, 2))
    end

    # Test BlockKernel.
    let
        rng, N1, N2, N1′, N2′, D = MersenneTwister(123456), 5, 6, 3, 4, 2
        X0, X0′ = ColsAreObs(randn(rng, D, N1)), ColsAreObs(randn(rng, D, N2))
        X1, X1′ = ColsAreObs(randn(rng, D, N1′)), ColsAreObs(randn(rng, D, N2′))
        X2, X2′ = ColsAreObs(randn(rng, D, N1)), ColsAreObs(randn(rng, D, N2))

        # Construct BlockKernel.
        k11, k22, k12 = EQ(), EQ(), ZeroKernel{Float64}()
        ks_off = Matrix{Stheno.CrossKernel}(undef, 2, 2)
        ks_off[1, 2] = k12
        k = BlockKernel([k11, k22], ks_off)

        @test size(k) == (Inf, Inf)
        @test size(k, 1) == Inf
        @test size(k, 2) == Inf

        row1 = hcat(pairwise(k11, X0), pairwise(k12, X0, X0′))
        row2 = hcat(Zeros{Float64}(N2, N1), pairwise(k22, X0′))
        @test pairwise(k, BlockData([X0, X0′])) == vcat(row1, row2)

        # Compute xcov for BlockKernel with infinite kernels.
        manual = vcat(hcat(pairwise(k11, X1, X1′), pairwise(k12, X1, X2′)),
                      hcat(pairwise(k12, X2, X1′), pairwise(k22, X2, X2′)),)
        @test pairwise(k, BlockData([X1, X2]), BlockData([X1′, X2′])) == manual

        # kernel_tests(k, [X0, X0′], [X2, X2′], [X1, X1′])

        # Tests for finite case.
        k11f, k22f = FiniteKernel(k11, X0), FiniteKernel(k22, X1)
        k12f = FiniteCrossKernel(k12, X0, X1)
        ks_off_f = Matrix{Stheno.CrossKernel}(undef, 2, 2)
        ks_off_f[1, 2] = k12f
        k = BlockKernel([k11f, k22f], ks_off_f)

        @test eachindex(k) == BlockData([eachindex(k11f), eachindex(k22f)])
        Σ11, Σ21 = AbstractMatrix(k11f), AbstractMatrix(k12f)'
        Σ12, Σ22 = AbstractMatrix(k12f), AbstractMatrix(k22f)
        @test AbstractMatrix(k) ==
            BlockMatrix(reshape([Σ11, Σ21, Σ12, Σ22], 2, 2))
    end
end


text:
/* void expression as statement */

int main() {
  foo();
  return 0 ;

}

void foo() {
   printString("foo");
   return;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
using GAPLib
using Test

@testset "LoadSymbol" begin
    data = loadGAP(:a05200)
    @test data.name == "a05200"
    @test data.agents[3].capacity == 671
    @test length(data.agents) == 5
    @test length(data.jobs) == 200
    @test data.lb == typemin(Int64)
    @test data.ub == typemax(Int64)
    @test_nowarn println(data)
    @test_nowarn println(data.agents[1])
    @test_nowarn println(data.jobs[1])
end

@testset "ErrorSymbol" begin
    data = loadGAP(:notaninstance)
    @test data === nothing
end

@testset "LoadString" begin
    data = loadGAP(joinpath(pkgdir(GAPLib), "test/data/a05100"))
    @test data.name == "a05100"
    @test data.agents[3].capacity == 342
    @test length(data.agents) == 5
    @test length(data.jobs) == 100
    @test data.lb == typemin(Int64)
    @test data.ub == typemax(Int64)
    @test_nowarn println(data)
end

@testset "ErrorString" begin
    data = loadGAP("notaninstance")
    @test data === nothing
end


text:
using UnicodePlots, ReferenceTests, Test
using ReferenceTests: BeforeAfterFull
using StatsBase
using SparseArrays
using Random: seed!
using Dates: Date, Day

tests = [
    "tst_common.jl",
    "tst_graphics.jl",
    "tst_canvas.jl",
    "tst_plot.jl",
    "tst_barplot.jl",
    "tst_histogram.jl",
    "tst_scatterplot.jl",
    "tst_lineplot.jl",
    "tst_spy.jl",
    "tst_boxplot.jl",
    "tst_heatmap.jl"
]

for test in tests
    @testset "$test" begin
        include(test)
    end
end


text:
using QuantEcon
import QuantEcon: simulate

mutable struct Model{TF <: AbstractFloat,
                    TM <: AbstractMatrix{TF},
                    TV <: AbstractVector{TF}}
    β::TF
    Π::TM
    G::TV
    Θ::TV
    transfers::Bool
    U::Function
    Uc::Function
    Ucc::Function
    Un::Function
    Unn::Function
    n_less_than_one::Bool
end


text:
@testset "Connectivity" begin
    g6 = smallgraph(:house)
    gx = PathGraph(4)
    add_vertices!(gx, 10)
    add_edge!(gx, 5, 6)
    add_edge!(gx, 6, 7)
    add_edge!(gx, 8, 9)
    add_edge!(gx, 10, 9)

    for g in testgraphs(gx)
        @test @inferred(!is_connected(g))
        cc = @inferred(connected_components(g))
        label = zeros(eltype(g), nv(g))
        @inferred(LightGraphs.connected_components!(label, g))
        @test label[1:10] == [1, 1, 1, 1, 5, 5, 5, 8, 8, 8]
        import LightGraphs: components, components_dict
        cclab = @inferred(components_dict(label))
        @test cclab[1] == [1, 2, 3, 4]
        @test cclab[5] == [5, 6, 7]
        @test cclab[8] == [8, 9, 10]
        @test length(cc) >= 3 && sort(cc[3]) == [8, 9, 10]
    end
    for g in testgraphs(g6)
        @test @inferred(is_connected(g))
    end

    g10 = DiGraph(4)
    add_edge!(g10, 1, 3)
    add_edge!(g10, 2, 4)
    for g in testdigraphs(g10)
        @test @inferred(is_bipartite(g))
    end
    add_edge!(g10, 1, 4)
    for g in testdigraphs(g10)
        @test @inferred(is_bipartite(g))
    end

    g10 = DiGraph(20)
    for g in testdigraphs(g10)
        for m = 1:50
            i = rand(1:10)
            j = rand(11:20)
            if rand() < 0.5
                i, j = j, i
            end
          if !has_edge(g, i, j)
              add_edge!(g, i, j)
              @test @inferred(is_bipartite(g))
          end
        end
    end

    # graph from https://en.wikipedia.org/wiki/Strongly_connected_component
    h = DiGraph(8)
    add_edge!(h, 1, 2); add_edge!(h, 2, 3); add_edge!(h, 2, 5);
    add_edge!(h, 2, 6); add_edge!(h, 3, 4); add_edge!(h, 3, 7);
    add_edge!(h, 4, 3); add_edge!(h, 4, 8); add_edge!(h, 5, 1);
    add_edge!(h, 5, 6); add_edge!(h, 6, 7); add_edge!(h, 7, 6);
    add_edge!(h, 8, 4); add_edge!(h, 8, 7)
    for g in testdigraphs(h)
        @test @inferred(is_weakly_connected(g))
        scc = @inferred(strongly_connected_components(g))
        wcc = @inferred(weakly_connected_components(g))

        @test length(scc) == 3 && sort(scc[3]) == [1, 2, 5]
        @test length(wcc) == 1 && length(wcc[1]) == nv(g)
    end

    function scc_ok(graph)
      """Check that all SCC really are strongly connected"""
      scc = @inferred(strongly_connected_components(graph))
      scc_as_subgraphs = map(i -> graph[i], scc)
      return all(is_strongly_connected, scc_as_subgraphs)
    end

    # the two graphs below are isomorphic (exchange 2 <--> 4)
    h = DiGraph(4);  add_edge!(h, 1, 4); add_edge!(h, 4, 2); add_edge!(h, 2, 3); add_edge!(h, 1, 3);
    for g in testdigraphs(h)
      @test scc_ok(g)
    end

    h2 = DiGraph(4); add_edge!(h2, 1, 2); add_edge!(h2, 2, 4); add_edge!(h2, 4, 3); add_edge!(h2, 1, 3);
    for g in testdigraphs(h2)
      @test scc_ok(g)
    end

    h = DiGraph(6)
    add_edge!(h, 1, 3); add_edge!(h, 3, 4); add_edge!(h, 4, 2); add_edge!(h, 2, 1)
    add_edge!(h, 3, 5); add_edge!(h, 5, 6); add_edge!(h, 6, 4)
    for g in testdigraphs(h)
      scc = @inferred(strongly_connected_components(g))
      @test length(scc) == 1 && sort(scc[1]) == [1:6;]
    end
    # tests from Graphs.jl
    h = DiGraph(4)
    add_edge!(h, 1, 2); add_edge!(h, 2, 3); add_edge!(h, 3, 1); add_edge!(h, 4, 1)
    for g in testdigraphs(h)
      scc = @inferred(strongly_connected_components(g))
      @test length(scc) == 2 && sort(scc[1]) == [1:3;] && sort(scc[2]) == [4]
    end

    h = DiGraph(12)
    add_edge!(h, 1, 2); add_edge!(h, 2, 3); add_edge!(h, 2, 4); add_edge!(h, 2, 5);
    add_edge!(h, 3, 6); add_edge!(h, 4, 5); add_edge!(h, 4, 7); add_edge!(h, 5, 2);
    add_edge!(h, 5, 6); add_edge!(h, 5, 7); add_edge!(h, 6, 3); add_edge!(h, 6, 8);
    add_edge!(h, 7, 8); add_edge!(h, 7, 10); add_edge!(h, 8, 7); add_edge!(h, 9, 7);
    add_edge!(h, 10, 9); add_edge!(h, 10, 11); add_edge!(h, 11, 12); add_edge!(h, 12, 10)

    for g in testdigraphs(h)
      scc = @inferred(strongly_connected_components(g))
      @test length(scc) == 4
      @test sort(scc[1]) == [7, 8, 9, 10, 11, 12]
      @test sort(scc[2]) == [3, 6]
      @test sort(scc[3]) == [2, 4, 5]
      @test scc[4] == [1]
    end
    # Test examples with self-loops from
    # Graph-Theoretic Analysis of Finite Markov Chains by J.P. Jarvis & D. R. Shier

    # figure 1 example
    fig1 = spzeros(5, 5)
    fig1[[3, 4, 9, 10, 11, 13, 18, 19, 22, 24]] = [.5, .4, .1, 1., 1., .2, .3, .2, 1., .3]
    fig1 = DiGraph(fig1)
    scc_fig1 = Vector[[2, 5], [1, 3, 4]]

    # figure 2 example
    fig2 = spzeros(5, 5)
    fig2[[3, 10, 11, 13, 14, 17, 18, 19, 22]] = 1
    fig2 = DiGraph(fig2)

    # figure 3 example
    fig3 = spzeros(8, 8)
    fig3[[1, 7, 9, 13, 14, 15, 18, 20, 23, 27, 28, 31, 33, 34, 37, 45, 46, 49, 57, 63, 64]] = 1
    fig3 = DiGraph(fig3)
    scc_fig3 = Vector[[3, 4], [2, 5, 6], [8], [1, 7]]
    fig3_cond = DiGraph(4);
    add_edge!(fig3_cond, 4, 3); add_edge!(fig3_cond, 2, 1)
    add_edge!(fig3_cond, 4, 1); add_edge!(fig3_cond, 4, 2)

    # construct a n-number edge ring graph (period = n)
    n = 10
    n_ring_m = spdiagm(ones(n - 1), 1, n, n); n_ring_m[end, 1] = 1
    n_ring = DiGraph(n_ring_m)

    n_ring_shortcut = copy(n_ring); add_edge!(n_ring_shortcut, 1, 4)

    # figure 8 example
    fig8 = spzeros(6, 6)
    fig8[[2, 10, 13, 21, 24, 27, 35]] = 1
    fig8 = DiGraph(fig8)

    @test Set(@inferred(strongly_connected_components(fig1))) == Set(scc_fig1)
    @test Set(@inferred(strongly_connected_components(fig3))) == Set(scc_fig3)

    @test @inferred(period(n_ring)) == n
    @test @inferred(period(n_ring_shortcut)) == 2

    @test @inferred(condensation(fig3)) == fig3_cond

    @test @inferred(attracting_components(fig1)) == Vector[[2, 5]]
    @test @inferred(attracting_components(fig3)) == Vector[[3, 4], [8]]

    g10 = StarGraph(10)
    for g in testgraphs(g10)
      @test @inferred(neighborhood(g, 1, 0)) == [1]
      @test length(@inferred(neighborhood(g, 1, 1))) == 10
      @test length(@inferred(neighborhood(g, 2, 1))) == 2
      @test length(@inferred(neighborhood(g, 1, 2))) == 10
      @test length(@inferred(neighborhood(g, 2, 2))) == 10
    end
    g10 = StarDiGraph(10)
    for g in testdigraphs(g10)
        @test @inferred(neighborhood(g10, 1, 0, dir=:out)) == [1]
        @test length(@inferred(neighborhood(g, 1, 1, dir=:out))) == 10
        @test length(@inferred(neighborhood(g, 2, 1, dir=:out))) == 1
        @test length(@inferred(neighborhood(g, 1, 2, dir=:out))) == 10
        @test length(@inferred(neighborhood(g, 2, 2, dir=:out))) == 1
        @test @inferred(neighborhood(g, 1, 0, dir=:in)) == [1]
        @test length(@inferred(neighborhood(g, 1, 1, dir=:in))) == 1
        @test length(@inferred(neighborhood(g, 2, 1, dir=:in))) == 2
        @test length(@inferred(neighborhood(g, 1, 2, dir=:in))) == 1
        @test length(@inferred(neighborhood(g, 2, 2, dir=:in))) == 2
    end
    @test @inferred(!isgraphical([1, 1, 1]))
    @test @inferred(isgraphical([2, 2, 2]))
    @test @inferred(isgraphical(fill(3, 10)))
end


text:
using Test
using HOI4CombatSimulator
using Statistics

macro rep(n::Int, exp::Expr)
    quote
        resp_array = Array{Any, 1}(undef, $n)
        for i in 1:$n
            resp_array[i] = $exp
        end
        m = length(resp_array)
        n = length(resp_array[1])
        mat = zeros(m, n)
        for i in 1:m
            for j in 1:n
                mat[i, j] = resp_array[i][j]
            end
        end
        mat
    end
end

@testset "HOI4CombatSimulator" begin
    @testset "get_unit_list" begin
        unit_list = HOI4CombatSimulator.get_unit_list()
        @test length(unit_list) > 0
    end

    @testset "wiki_example" begin
        HP_loss_Org_loss_mat = @rep 1000 HOI4CombatSimulator.fire(
            evade=500, HPPercent=0.99, ShooterHardAtk=0,
            ShooterSoftAtk=1000, VictimHardness=0,
            ShooterPierce=20, ShooterArmor=10,
            VictimPierce=1, VictimArmor=0
        ) # (n, 2)
        HP_loss_arr = HP_loss_Org_loss_mat[:, 1]
        Org_loss_arr = HP_loss_Org_loss_mat[:, 2]
        diff = abs(mean(Org_loss_arr) - 3.9)
        tol = (std(Org_loss_arr) / sqrt(1000))*5 # 5 sigma
        @test  diff < tol
    end

    @testset "conv" begin
        x = HOI4CombatSimulator.Inference.fast_conv(ones(2)/2, 2)
        @test maximum(abs.(x .- [0.25, 0.5, 0.25])) < 1e-6
    end

    @testset "DiceSumDistribution" begin
        d6x2 = HOI4CombatSimulator.Inference.dice_sum_distribution(6, 2)
        @test d6x2.support[1] == 2
        @test d6x2.support[end] == 12
        @test length(d6x2.support) == 11
        @test abs(d6x2.p[1] - 1/36) < 1e-6

        d6x4 = HOI4CombatSimulator.Inference.dice_sum_distribution(6, 4)
        @test d6x4.p[end] == (1/6)^4
        @test abs(sum(d6x4.p)- 1.) < 1e-6
    end

    @testset "prob_fire_loss" begin
        HP_Org_loss_mat = HOI4CombatSimulator.Inference.prob_fire(104.0, 160.0, 6)
        @test size(HP_Org_loss_mat, 1) > 1
        @test size(HP_Org_loss_mat, 2) > 1
        # TODO: detailed test
    end
end


text:
using Pumas, StatsFuns
data = read_pumas(example_nmtran_data("pain_remed"),
    cvs = [:arm, :dose, :conc, :painord,:remed],event_data=false)

import StatsFuns.logistic
binary_model = @model begin
    @param begin
        intercept ∈ RealDomain(init=0.001)
        tvslope ∈ RealDomain(init=0.0001)
        Ω ∈ VectorDomain(1)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates arm dose

    @pre begin
        rx = dose > 0 ? 1 : 0
        slope = tvslope*rx
        logit = intercept + slope + η[1]
    end

    @derived begin
        pain = logistic(logit)
        dv ~ Bernoulli(logistic(logit))
    end
end

@pre begin
    logit = intercept + tvslope*(dose > 0 ? 1 : 0) + η[1]
    logit = intercept + tvslope*Int(dose > 0) + η[1]
end

param = (
    intercept = 0.001,
    tvslope = 0.0001,
    Ω = [1.0]
    )
res = fit(binary_model,data,param,Pumas.LaplaceI())

sim = simobs(binary_model,data,res.param)
simdf = DataFrame(sim, include_events=false)
first(simdf,6) # Print only the first 6 rows

pop = Population(map(i -> Subject(id=i,cvs=(dose = i.*10,),time=[0.0]),1:10))

poisson_model = @model begin
  @param begin
    tvbase ∈ RealDomain(init=3.0, lower=0.1)
    d50 ∈ RealDomain(init=50, lower=0.1)
    Ω  ∈ PSDDomain(fill(0.1, 1, 1))
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    baseline = tvbase*exp(η[1])
  end

  @covariates dose

  @derived begin
    dv ~ @. Poisson(baseline*(1-dose/(dose + d50)))
  end
end

sim = simobs(poisson_model,pop)
simdf = DataFrame(sim, include_events=false)

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
export @path_str

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

abstract type AbstractPath ; end

"""
    path"relative/path/to/resource"

A `RelPath` is a *key* into a hierarchical string-indexed tree datastructure,
with each component indexing one level of the hierarchy.

As a key, the resource referred to by a path may or may not exist.
Conversely, `BlobTree` and `Blob` refer to the actual data stored with a given
key.
"""
struct RelPath <: AbstractPath
    components::Vector{String}
end

RelPath(::AbstractString) = error("RelPath(::String) is not defined to avoid ambiguities between operating systems. Use  the `path\"...\"` string macro for path literals.")
RelPath(components::AbstractVector=String[]) = RelPath(convert(Vector{String}, components))

# Path manipulation.
function Base.joinpath(path::RelPath, xs::AbstractString...)
    for x in xs
        if '/' in x || '\\' in x
            error("Path components cannot contain '/' or '\\' (got $(repr(x)))")
        end
    end
    RelPath(vcat(path.components, xs...))
end
Base.joinpath(path::RelPath, xs::RelPath...) =
    RelPath(vcat(path.components, [x.components for x in xs]...))

Base.isempty(p::RelPath) = isempty(p.components)

Base.basename(path::RelPath) = isempty(path) ? "" : last(path.components)
Base.dirname(path::RelPath) = RelPath(path.components[1:end-1])

Base.print(io::IO, p::RelPath) = print(io, join(p.components, '/'))
Base.show(io::IO,  p::RelPath) = print(io, "path", repr(string(p)))

function Base.startswith(a::RelPath, b::RelPath)
    return length(a.components) >= length(b.components) &&
           a.components[1:length(b.components)] == b.components
end

function Base.:(==)(p1::RelPath, p2::RelPath)
    return p1.components == p2.components
end

function Base.isless(p1::RelPath, p2::RelPath)
    for (a,b) in zip(p1.components, p2.components)
        isless(a, b) && return true
        isless(b, a) && return false
    end
    return length(p1.components) < length(p2.components)
end

"""
    A tool to write relative path literals succinctly
"""
macro path_str(str)
    # For OS-independence, just allow / and \ as equivalent; normalize to '/' when printing.
    # As a side effect, this disallows use of \ in individual *path literal
    # components* on unix. But that would be perverse, right?
    components = isempty(str) ? [] : split(str, ('/', '\\'))
    RelPath(components)
end

#-------------------------------------------------------------------------------
"""
    An AbsPath is the *key* into a hierarchical tree index, relative to some root.

As a *key*, the resource pointed to by this key may or may not exist.
"""
struct AbsPath{Root} <: AbstractPath
    root::Root
    path::RelPath
end

AbsPath(root, path::AbstractString) = AbsPath(root, RelPath(path))

function Base.show(io::IO, ::MIME"text/plain", path::AbsPath)
    print(io, "AbsPath ", path.path, " @ ", path.root)
end

Base.mkdir(p::AbsPath; kws...) = mkdir(p.root, p.path; kws...)
Base.rm(p::AbsPath; kws...) = rm(p.root, p.path; kws...)

Base.read(p::AbsPath, ::Type{T}) where {T} = read(p.root, p.path, T)
Base.abspath(p::AbsPath) = p


text:
struct PlatonicSolid{N, T} <: AbstractChargeCloud
    locations::SVector{N, CartesianPoint{T}}
end

function get_shell_structure_locations(structure::NTuple{A,Tuple{Int,T,T}}, center::CartesianPoint{T}, length::T = T(1))::SVector where {A,T <: SSDFloat}
    N::Int = sum(map(p -> p[1], structure))
    locations::Vector{CartesianPoint{T}} = Vector{CartesianPoint{T}}(undef, N)
    i::Int = 1
    for s in structure
        n::Int = s[1]
        z::T = length * s[2]
        ρ::T = length * sqrt(1 - s[2]^2)
        for φ in range(s[3], step = T(2π/n), length = n)
            sinφ::T, cosφ::T = sincos(φ)
            locations[i] = center + CartesianVector{T}(ρ*cosφ,ρ*sinφ,z)
            i+=1
        end
    end
    SVector{N, CartesianPoint{T}}(locations)
end

get_vertices(::Type{PlatonicSolid{N, T}}) where {N, T} = N

"""
    struct PointCharge{T} <: AbstractChargeCloud

Struct which defines a single point-like charge carrier.

## Fields
* `locations::SVector{1, CartesianPoint{T}}`: Position of the charge carrier, saved as single entry of a `Vector`.

## Example
```julia
center = CartesianPoint{T}(0,0,0)
pc = PointCharge(center) # Constructor: creates a PointCharge around (0,0,0)
pc.locations             # Array that only contains the center point
```
"""
const PointCharge{T} = PlatonicSolid{1, T}
PointCharge(center::Vector{CartesianPoint{T}}, args...) where {T} = PointCharge{T}(center)
PointCharge(center::CartesianPoint{T}, args...) where {T} = PointCharge{T}([center])

const Tetrahedron{T} = PlatonicSolid{4, T}
function Tetrahedron(center::CartesianPoint{T}, length::Real = T(1))::Tetrahedron{T} where {T}
    structure::NTuple{2, Tuple{Int,T,T}} = ((1,T(1),T(0)), (3,T(-1/3),T(0)))
    PlatonicSolid(get_shell_structure_locations(structure, center, T(length)))
end

const Hexahedron{T} = PlatonicSolid{8, T}
function Hexahedron(center::CartesianPoint{T}, length::Real = T(1))::Hexahedron{T} where {T}
    structure::NTuple{4, Tuple{Int,T,T}} = ((1,T(1),T(0)),(3,T(1/3),T(0)),(3,T(-1/3),T(π/3)),(1,T(-1),T(0)))
    PlatonicSolid(get_shell_structure_locations(structure, center, T(length)))
end

const Octahedron{T} = PlatonicSolid{6, T}
function Octahedron(center::CartesianPoint{T}, length::Real = T(1))::Octahedron{T} where {T}
    structure::NTuple{3, Tuple{Int,T,T}} = ((1,T(1),T(0)), (4,T(0),T(0)), (1,T(-1),T(0)))
    PlatonicSolid(get_shell_structure_locations(structure, center, T(length)))
end

const Icosahedron{T} = PlatonicSolid{12, T}
function Icosahedron(center::CartesianPoint{T}, length::Real = T(1))::Icosahedron{T} where {T}
    structure::NTuple{4, Tuple{Int,T,T}} = ((1,T(1),T(0)), (5,T(sqrt(1/5)),T(0)), (5,T(-sqrt(1/5)),T(π/5)), (1,T(-1),T(0)))
    PlatonicSolid(get_shell_structure_locations(structure, center, T(length)))
end

const Dodecahedron{T} = PlatonicSolid{20, T}
function Dodecahedron(center::CartesianPoint{T}, length::Real = T(1))::Dodecahedron{T} where {T}
    zmax::T = sqrt((5 + 2*sqrt(5))/15)
    tmp::T = 1 - 2*(1 - zmax^2) * sin(π/5)^2
    zmin::T = zmax * tmp - sqrt(1 - zmax^2) * sqrt(1 - tmp^2)
    structure::NTuple{4, Tuple{Int,T,T}} = ((5,zmax,T(0)),(5,zmin,T(0)),(5,-zmin,T(π/5)),(5,-zmax,T(π/5)))
    PlatonicSolid(get_shell_structure_locations(structure, center, T(length)))
end


text:
# Diagonal operations

using NumericExtensions
using Base.Test

function zero_diag(x::Matrix)
    xc = copy(x)
    for i in 1 : min(size(x,1), size(x,2))
        xc[i,i] = 0
    end
    xc
end

@test diagm(3, 2.) == [2. 0. 0.; 0. 2. 0.; 0. 0. 2.]
@test eltype(diagm(3, 2.)) == Float64

x = rand(3, 3)
xc = copy(x)

@test set_diag(x, 2) == zero_diag(x) + diagm(3, 2.)
@test set_diag(x, [1. 2. 3.]) == zero_diag(x) + diagm([1., 2., 3.])
@test set_diag(x, [1. 2. 3.], 2.) == zero_diag(x) + diagm([2., 4., 6.])

@test xc == x

@test add_diag(x, 2) == x + diagm(3, 2.)
@test add_diag(x, [1. 2. 3.]) == x + diagm([1., 2., 3.])
@test add_diag(x, [1. 2. 3.], 2.) == x + diagm([2., 4., 6.])

@test xc == x


text:
module Utils

using Distributions
using SpecialFunctions
import SpecialFunctions: logbeta, logfactorial
using LinearAlgebra
using Clustering
using RCall
using ProgressMeter
using Combinatorics

function one_hot_encode(assignments::Vector{Int64})
    n = length(assignments)
    youngest_cluster = maximum(assignments)
    encoding_matrix = zeros(Int64, youngest_cluster, n)
    for observation = 1:n
        cluster = assignments[observation]
        encoding_matrix[cluster, observation] = 1
    end
    return encoding_matrix
end

function compute_co_occurrence_matrix(markov_chain::Matrix{Int64})
    num_instances = size(markov_chain)[2]
    weights = ones(Float64, num_instances)
    return compute_co_occurrence_matrix(markov_chain, weights)
end

function compute_co_occurrence_matrix(markov_chain::Matrix{Int64}, weights::Vector{Float64})
    n = size(markov_chain)[1]
    co_occurrence_matrix = zeros(Int64, n, n)
    num_instances = size(markov_chain)[2]
    @showprogress for i = 1:num_instances
        assignment = markov_chain[:, i]
        ohe_assignment = one_hot_encode(assignment)
        instance_co_occurrence_matrix = transpose(ohe_assignment) * ohe_assignment
        weight = weights[i]
        co_occurrence_matrix += weight*instance_co_occurrence_matrix
    end
    max_value = maximum(co_occurrence_matrix)
    co_occurrence_matrix /= max_value
    return co_occurrence_matrix
end

eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))

logbeta(parameters::Vector{T}) where {T<:Real} = logbeta(parameters[1], parameters[2])

function logmvbeta(parameters::Vector{T}) where {T <: Real}
    return sum(loggamma.(parameters)) - loggamma(sum(parameters))
end

function log_multinomial_coeff(counts::Vector{Int64})
    return logfactorial(sum(counts)) - sum(logfactorial.(counts))
end

logfactorial(x::Float64) = logfactorial(Int(x))

function gumbel_max(objects::Vector{Int64}, log_weights::Vector{Float64})
    gumbel_values = rand(Gumbel(0, 1), length(log_weights))
    index = argmax(gumbel_values + log_weights)
    return (objects[index], log_weights[index])
end

function gumbel_max(num_draws::Int64, log_weights::Vector{Float64})
    indices = Vector{Int64}(1:length(log_weights))
    resampled_indices = Vector{Int64}(undef, num_draws)
    for i = 1:num_draws
        (sampled_index, log_weight) = gumbel_max(indices, log_weights)
        resampled_indices[i] = sampled_index
    end
    return resampled_indices
end

isnothing(::Any) = false
isnothing(::Nothing) = true

function compute_normalized_weights(log_weights)
    num_weights = length(log_weights)
    max_value = maximum(log_weights)
    shifted_log_weights = log_weights .- max_value
    weights = exp.(shifted_log_weights)
    normalized_weights = weights ./ sum(weights)
    @assert !any(isnan.(normalized_weights))
    return normalized_weights
end

function effective_sample_size(log_weights)
    normalized_weights = compute_normalized_weights(log_weights)
    ess = 1/sum(normalized_weights.^2)
    @assert !isnan(ess)
    return ess
end

function map_estimate(instances, log_likelihoods; num_burn_in=0)
    @assert size(instances)[2] === length(log_likelihoods)
    burned_in_log_likelihoods = log_likelihoods[num_burn_in+1:end]
    max_value = maximum(burned_in_log_likelihoods)
    burned_in_instances = instances[:, num_burn_in+1:end]
    map_index = findfirst(burned_in_log_likelihoods .=== max_value)
    return burned_in_instances[:, map_index]
end

function num_clusters(instances)
    return vec(mapslices(a -> length(unique(a)), instances; dims=[1]))
end

function mean_num_clusters(instances, weights)
    number_of_clusters = num_clusters(instances)
    return sum(weights .* number_of_clusters)/sum(weights)
end

function mean_num_clusters(instances::Matrix)
    return mean(num_clusters(instances))
end

function mean_num_clusters(instances::Vector)
    return mean_num_clusters(instances, [1])
end

function hist(v::AbstractVector, edg::AbstractVector)
    n = length(edg)-1
    h = zeros(Int, n)
    for x in v
        i = searchsortedfirst(edg, x)-1
        if 1 <= i <= n
            h[i] += 1
        end
    end
    return h
end

function mvt_logpdf(location, scale, dof, datum)
    log_det = logdet(scale)
    dim = length(datum)
    log_likelihood = -(dof + dim)log(1 + transpose(datum - location)*inv(scale)*(datum - location)/dof)/2
    log_likelihood += loggamma((dof + dim)/2) - loggamma(dof/2) - dim*log(dof*pi)/2 - log_det/2
    return log_likelihood
end

function adjusted_rand_index(c1::Vector{Int64}, c2::Vector{Int64})
    return Clustering.randindex(c1,c2)[1]
end

function ari_over_markov_chain(true_clustering::Vector{Int64}, markov_chain::Matrix{Int64})
    posterior_ari = vec(mapslices(z->adjusted_rand_index(true_clustering, z), markov_chain, dims=1))
    return posterior_ari
end

function minbinder(posterior_similarity_matrix::Matrix{Float64}, iterations::Matrix{Int64})
    draws = Matrix(transpose(iterations))
    R"library('mcclust')"
    R"library('mcclust.ext')"
    @rput posterior_similarity_matrix
    @rput draws
    R"binder_clustering <- minbinder.ext(posterior_similarity_matrix, draws, method='draws')"
    clustering = vec(rcopy(R"binder_clustering$cl"))
    R"rm(list = ls())"
    return clustering
end

function minbinder(markov_chain::Matrix{Int64}; num_burn_in=0)
    psm = compute_co_occurrence_matrix(markov_chain[:, (num_burn_in+1):end])
    return minbinder(psm, markov_chain)
end

function minVI(posterior_similarity_matrix, iterations)
    draws = Matrix(transpose(iterations))
    R"library('mcclust')"
    R"library('mcclust.ext')"
    @rput posterior_similarity_matrix
    @rput draws
    R"vi_clustering <- minVI(posterior_similarity_matrix, draws, method='draws')"
    clustering = vec(rcopy(R"vi_clustering$cl"))
    R"rm(list = ls())"
    return clustering
end

function minVI(markov_chain::Matrix{Int64}; num_burn_in=0)
    psm = compute_co_occurrence_matrix(markov_chain[:, (num_burn_in+1):end])
    return minVI(psm, markov_chain[:, (num_burn_in+1):end])
end

hasproperty(x, s::Symbol) = s in fieldnames(typeof(x))

function generate_all_clusterings(n::Int64)
    clusterings = Vector{Vector{Vector{Int64}}}(undef, 0)
    for i=1:n
        clusterings = vcat(clusterings, collect(partitions(1:n, i)))
    end
    return clusterings
end

end

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
function K_TMGT(elm::Tria)::Matrix{Float64}
    return K_GT9(elm)+K_TMT(elm)
end


text:
function triple02()
#
#  test polynomial suggested by Goedecker
#
    p = triple(10, 10, 10)
    z = [[0.9;1;1.1]  [10; 10; 10]]

    p, PolyZeros(z)
end


text:
# # Getting Started
#
# ## Installation
# The FinancialFFT package is available through the Julia package system
# by running `]add FinancialFFT`.
# Throughout, we assume that you have installed the package.

# ## Basic syntax
# The basic syntax for FinancialFFT is simple. Here is an example of pricing a European Option with a Black Scholes model:

using FinancialFFT, FinancialMonteCarlo
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;
# Define FinancialFFT Parameters:
Nsim = 10000;
Nstep = 30;
# Define Model Parameters.
σ = 0.2;
# Build the Lewis configuration and the zero rate.
method_lewis = LewisMethod(700.0, 200000);
rfCurve = ZeroRate(r);
# Define The Option
EuOption = EuropeanOption(T, K)
# Define the Model of the Underlying
Model = BlackScholesProcess(σ, Underlying(S0, d));
# Call the pricer metric
@show EuPrice = pricer(Model, rfCurve, method_lewis, EuOption);
# ## Using Other Models and Options
# The package contains a large number of models of three main types:
# * `Ito Process`
# * `Jump-Diffusion Process`
# * `Infinity Activity Process`

# The usage is analogous to the above example.


text:
#=
Code generation functions for Julia.
=#

# Extract the input args. This must match the arguments passed by the solver.
function handle_input_args_cg_julia(lorr, vors)
    code =
"var =       args[1];
eid =       args[2];
fid =       args[3];
grid =      args[4];
geo_facs =  args[5];
refel =     args[6]
time =      args[7];
dt =        args[8]
"
    # A trick for uniform grids to avoid repeated work
    if vors == "volume" && config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
        code *= "stiffness = args[9]; # set of stiffness matrices for each dimension\n";
        code *= "mass =      args[10]; # mass matrix\n";
    end

    return code;
end

# Allocate, compute, or fetch all needed values
function prepare_needed_values_cg_julia(entities, var, lorr, vors)
    # Only gather the pieces that are used. See the end of this function.
    if vors == "volume"
        # el, loc2glb, nodex, detj, J
        piece_needed = zeros(Bool, 5);
    else # surface
        # (els, loc2glb, nodex, frefelind, face2glb, facex, normal, fdetJ, vol_J)
        piece_needed = zeros(Bool, 9);
    end
    need_deriv_matrix = false;

    code = "";

    # First label the index values
    for i=1:length(indexers)
        if indexers[i].symbol === :elements || indexers[i].symbol === :cells
            # not needed
        else
            code *= "INDEX_VAL_"*string(indexers[i].symbol)*" = kwargs["*string(i)*"]\n";
        end
    end

    for i=1:length(entities)
        cname = make_coef_name(entities[i]);
        if is_test_function(entities[i])
            # Assign it a transpose quadrature matrix
            if vors == "volume"
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        code *= cname * " = TRQ"*string(entities[i].derivs[di])*"; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function\n";
                    end
                    need_deriv_matrix = true;
                else
                    code *= cname * " = refel.Q'; # test function.\n";
                end
            else
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        code *= cname * " = TRQ"*string(entities[i].derivs[di])*"; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function\n";
                        need_deriv_matrix = true;
                    end
                else
                    code *= cname * " = (refel.surf_Q[frefelind[1]]); # test function.\n";
                end
            end

        elseif is_unknown_var(entities[i], var) && lorr == LHS
            if vors == "volume"
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        code *= cname * " = RQ"*string(entities[i].derivs[di])*"; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function\n";
                    end
                    need_deriv_matrix = true;
                else
                    code *= cname * " = refel.Q; # trial function.\n";
                end
            else
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        code *= cname * " = RQ"*string(entities[i].derivs[di])*"; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function\n";
                        need_deriv_matrix = true;
                    end
                else
                    code *= cname * " = refel.surf_Q[frefelind[1]]; # trial function.\n";
                end
            end

        else
            # Is coefficient(number or function) or variable(array)?
            (ctype, cval) = get_coef_val(entities[i]);
            if ctype == -1
                # It was a special symbol like dt
                if vors == "surface"
                    if entities[i].name == "FACENORMAL1"
                        code *= cname * " = normal["*string(entities[i].index)*"]; # normal vector component\n"
                        piece_needed[7] = true;
                    else entities[i].name == "FACENORMAL2"
                        code *= cname * " = -normal["*string(entities[i].index)*"]; # reverse normal vector component\n"
                        piece_needed[7] = true;
                    end
                end
            elseif ctype == 0
                # It was a number, do nothing?
            elseif ctype == 1 # a constant wrapped in a coefficient
                # This generates something like: coef_k_i = 4;
                if length(entities[i].derivs) > 0
                    code *= cname * " = 0; # NOTE: derivative applied to constant coefficient = 0\n";
                else
                    code *= cname * " = " * string(cval) * ";\n";
                end

            elseif ctype == 2 # a coefficient function
                # This generates something like:
                ######################################
                # coef_n_i = zeros(refel.Np);
                # for coefi = 1:refel.Np
                #     coef_k_i[coefi] = (Femshop.genfunctions[cval]).func(x[1,coefi], x[2,coefi],x[3,coefi],time);
                # end
                ######################################
                # Note: for LHS, this is only Nfp values. For RHS, this is Np values, but only evaluated at face
                cargs = "(nodex[coefi], 0, 0, time)";
                if config.dimension == 2
                    cargs = "(nodex[1, coefi], nodex[2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargs = "(nodex[1, coefi], nodex[2, coefi], nodex[3, coefi], time)";
                end
                cargs1 = "(nodex[1][coefi], 0, 0, time)";
                if config.dimension == 2
                    cargs1 = "(nodex[1][1, coefi], nodex[1][2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargs1 = "(nodex[1][1, coefi], nodex[1][2, coefi], nodex[1][3, coefi], time)";
                end
                cargs2 = "(nodex[2][coefi], 0, 0, time)";
                if config.dimension == 2
                    cargs2 = "(nodex[2][1, coefi], nodex[2][2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargs2 = "(nodex[2][1, coefi], nodex[2][2, coefi], nodex[2][3, coefi], time)";
                end
                cargsf = "(facex[coefi], 0, 0, time)";
                if config.dimension == 2
                    cargsf = "(facex[1, coefi], facex[2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargsf = "(facex[1, coefi], facex[2, coefi], facex[3, coefi], time)";
                end
                if vors == "volume"
                    code *= cname * " = zeros(refel.Np);\n";
                    code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*" and interpolate at quadrature points.\n";
                        end
                        need_deriv_matrix = true;
                    else
                        code *= cname * " = refel.Q * " * cname * "; # Interpolate at quadrature points.\n";
                    end
                    piece_needed[3] = true;

                else # surface
                    # If derivatives are needed, must evaluate at all volume nodes.
                    if length(entities[i].derivs) > 0
                        code *= cname * " = zeros(refel.Np);\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.genfunctions["*string(cval)*"]).func" * cargs1 * " end\n";
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RD"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]]; # extract face values only.";
                        need_deriv_matrix = true;
                        piece_needed[3] = true;

                    else # no derivatives, only need surface nodes
                        code *= cname * " = zeros(Nfp);\n";
                        code *= "for coefi = 1:Nfp " * cname * "[coefi] = (Femshop.genfunctions["*string(cval)*"]).func" * cargsf * " end\n";
                        piece_needed[6] = true;
                    end

                    # Interpolate at surface quadrature points
                    code *= cname * " = refel.surf_Q[frefelind[1]][:,refel.face2local[frefelind[1]]] * " * cname * "; # Interpolate at quadrature points.\n";
                end

            elseif ctype == 3 # a known variable value
                # This generates something like: coef_u_1 = copy((Femshop.variables[1]).values[1, loc2glb])
                if typeof(entities[i].index) <: Array
                    # It is an indexed variable
                    if length(entities[i].index) == 1
                        indstr = "INDEX_VAL_"*entities[i].index[1];
                    else
                        # There is more than one index. Need to form an expression for it.
                        indstr = "(INDEX_VAL_"*entities[i].index[1];
                        indices = variables[cval].indexer;
                        for indi=2:length(entities[i].index)
                            indstr *= " + ("*string(length(indices[indi-1].range))*"*(INDEX_VAL_"*entities[i].index[indi]*"-1)";
                        end
                        for indi=1:length(entities[i].index)
                            indstr *= ")";
                        end
                    end

                else
                    indstr = string(entities[i].index);
                end

                if vors == "volume"
                    code *= cname * " = copy((Femshop.variables["*string(cval)*"]).values["*indstr*", loc2glb]);\n";
                    # Apply any needed derivative operators.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*"\n";
                        end
                        need_deriv_matrix = true;
                    else
                        code *= cname * " = refel.Q * " * cname * "; # Interpolate at quadrature points.\n";
                    end
                    piece_needed[2] = true;

                else # surface
                    # If derivatives are needed, must evaluate at all volume nodes.
                    if length(entities[i].derivs) > 0
                        code *= cname * " = copy((Femshop.variables["*string(cval)*"]).values["*indstr*", loc2glb[1]]);\n";
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RD"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        #     code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]; # extract face values only.";
                        need_deriv_matrix = true;
                        piece_needed[2] = true;

                    else # no derivatives
                        code *= cname * " = copy((Femshop.variables["*string(cval)*"]).values["*indstr*", loc2glb[1]]);\n";
                        # piece_needed[5] = true;
                        piece_needed[2] = true;
                    end

                    #     code *= cname * " = refel.surf_Q[frefelind[1]][:,refel.face2local[frefelind[1]]] * " * cname * "; # Interpolate at quadrature points.\n";
                end

            elseif ctype == 4 # an indexed coefficient
                # This generates something like:
                ######################################
                # coef_n_1 = zeros(refel.Np); # allocate
                # for coefi = 1:refel.Np
                #     coef_k_1[coefi] = (Femshop.coefficients[cval]).value[INDEX_VAL_i].func(x[1,coefi], x[2,coefi],x[3,coefi],time); # evaluate at nodes
                # end
                ######################################
                if vors == "surface" && length(entities[i].derivs) == 0
                    nodesymbol = "facex"
                    piece_needed[6] = true;
                else
                    nodesymbol = "nodex"
                    piece_needed[3] = true;
                end
                cargs = "("*nodesymbol*"[coefi], 0, 0, time)";
                if config.dimension == 2
                    cargs = "("*nodesymbol*"[1, coefi], "*nodesymbol*"[2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargs = "("*nodesymbol*"[1, coefi], "*nodesymbol*"[2, coefi], "*nodesymbol*"[3, coefi], time)";
                end

                indstr = "";
                for indi=1:length(entities[i].index)
                    if indi>1
                        indstr *= ",";
                    end
                    indstr *= "INDEX_VAL_"*entities[i].index[indi];
                end

                if vors == "volume"
                    code *= cname * " = zeros(refel.Np);\n";
                    code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*" and interpolate at quadrature points.\n";
                        end
                        need_deriv_matrix = true;
                    else
                        code *= cname * " = refel.Q * " * cname * "; # Interpolate at quadrature points.\n";
                    end

                else # surface
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        code *= cname * " = zeros(refel.Np);\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]]; # extract face values only.";
                        need_deriv_matrix = true;

                    else # no derivatives, only need surface nodes
                        code *= cname * " = zeros(refel.Nfp[frefelind[1]]);\n";
                        code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = (Femshop.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                    end
                    # Interpolate at surface quadrature points
                    code *= cname * " = refel.surf_Q[frefelind[1]][:,refel.face2local[frefelind[1]]] * " * cname * "; # Interpolate at quadrature points.\n";
                end

            end
        end # if coefficient
        #code *= "println(\"did \"*string("*string(i)*"))\n";
    end # entity loop

    # Add the needed pieces
    if vors == "volume"
        # el, loc2glb, nodex, detj, J
        piece_needed[4] = true; #I'm basically always going to need this.
        needed_pieces = "";
        if piece_needed[1]
            needed_pieces *= "el = eid;\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = grid.loc2glb[:,eid];    # local to global map\n"
        end
        if piece_needed[3]
            needed_pieces *= "nodex = grid.allnodes[:,loc2glb]; # node coordinates\n"
        end
        if piece_needed[4]
            needed_pieces *= "detj = geo_factors.detJ[eid];     # geometric factors\n"
            needed_pieces *= "wdetj = refel.wg .* detj;         # quadrature weights\n"
        end
        if piece_needed[5] || need_deriv_matrix
            needed_pieces *= "J = geo_factors.J[eid];           # geometric factors\n"
        end

        if need_deriv_matrix
            needed_pieces *=
"
# Note on derivative matrices:
# RQn are vandermond matrices for the derivatives of the basis functions
# with Jacobian factors. They are made like this.
# |RQ1|   | rx sx tx || Qx |
# |RQ2| = | ry sy ty || Qy |
# |RQ3|   | rz sz tz || Qz |

"
            use_full_deriv_mat = lorr==RHS;
            if config.dimension == 1
                needed_pieces *= "(RQ1, RD1) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "TRQ1 = RQ1';\n"

            elseif config.dimension == 2
                needed_pieces *= "(RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "(TRQ1, TRQ2) = (RQ1', RQ2');\n"

            elseif config.dimension == 3
                needed_pieces *= "(RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "(TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');\n"
            end
        end

    else # surface
        piece_needed[4] = true; #I'm basically always going to need this.
        piece_needed[8] = true; #I'm basically always going to need this.
        needed_pieces =
"if grid.face2element[1, fid] == eid # The normal points out of e"*(piece_needed[7] ? "\n\tnormal = grid.facenormals[:, fid];" : "")*"
    neighbor = grid.face2element[2, fid];
    frefelind = [grid.faceRefelInd[1,fid], grid.faceRefelInd[2,fid]]; # refel based index of face in both elements
else # The normal points into e"*(piece_needed[7] ? "\n\tnormal = -grid.facenormals[:, fid];" : "")*"
    neighbor = grid.face2element[1, fid];
    frefelind = [grid.faceRefelInd[2,fid], grid.faceRefelInd[1,fid]]; # refel based index of face in both elements
end
if neighbor == 0 # This is a boundary face. For now, just compute as if neighbor is identical. BCs handled later.
    neighbor = eid;
    frefelind[2] = frefelind[1];
end
Nfp = refel.Nfp[frefelind[1]]; # number of face nodes\n
";
        # (els, loc2glb, nodex, frefelind, face2glb, facex, normal, fdetJ, vol_J)
        if piece_needed[1]
            needed_pieces *= "els = (eid, neighbor); # indices of elements on both sides\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = (grid.loc2glb[:,eid], grid.loc2glb[:, neighbor]); # volume local to global\n"
        end
        if piece_needed[3]
            needed_pieces *= "nodex = (grid.allnodes[:,loc2glb[1][:]], grid.allnodes[:,loc2glb[2][:]]); # volume node coordinates\n"
        end
        if piece_needed[5] || piece_needed[6]
            needed_pieces *= "face2glb = grid.face2glb[:,:,fid];         # global index for face nodes for each side of each face\n"
        end
        if piece_needed[6]
            needed_pieces *= "facex = grid.allnodes[:, face2glb[:, 1]];  # face node coordinates\n"
        end
        if piece_needed[8]
            needed_pieces *= "face_detJ = geo_facs.face_detJ[fid];       # detJ on face\n"
            needed_pieces *= "wdetj = refel.surf_wg[1] .* face_detJ;     # quadrature weights\n"
        end
        if piece_needed[9] || need_deriv_matrix
            needed_pieces *= "vol_J = (geo_facs.J[eid], geo_facs.J[neighbor]);\n"
        end

        if need_deriv_matrix
            needed_pieces *=
"
# Note on derivative matrices:
# RQn are vandermond matrices for the derivatives of the basis functions
# with Jacobian factors. They are made like this.
# |RQ1|   | rx sx tx || Qx |
# |RQ2| = | ry sy ty || Qy |
# |RQ3|   | rz sz tz || Qz |

"
            use_full_deriv_mat = lorr==RHS;
            if config.dimension == 1
                needed_pieces *= "RQ1 = refel.surf_Qr[frefelind[1]] .* vol_J[1].rx[1]";
                needed_pieces *= "TRQ1 = RQ1'";

            elseif config.dimension == 2
                needed_pieces *= "(RQ1,RQ2,RD1,RD2) = build_face_deriv_matrix(refel, frefelind[1], vol_J[1], $use_full_deriv_mat)";
                needed_pieces *= "(TRQ1,TRQ1) = (RQ1',RQ1')";

            elseif config.dimension == 3
                needed_pieces *= "(RQ1,RQ2,RQ3,RD1,RD2,RD3) = build_face_deriv_matrix(refel, frefelind[1], vol_J[1], $use_full_deriv_mat)";
                needed_pieces *= "(TRQ1,TRQ2,TRQ3) = (RQ1',RQ2',RQ3')";
            end
        end
    end

    code = needed_pieces * "\n" * code;

    return code;
end

function make_elemental_computation_cg_julia(terms, var, dofsper, offset_ind, lorr, vors)
    # Here is where I make some assumption about the form of the expression.
    # Since it was expanded by the parser it should look like a series of terms: t1 + t2 + t3...
    # Where each term is multiplied by one test function component, and if LHS, involves one unknown component.
    # The submatrix modified by a term is determined by these, so go through the terms and divide them
    # into their submatrix expressions.
    # Each term will look something like
    # LHS: test_part * diagm(weight_part .* coef_part) * trial_part
    # RHS: test_part * (weight_part .* coef_part)
    code = "";

    # Allocate the vector or matrix to be returned if needed
    if dofsper > 1
        if lorr == RHS
            code *= "element_vector = zeros(refel.Np * "*string(dofsper)*"); # Allocate the returned vector.\n"
        else
            code *= "element_matrix = zeros(refel.Np * "*string(dofsper)*", refel.Np * "*string(dofsper)*"); # Allocate the returned matrix.\n"
        end
    end

    # Separate the factors of each term into test, trial, coef and form the calculation
    if dofsper > 1
        # Submatrices or subvectors for each component
        if lorr == LHS
            submatrices = Array{String, 2}(undef, dofsper, dofsper);
        else # RHS
            submatrices = Array{String, 1}(undef, dofsper);
        end
        for smi=1:length(submatrices)
            submatrices[smi] = "";
        end

        if typeof(var) <: Array
            for vi=1:length(var) # variables
                # Process the terms for this variable
                for ci=1:length(terms[vi]) # components
                    for i=1:length(terms[vi][ci])
                        (term_result, test_ind, trial_ind) = generate_term_calculation_cg_julia(terms[vi][ci][i], var, lorr, vors);

                        # println(terms)
                        # println(terms[vi])
                        # println(terms[vi][ci])
                        # println(terms[vi][ci][i])
                        # println(term_result * " : "*string(test_ind)*", "*string(trial_ind))

                        # Find the appropriate submatrix for this term
                        submati = offset_ind[vi] + test_ind;
                        submatj = trial_ind;
                        if lorr == LHS
                            submat_ind = submati + dofsper * (submatj-1);
                        else
                            submat_ind = submati;
                        end

                        if length(submatrices[submat_ind]) > 1
                            submatrices[submat_ind] *= " .+ " * term_result;
                        else
                            submatrices[submat_ind] = term_result;
                        end
                    end
                end

            end # vi

        else # only one variable
            # Process the terms for this variable
            for ci=1:length(terms) # components
                for i=1:length(terms[ci])
                    (term_result, test_ind, trial_ind) = generate_term_calculation_cg_julia(terms[ci][i], var, lorr, vors);

                    # Find the appropriate submatrix for this term
                    if lorr == LHS
                        submat_ind = test_ind + dofsper * (trial_ind-1);
                    else
                        submat_ind = test_ind;
                    end

                    if length(submatrices[submat_ind]) > 1
                        submatrices[submat_ind] *= " .+ " * term_result;
                    else
                        submatrices[submat_ind] = term_result;
                    end
                end
            end

        end

        # Put the submatrices together into element_matrix or element_vector
        if lorr == LHS
            for emi=1:dofsper
                for emj=1:dofsper
                    if length(submatrices[emi, emj]) > 1
                        rangei = "("*string(emi-1)*"*refel.Np + 1):("*string(emi)*"*refel.Np)";
                        rangej = "("*string(emj-1)*"*refel.Np + 1):("*string(emj)*"*refel.Np)";
                        code *= "element_matrix["*rangei*", "*rangej*"] = " * submatrices[emi,emj] * "\n";
                    end
                end
            end
            code *= "return element_matrix;\n"

        else # RHS
            for emi=1:dofsper
                if length(submatrices[emi]) > 1
                    rangei = "("*string(emi-1)*"*refel.Np + 1):("*string(emi)*"*refel.Np)";
                    code *= "element_vector["*rangei*"] = " * submatrices[emi] * "\n";
                end
            end
            code *= "return element_vector;\n"
        end

    else # one dof
        terms = terms[1];
        if lorr == LHS
            result = "zeros(refel.Np, refel.Np)";
        else
            result = "zeros(refel.Np)";
        end

        #process each term
        for i=1:length(terms)
            (term_result, test_ind, trial_ind) = generate_term_calculation_cg_julia(terms[i], var, lorr, vors);

            if i > 1
                result *= " .+ " * term_result;
            else
                result = term_result;
            end
        end
        code *= "return " * result * ";\n";
    end

    return code;
end

function generate_term_calculation_cg_julia(term, var, lorr, vors)
    precomputed_mass_stiffness = config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE;
    function extract_one_entity(ex)
        if typeof(ex) == Expr
            if ex.head == :call && (ex.args[1] == :- || ex.args[1] == :.-) && length(ex.args) == 2 && typeof(ex.args[2]) == SymEntity
                return (ex.args[2], true);
            else
                return nothing;
            end
        elseif typeof(ex) == SymEntity
            return (ex, false);
        end
        # return nothing
    end

    result = "";

    if lorr == LHS
        (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(term, var);
        # LHS: test_part * diagm(weight_part .* coef_part) * trial_part
        if !(coef_part === nothing)
            result = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                    string(replace_entities_with_symbols(coef_part)) * ") * " *
                    string(replace_entities_with_symbols(trial_part));
        else # no coef_part
            using_precomputed = false;
            if precomputed_mass_stiffness
                test_ent= extract_one_entity(test_part); # defined above
                trial_ent = extract_one_entity(trial_part);
                if !(test_ent === nothing || trial_ent === nothing)
                    if length(test_ent[1].derivs) == 0 && length(trial_ent[1].derivs) == 0
                        if test_ent[2] || trial_ent[2]
                            result = "-mass";
                        else
                            result = "mass";
                        end
                        using_precomputed = true;
                    elseif length(test_ent[1].derivs) == 1 && length(trial_ent[1].derivs) == 1 && test_ent[1].derivs[1] == trial_ent[1].derivs[1]
                        if test_ent[2] || trial_ent[2]
                            result = "-stiffness["*string(test_ent[1].derivs[1])*"]";
                        else
                            result = "stiffness["*string(test_ent[1].derivs[1])*"]";
                        end
                        using_precomputed = true;
                    end
                end
            end
            if !using_precomputed
                result = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                    string(replace_entities_with_symbols(trial_part));
            end
        end
    else
        (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(term);
        # RHS: test_part * (weight_part .* coef_part)
        if !(coef_part === nothing)
            result = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* " *
                    string(replace_entities_with_symbols(coef_part)) * ")";
        else
            result = string(replace_entities_with_symbols(test_part)) * " * (wdetj)";
        end
    end

    return (result, test_ind, trial_ind);
end

# Generate the assembly loop structures and insert the content
function generate_assembly_loop_cg_julia(indices)
    # Each of the indices must be passed to the functions in a named tuple.
    # Pass all defined indexers.
    index_args = "";
    for i=1:length(indexers)
        if i > 1
            index_args *= ", ";
        end
        index_args *= "index_val_"*string(indexers[i].symbol)*"=indexing_variable_"*string(indexers[i].symbol);
    end
    if length(indexers) == 0
        index_args = "empty_index = 1";
    end

    # Make sure the elemental loop is included
    elements_included = false;
    for i=1:length(indices)
        if indices[i] == "elements" || indices[i] == "cells"
            elements_included = true;
        end
    end
    # If elements were not included, make them the outermost loop
    if !elements_included
        indices = ["elements"; indices];
    end

    code =
"# Label things that were allocated externally
b = allocated_vecs[1];
if !rhs_only
    AI = allocated_vecs[2];
    AJ = allocated_vecs[3];
    AV = allocated_vecs[4];
end
# zero b
b .= 0;

Np = refel.Np;
nel = mesh_data.nel;

# Stiffness and mass are precomputed for uniform grid meshes
precomputed_mass_stiffness = config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
if precomputed_mass_stiffness
    wdetj = refel.wg .* geo_factors.detJ[1];
    J = geo_factors.J[1];

    if config.dimension == 1
        (RQ1, RD1) = build_deriv_matrix(refel, J);
        TRQ1 = RQ1';
        stiffness = [(TRQ1 * diagm(wdetj) * RQ1)];
    elseif config.dimension == 2
        (RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);
        (TRQ1, TRQ2) = (RQ1', RQ2');
        stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2)];
    else
        (RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);
        (TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');
        stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2) , (TRQ3 * diagm(wdetj) * RQ3)];
    end
    mass = (refel.Q)' * diagm(wdetj) * refel.Q;
else
    stiffness = 0;
    mass = 0;
end

loop_time = Base.Libc.time();

"
    # generate the loop structures
    code *= "# Loops\n"
    loop_start = ""
    loop_end = ""
    ind_offset = ""
    prev_ind = ""
    for i=1:length(indices)
        if indices[i] == "elements" || indices[i] == "cells"
            loop_start *= "for eid in elemental_order\n";
            # loop_start *= "face_done .= 0; # Reset face_done in elemental loop\n";
            loop_end *= "end # loop for elements\n";
        else
            loop_start *= "for indexing_variable_"*string(indices[i].symbol)*" in "*string(indices[i].range)*"\n";
            loop_end *= "end # loop for "*string(indices[i].symbol)*"\n";
            if length(prev_ind) == 0
                ind_offset *= "(indexing_variable_"*string(indices[i].symbol);
                prev_ind = string(length(indices[i].range));
            else
                ind_offset *= " + "*prev_ind*" * (indexing_variable_"*string(indices[i].symbol)*" - 1";
                prev_ind = string(length(indices[i].range));
            end
        end
    end
    for i=2:length(indices)
        ind_offset *= ")";
    end
    # index_offset = "dofs_per_loop * ("*ind_offset*" - 1 + dofs_per_node * (eid - 1)) + 1"
    # face_index_offset = "dofs_per_loop * ("*ind_offset*" - 1 + dofs_per_node * (fid - 1)) + 1"

    code *= loop_start;

    # insert the content
    code *= "
    loc2glb = grid_data.loc2glb[:,eid]; # global indices of this element's nodes

    # determine index in global vector
    dof_index = (loc2glb .- 1) .* dofs_per_node .+ "*ind_offset*";
    #println(\"dof index \"*string(dof_index))

    if !rhs_only
        Astart = (eid-1)*Np*Np*dofs_per_node + 1; # The segment of AI, AJ, AV for this element
        #println(\"Astart \"*string(Astart))
    end

    volargs = (var, eid, 0, grid_data, geo_factors, refel, t, dt, stiffness, mass);
    linchunk = linear.func(volargs; "*index_args*");  # get the elemental linear part
    if !rhs_only
        bilinchunk = bilinear.func(volargs; "*index_args*"); # the elemental bilinear part
    end

    if dofs_per_node == 1
        b[loc2glb] += linchunk;
        if !rhs_only
            for jj=1:Np
                offset = Astart - 1 + (jj-1)*Np;
                for ii=1:Np
                    AI[offset + ii] = loc2glb[ii];
                    AJ[offset + ii] = loc2glb[jj];
                    AV[offset + ii] = bilinchunk[ii, jj];
                end
            end
        end

    else # more than one DOF per node
        b[dof_index] += linchunk;
        if !rhs_only
            #print(\"offsets \")
            for jj=1:Np
                offset = Astart + Np*Np*("*ind_offset*"-1) + (jj-1)*Np - 1;
                #print(\", \"*string(offset))
                for ii=1:Np
                    AI[offset + ii] = dof_index[ii];
                    AJ[offset + ii] = dof_index[jj];
                    AV[offset + ii] = bilinchunk[ii,jj];
                end
            end
            #println()
        end
    end
"
    # close loops
    code *= loop_end;

    # finish
    code *=
"
loop_time = Base.Libc.time() - loop_time;

if !rhs_only
    # Build the sparse A. Uses default + to combine overlaps
    #println(\"AI \"*string(AI))
    #println(\"AJ \"*string(AJ))
    #println(\"AV \"*string(AV))
    A = sparse(AI[1:nel*Np*Np*dofs_per_node], AJ[1:nel*Np*Np*dofs_per_node], AV[1:nel*Np*Np*dofs_per_node]);
end

# Boundary conditions
bc_time = Base.Libc.time();
if rhs_only
    b = CGSolver.apply_boundary_conditions_rhs_only(var, b, t);
else
    (A, b) = CGSolver.apply_boundary_conditions_lhs_rhs(var, A, b, t);
end
bc_time = Base.Libc.time() - bc_time;

if rhs_only
    return b;

else
    log_entry(\"Elemental loop time:     \"*string(loop_time));
    log_entry(\"Boundary condition time: \"*string(bc_time));
    #println(size(A))
    #println(size(b))
    #display(Array(A))
    #println(b)
    return (A, b);
end
"
    # println(code)
    return code;
end


text:
facts("Random sampling on unit, n-dimensional sphere") do

  uv1 = BlackBoxOptim.sample_unit_hypersphere(1)
  @fact size(uv1, 1) => 1

  uv2 = BlackBoxOptim.sample_unit_hypersphere(2)
  @fact size(uv2, 1) => 2

  uv3 = BlackBoxOptim.sample_unit_hypersphere(3, 4)
  @fact size(uv3, 1) => 3
  @fact size(uv3, 2) => 4

  for(i in 1:100)
    n = rand(1:100)
    num = rand(2:10)
    u = BlackBoxOptim.sample_unit_hypersphere(n, num)
    @fact size(u, 1) => n
    @fact size(u, 2) => num
    for(j in 1:num)
      @fact isapprox( norm(u[:,j]), 1.0 ) => true
    end
  end

end

facts("Random direction generator") do

  rdg1 = BlackBoxOptim.RandomDirectionGen(2, 3)
  ds1 = BlackBoxOptim.directions_for_k(rdg1, 1)
  @fact size(ds1) => (2, 3)

  rdg2 = BlackBoxOptim.RandomDirectionGen(10, 17)
  ds2 = BlackBoxOptim.directions_for_k(rdg2, 1)
  @fact size(ds2) => (10, 17)

end

facts("Mirrored random direction generator") do

  mrdg1 = BlackBoxOptim.MirroredRandomDirectionGen(2, 4)
  ds1 = BlackBoxOptim.directions_for_k(mrdg1, 1)
  @fact size(ds1) => (2, 4)
  @fact ds1[:,3] == -ds1[:,1] => true
  @fact ds1[:,4] == -ds1[:,2] => true

  mrdg2 = BlackBoxOptim.MirroredRandomDirectionGen(10, 6)
  ds2 = BlackBoxOptim.directions_for_k(mrdg2, 1)
  @fact size(ds2) => (10, 6)
  @fact ds2[:,4] == -ds2[:,1] => true
  @fact ds2[:,5] == -ds2[:,2] => true
  @fact ds2[:,6] == -ds2[:,3] => true

  # Must be even number of directions
  @fact_throws BlackBoxOptim.MirroredRandomDirectionGen(10, 1)
  @fact_throws BlackBoxOptim.MirroredRandomDirectionGen(10, 3)
  @fact_throws BlackBoxOptim.MirroredRandomDirectionGen(10, 7)

end


text:
# ----------------------------------------------------------------------------
# Sep24 | more issues with lists and sub
md = """
    1. A
    1. B
    """
g = FP.md_partition(md) |> FP.md_grouper
@test g[1].role == :LIST
@test g[1] // md
md = """
    * A **B** `C`
    * D _E_
    """
mdss = FP.subs(md, 3, 13)
p = FP.md_partition(mdss) |> FP.md_grouper
@test p[1].ss // "A **B** `C`"
# ----------------------------------------------------------------------------
# Sep21 | Empty first block (P_BREAK)
md = """

    * A
    """
mdss = FP.subs(md, 4, 5)
p = FP.md_partition(mdss)
@test p[1].name == :TEXT
@test p[1] // "A"


text:
using LinearAlgebra
using DoubleFloats
# using SemiLagMPI
using MPI
using SemiLagrangian

function printout(advd::AdvectionData{T,N,timeopt}, str) where {T,N,timeopt}
    if timeopt != MPIOpt || advd.adv.mpid.ind == 1
        println(str)
    end
end
printout(str) = println(str)

function sqg2(
    t_max::T,
    timeopt,
    dt::T,
    sz,
    interp::AbstractInterpolation,
    typealg::TimeAlgorithm,
    ordalg::Int,
) where {T}
    nbdt = Int(round(t_max / dt))
    spmin, spmax, nsp = T(0), T(1e6), sz[1]
    vmin, vmax, nv = T(0), T(1e6), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    tabst = [([1, 2], 2, 1, false)]

    adv = Advection(
        (mesh_sp, mesh_v),
        [interp, interp],
        dt,
        tabst,
        tab_coef = nosplit(dt),
        timeopt = timeopt,
        timealg = typealg,
        ordalg = ordalg,
    )

    pvar = getgeovar(adv)

    result = Array{T,2}[]
    initdatas = missing
    if typealg == ABTimeAlg_init
        nbfordt = Int(floor(sqrt(nbdt))) + ordalg + 2
        newnb = nbfordt * (3ordalg - 1)
        newdt = dt / nbfordt
        newt_max = dt * (3ordalg - 1)
        newtypealg = ordalg == 2 ? NoTimeAlg : typealg
        newordalg = ordalg == 2 ? 0 : (ordalg - 1)
        res = sqg2(newt_max, timeopt, newdt, sz, interp, newtypealg, newordalg)
        initdatas = res[nbfordt:nbfordt:end]
        append!(result, initdatas)
    end

    advd = AdvectionData(adv, zeros(T, sz), pvar, initdatas = initdatas)

    initdata!(pvar, advd)

    borne_t = t_max - dt / 2

    while advd.time_cur < borne_t
        while advection!(advd)
        end
        push!(result, copy(advd.data))
    end
    return result
end

function run_mesure(t_max::T, timeopt, sz, interp) where {T}

    tabtypealg = [
        NoTimeAlg,
        ABTimeAlg_ip,
        ABTimeAlg_init,
        ABTimeAlg_init,
        ABTimeAlg_init,
        ABTimeAlg_init,
    ]
    # ta

    tabordalg = [0, 2, 2, 3, 4, 5]
    # tabtxtsplit = ["stdsplit", "strangsplit", "triplejumpsplit", "order6split", "fernandosplit"]
    # tabtxt = ["stdsplit", "strangsplit", "fernandosplit", "std2d", "stdAB_2", "stdAB_3", "stdAB_4", "stdAB_5", "stdAB_6",]
    # tabtxt = ["stdsplit", "strangsplit", "fernandosplit", "stdRK4",]
    tabtxt = ["NO", "AB_ip2", "AB_init2", "AB_init3", "AB_init4", "AB_init5"]
    tabnbdt = [
        50,
        100,
        200,
        500,
        1000,
        2000,
        5000,
        10000,
        20000,
        50000,
        100000,
        200000,
        500000,
        1000000,
    ]

    #        if MPI.Comm_rank(MPI.COMM_WORLD) == 1
    #        end

    resdata = [zeros(T, sz) for i in CartesianIndices((length(tabnbdt), length(tabtxt)))]

    lastind = zeros(Int, length(tabtxt))
    t_loc = time_ns()
    for inbdt = 1:length(tabnbdt), itc = 1:length(tabtxt)
        lastind[itc] = inbdt
        tabres = sqg2(
            t_max,
            timeopt,
            t_max / tabnbdt[inbdt],
            sz,
            interp,
            tabtypealg[itc],
            tabordalg[itc],
        )
        resdata[inbdt, itc] .= tabres[end]
        # to run with mpi here is the command for n process (n must be a number) from SemiLagrangian.jl :
        # ./mpiexec.loc n notebooks/run_mesure_sqg.jl
        # else run like this :
        # julia --project=.  notebooks/run_mesure_sqg.jl

        if timeopt != MPIOpt || MPI.Comm_rank(MPI.COMM_WORLD) == 1
            nb = timeopt != MPIOpt ? (timeopt == NoTimeOpt ? 1 : Threads.nthreads()) : MPI.Comm_size(MPI.COMM_WORLD)
            t_now = time_ns()
            t = (t_now - t_loc)*1e-9
            t_loc = t_now
            for k = 1:length(tabtxt)
                if lastind[k] > 0 && (k != 1 || lastind[k] > 1)
                    println(
                        "# sz=$sz t_max=$t_max interp=$interp ref=$(tabtxt[k]) timeopt=$timeopt nb=$nb t=$t\n# t",
                    )
                    for txt in tabtxt
                        print("\t$txt")
                    end
                    println("")
                    for j = 1:length(tabnbdt), i = 0:length(tabtxt)
                        if i == 0
                            res = Float32(t_max / tabnbdt[j])
                        else
                            res = Float64(
                                (j < lastind[k] || (j == lastind[k] && i < k)) ?
                                norm(resdata[j, i] - resdata[lastind[k], k]) : 0,
                            )
                        end
                        print("$res")
                        if i == length(tabtxt)
                            print("\n")
                        else
                            print("\t")
                        end
                    end
                end
            end

            println("free memory : $(Sys.free_memory()/2^30)")
            flush(stdout)
        end
    end
end

T = Double64
# run_mesure(T(100000), MPIOpt, (128, 128), Lagrange(11, T))
run_mesure(T(100000), SimpleThreadsOpt, (128, 128), Lagrange(11, T))


text:
module RandomBasedArrays

export RandomBasedArray

struct RandomBasedArray{T,N,P<:AbstractArray} <: AbstractArray{T,N}
    parent::P
end
RandomBasedArray(p::P) where {T,N,P<:AbstractArray{T,N}} =
    RandomBasedArray{T,N,P}(p)

index(p::AbstractArray) = rand(eachindex(p))

Base.parent(A::RandomBasedArray) = A.parent
Base.size(A::RandomBasedArray) = size(parent(A))
Base.getindex(A::RandomBasedArray, ::Int) =
    @inbounds getindex(parent(A), index(parent(A)))
Base.getindex(A::RandomBasedArray{T,N}, ::Vararg{Int,N}) where {T,N} =
    getindex(parent(A), Base._ind2sub(parent(A), index(parent(A)))...)
Base.setindex!(A::RandomBasedArray, v, ::Int) =
    @inbounds setindex!(parent(A), v, index(parent(A)))
Base.setindex!(A::RandomBasedArray{T,N}, v, ::Vararg{Int,N}) where {T,N} =
    setindex!(parent(A), v, Base._ind2sub(parent(A), index(parent(A)))...)
# Show a fixed copy of the array, otherwise the alignment algorithm would get
# crazy trying to show a "mutating" array
Base.show(io::IO, m::MIME"text/plain", A::RandomBasedArray) = show(io, m, copy(A))

end # module


text:
function set_api_key(api_key)

    if "API_FRED" in keys(ENV)
        if length(ENV["API_FRED"]) == 0
            ENV["API_FRED"] = api_key
            println("Your API key is: " * ENV["API_FRED"])
        end
    else
        ENV["API_FRED"] = api_key
        println("local API key is set.")
    end

    return(nothing)
end


text:
"""
    LinearKoopman(K, B, Q, P, discrete)

An approximation of the koopman operator which is linear in the states.

`K` is the array representing the operator, `B` is the (possible present) array
representing the influence of exogenuos inputs on the evolution.

`Q` and `P` are matrices used for updating the operator with new measurements.
`discrete` indicates if the operator is discrete or continuous.

The koopman operator is callable with the typical signature of `f(u,p,t)` and `f(du,u,p,t)` respectively.

# Example

```julia
k = LinearKoopman([1.0 0; 0 0.3], [], [], [], true)

u = k([2.0; 0.5], nothing, nothing)
du = similiar(u)
k(du, u, nothing, nothing)
```
"""
mutable struct LinearKoopman <: AbstractKoopmanOperator
    operator::AbstractArray
    input::AbstractArray

    Q::AbstractArray
    P::AbstractArray

    discrete::Bool
end

outputmap(k::LinearKoopman) = throw(AssertionError("Linear Koopman Operator has no output map."))

function (k::LinearKoopman)(u, p, t)
    return k.operator*u
end

function (k::LinearKoopman)(du, u, p, t)
    mul!(du, k.operator, u)
end

"""
    update!(k, X, Y; threshold = eps())

Update the Koopman `k` given new data `X` and `Y`. The operator is updated in place, if
the L2 error of the prediction exceeds the `threshold`.

"""
function update!(k::LinearKoopman, X::AbstractArray, Y::AbstractArray; threshold::T = eps()) where {T <: Real}
    @assert updateable(k) "Linear Koopman is not updateable."

    ϵ = norm(Y-operator(k)*X, 2)

    if ϵ < threshold
        return
    end

    k.Q += Y*X'
    k.P += X*X'
    k.operator .= k.Q*inv(k.P)
    return
end


text:
function getsites(g)
  x = filter(x -> x isa Tuple, g)
  isempty(x) && return x
  return only(x)
end

sort_gates_by(g) =
  TupleTools.sort(getsites(g))

function sort_gates_lt(g1, g2)
  if length(g1) ≠ length(g2)
    return length(g1) > length(g2)
  end
  return g1 < g2
end

sort_gates(gates) =
  sort(gates; by=sort_gates_by, lt=sort_gates_lt)

function trotter1(H::OpSum, δτ::Number)
  onequbitgates = Tuple[]
  multiqubitgates = Tuple[]

  n = 1
  for k in 1:length(H)
    coupling = ITensors.coef(H[k])
    O = ITensors.ops(H[k])
    length(O) > 1 && error("only a single operator allowed per term")
    localop = ITensors.name(O[1])
    support = ITensors.sites(O[1])
    params = ITensors.params(O[1])

    # single-qubit gate
    if length(support) == 1
      g = (localop, support[1], (params..., f = x -> exp(-δτ * coupling * x),))
      push!(onequbitgates, g)
      n = support[1] ≥ n ? support[1] : n
    # multi-qubit gate
    else
      g = (localop, support, (params..., f = x -> exp(-δτ * coupling * x),))
      push!(multiqubitgates, g)
      n = maximum(support) ≥ n ? maximum(support) : n
    end
  end

  sorted_multi_qubit = sort_gates(multiqubitgates)
  sorted_one_qubit = onequbitgates[sortperm([s[2] for s in onequbitgates])]

  # TODO: simplify this unison sorting loop
  tebd1 = Tuple[]
  g1_counter = 1; nmq = length(sorted_multi_qubit)
  gm_counter = 1; n1q = length(onequbitgates)
  for j in 1:n
    while (gm_counter ≤ nmq) && any(sorted_multi_qubit[gm_counter][2] .== j)
      push!(tebd1, sorted_multi_qubit[gm_counter])
      gm_counter += 1
    end
    while (g1_counter ≤ n1q) && any(sorted_one_qubit[g1_counter][2] .== j)
      push!(tebd1, sorted_one_qubit[g1_counter])
      g1_counter += 1
    end
  end
  #tebd1 = vcat(sorted_multi_qubit, sorted_one_qubit)
  return tebd1
end

"""
    trotter2(H::OpSum; δt::Float64=0.1, δτ=im*δt)

Generate a single layer of gates for one step of 2nd order TEBD.
"""
function trotter2(H::OpSum, δτ::Number)
  tebd1 = trotter1(H, δτ/2)
  tebd2 = copy(tebd1)
  append!(tebd2, reverse(tebd1))
  return tebd2
end

function trotter4(H::OpSum, δτ::Number)
  δτ1 = δτ / (4 - 4^(1/3))
  δτ2 = δτ - 4 * δτ1

  tebd2_δ1 = trotter2(H, δτ1)
  tebd2_δ2 = trotter2(H, δτ2)

  tebd4 = vcat(tebd2_δ1,tebd2_δ1)
  append!(tebd4, tebd2_δ2)
  append!(tebd4, vcat(tebd2_δ1,tebd2_δ1))
  return tebd4
end

"""
    trotterlayer(H::OpSum; order::Int = 2, kwargs...)

Generate a single layer of gates for one step of TEBD.
"""
function trotterlayer(H::OpSum, δτ::Number; order::Int = 2)
  order == 1 && return trotter1(H, δτ)
  order == 2 && return trotter2(H, δτ)
  error("Automated Trotter circuits with order > 2 not yet implemented")
  # TODO: understand weird behaviour of trotter4
  #order == 4 && return trotter4(H, δτ)
  #error("Automated Trotter circuits with order > 2 not yet implemented")
end

trottercircuit(H; kwargs...) =
  _trottercircuit(H, get_times(;kwargs...); kwargs...)

function _trottercircuit(H::Vector{<:OpSum}, τs::Vector; order::Int = 2, layered::Bool = true, kwargs...)
  @assert length(H) == (length(τs) -1) || length(H) == length(τs)
  δτs = diff(τs)
  circuit = [trotterlayer(H[t], δτs[t]; order = order) for t in 1:length(δτs)]
  layered && return circuit
  return vcat(circuit...)
end

_trottercircuit(H::OpSum, τs::Vector; kwargs...) =
  _trottercircuit(repeat([H], length(τs)-1), τs; kwargs...)

get_times(; δt=nothing, δτ=nothing, t=nothing, τ=nothing, ts=nothing, τs=nothing, kwargs...) = get_times(δt, δτ, t, τ, ts, τs)

get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::Vector,  τs::Nothing)   = im .* ts
get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::Nothing, τs::Vector)    = τs

function get_times(δt::Nothing, δτ::Number, t::Nothing, τ::Number, ts::Nothing, τs::Nothing)
  depth = abs(τ / δτ)
  (depth-Int(floor(depth)) > 1e-5) && @warn "Incommensurate Trotter step!"
  return collect(0.0:δτ:τ)
end

get_times(δt::Number, δτ::Nothing, t::Number, τ::Nothing, ts::Nothing, τs::Nothing) =
  im .* get_times(δτ, δt, τ, t, ts, τs)

get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::AbstractRange, τs::Nothing)       =
  get_times(δt, δτ, t, τ, collect(ts), τs)

get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::Nothing,       τs::AbstractRange) =
  get_times(δt, δτ, t, τ, ts, collect(τs))

get_times(δt::Number,  δτ::Nothing, t::Nothing, τ::Nothing, ts::Nothing, τs::Nothing)   =
  error("Total time `t` not set.")

get_times(δt::Nothing, δτ::Number,  t::Nothing, τ::Nothing, ts::Nothing, τs::Nothing)   =
  error("Total imaginary time `τ` not set")

get_times(δt::Nothing,  δτ::Nothing, t::Number, τ::Nothing, ts::Nothing, τs::Nothing)   =
  error("Trotter step `δt` not set.")

get_times(δt::Nothing,  δτ::Nothing, t::Nothing, τ::Number, ts::Nothing, τs::Nothing)   =
  error("Imaginary Trotter step `δτ` not set.")

#get_times(; δt=nothing, δτ=nothing, t=nothing, τ=nothing, ts=nothing, τs=nothing, kwargs...) = get_times(δt, δτ, t, τ, ts, τs)
#
#get_times(δt::Number,  δτ::Nothing, t::Nothing, τ::Nothing, ts::Nothing, τs::Nothing)   = im * δt
#get_times(δt::Nothing, δτ::Number,  t::Nothing, τ::Nothing, ts::Nothing, τs::Nothing)   = δτ
#get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::Vector,  τs::Nothing)   = im .* ts
#get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::Nothing, τs::Vector)    = τs
#
#function get_times(δt::Nothing, δτ::Number, t::Nothing, τ::Number, ts::Nothing, τs::Nothing)
#  depth = abs(τ / δτ)
#  (depth-Int(floor(depth)) > 1e-5) && @warn "Incommensurate Trotter step!"
#  return collect(0.0:δτ:τ)
#end
#
#get_times(δt::Number, δτ::Nothing, t::Number, τ::Nothing, ts::Nothing, τs::Nothing) =
#  im .* get_times(δτ, δt, τ, t, ts, τs)
#
#get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::AbstractRange, τs::Nothing)       = get_times(δt, δτ, t, τ, collect(ts), τs)
#get_times(δt::Nothing, δτ::Nothing, t::Nothing, τ::Nothing, ts::Nothing,       τs::AbstractRange) = get_times(δt, δτ, t, τ, ts, collect(τs))
#
#