text:
@testset "ClosestVector.jl" begin
    #=
    The quadratic triple QT = [Q, L, c] is an n-dimensional ellipsoid with rational coefficients.
    Where,
    Q is an nxn symmetric matrix of the quadratic form ,
    L is a rational column vector of length n &
    c is a rational number
    =#
    Q = matrix(QQ,3,3,[1 0 0; 0 1 0; 0 0 1]);
    L = matrix(QQ,3,1,[1,0,0]);
    c = fmpq(2//5);
    Lat = Zlattice(gram = matrix(QQ,3,3,[1,0,0,0,1,0,0,0,1]));
    v = convert(Array{Union{RingElem, AbstractFloat, Integer, Rational},1},[-1,0,0]);
    u = fmpq(3//5);
    @test Hecke.closest_vectors(Q,L,c)[1] == Hecke.closest_vectors(Lat, v, u)[1]

    function compare_functions(Q::MatrixElem, K::MatrixElem, d::RingElement)
        List1 = Hecke.closest_vectors(Q, K, d);
        L, v, c = Hecke.convert_type(Q, K, d)
        List2 = Hecke.closest_vectors(L, v, c);
        for i in 1:length(List1)
            if List1[i] == List2[i]
                return true
            else
                error("$i-th entry of the two lists don't match.")
            end
        end
    end

    function compare_functions(L::ZLat, v::Vector{RingElement} , upperbound::RingElement)
        List1 = Hecke.closest_vectors(L, v, upperbound);
        Q, K, d = Hecke.convert_type(L, v, upperbound);
        List2 = Hecke.closest_vectors(Q, K, d);
        for i in 1:length(List1)
            if List1[i] == List2[i]
                return true
            else
                error("$i-th entry of the two lists don't match.")
            end
        end
    end

    #EXAMPLE 1a: 3-dimensional sphere
    Q1 = matrix(QQ,3,3,[1 0 0; 0 1 0; 0 0 1]);
    L1 = matrix(QQ,3,1,[1,1,1]);
    c1 = 1;

    #EXAMPLE 1b: 4-dimensional sphere
    Q2 = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
    L2 = matrix(QQ, 4, 1 ,[1,1,1,1]);
    c2 = 3;

    #EXAMPLE 1c: 6-dimensional sphere
    Q3 = matrix(QQ,6,6,[1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1]);
    L3 = matrix(QQ,6,1,[1,0,1,0,2,0]);
    c3 = fmpq(4//3);

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

    #Examples for closest_vectors()
    #Example 1:
    Lat1 = Zlattice(gram = matrix(QQ,3,3,[1,0,0,0,1,0,0,0,1]));
    v1 = convert(Array{Union{RingElem, AbstractFloat, Integer, Rational},1},[-1,0,0]);
    u1 = fmpq(3//5);

    #Example 1:
    Lat2 = Zlattice(gram = matrix(QQ,4,4,[1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1]));
    v2 = convert(Array{Union{RingElem, AbstractFloat, Integer, Rational},1},[-1,-1,-1,-1]);
    u2 = fmpq(41//11);

    #Example 1:
    Lat3 = Zlattice(gram = matrix(QQ,6,6,[1,0,0,0,0,0, 0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1]));
    v3 = convert(Array{Union{RingElem, AbstractFloat, Integer, Rational},1},[-1,0,-1,0,-2,0]);
    u3 = fmpq(14//3);

    @testset "Comparing closest_vectors(Q, K, d) to closest_vectors(Lat, vec, upperbound) list" begin
        @test compare_functions(Q1, L1, c1) == true
        @test compare_functions(Q2, L2, c2) == true
        @test compare_functions(Q3, L3, c3) == true
    end

    @testset "Comparing closest_vectors(Lat, vec, upperbound) list to closest_vectors(Q, K, d)" begin
        @test compare_functions(Lat1, v1, u1) == true
        @test compare_functions(Lat2, v2, u2) == true
        @test compare_functions(Lat3, v3, u3) == true
    end

    @testset "Conversion from Quadratic Triple set to Closest Vector Enumeration input parameters" begin
        Q3 = matrix(QQ,6,6,[1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1]);
        L3 = matrix(QQ,6,1,[1,0,1,0,2,0]);
        c3 = fmpq(4//3);
        cvp = Hecke.convert_type(Q3,L3,c3);
        @test  det(basis_matrix(cvp[1])) == 1
        @test cvp[2] == [-1, 0, -1, 0, -2, 0]
        @test cvp[3] == 14//3
    end

    @testset "Conversion from Closest Vector Enumeration to Quadratic Triple set input parameters" begin
        Lat = Zlattice(gram = matrix(QQ,3,3,[1,0,0,0,1,0,0,0,1]));
        vec = convert(Array{Union{RingElem, AbstractFloat, Integer, Rational},1},[-1,-1,-1]);
        upbound = 2;
        qt = Hecke.convert_type(Lat, vec, upbound);
        @test qt[1] == matrix(QQ,3,3,[1 0 0; 0 1 0; 0 0 1])
        @test qt[2] == matrix(QQ,3,1,[1,1,1])
        @test qt[3] == 1
    end

  @testset "Quadratic Triple Tests" begin
      @testset "Projective Quadratic Triple" begin
            Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
            L = matrix(QQ, 4, 1 ,[1,1,1,1]);
            c = 3;
            @test Hecke.pojective_quadratic_triple(Q, L, c, 1) == (Q[1:1, 1:1], L[1:1,1:1], 0)
            @test Hecke.pojective_quadratic_triple(Q, L, c, 4) == (Q, L, c)
        end

        @testset "Finite set of integers for a 1-dimensional quadratic triple" begin
            Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
            L = matrix(QQ, 4, 1 ,[1,1,1,1]);
            c = 3;
            Q1 = Hecke.pojective_quadratic_triple(Q, L, c, 1)[1];
            L1 = Hecke.pojective_quadratic_triple(Q, L, c, 1)[2];
            c1 = Hecke.pojective_quadratic_triple(Q, L, c, 1)[3];
            @test Hecke.range_ellipsoid_dim1(Q1, L1, c1) == -2:0
        end

        @testset "Different Types of Extended Quadratic Triple" begin
            @testset "Quadratic triple extension from a 1-dimensional quadratic triple" begin
                Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
                L = matrix(QQ, 4, 1 ,[1,1,1,1]);
                c = 3;
                @test Hecke.positive_quadratic_triple(fmpz(1), Q, L, c)[1] == Q[2:4, 2:4]
                @test Hecke.positive_quadratic_triple(fmpz(1), Q, L, c)[2] == L[2:4,1:1]
                @test Hecke.positive_quadratic_triple(fmpz(1), Q, L, c)[3] == 6
            end

            @testset "Quadratic triple extension from a m-dimensional quadratic triple to m+n-1 quadratic triple" begin
                Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
                L = matrix(QQ, 4, 1 ,[1,1,1,1]);
                c = 3;
                @test Hecke.positive_quadratic_triple2(fmpz[-2,0,1], Q,L,c)[1][1] == 1
                @test Hecke.positive_quadratic_triple2(fmpz[-2,0,1], Q,L,c)[2][1] == 1
                @test Hecke.positive_quadratic_triple2(fmpz[-2,0,1], Q,L,c)[3] == 6
            end
        end
    end

    @testset "Finite set tests" begin
        Q = matrix(QQ, 4,4,[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]);
        Q1 = -Q;
        L = matrix(QQ, 4, 1 ,[1,1,1,1]);
        c = 3;
        @test Hecke.closest_vectors(Q, L, c)[1] == [-2, -1, -1, -1]
        @test size(Hecke.closest_vectors(Q, L, c), 1) == 9
        @test Hecke.closest_vectors(Q, L, c)[1] == Hecke.closest_vectors(Q1,L,c)[1]

        x = matrix(QQ, 4, 1, Hecke.closest_vectors(Q, L, c)[5]);
        xt = transpose(matrix(QQ, 4, 1, Hecke.closest_vectors(Q, L, c)[5]));
        R = xt * Q * x + 2 * xt * L + c;
        @test  R[1] <= 0

        y = matrix(QQ, 4, 1, Hecke.closest_vectors(Q1, L, c)[5]);
        yt = transpose(matrix(QQ, 4, 1, Hecke.closest_vectors(Q1, L, c)[5]));
        S = yt * Q * y + 2 * yt * L + c;
        @test  S[1] <= 0
    end
end


text:
@testset "scc graph wikipedia" begin
    n = 8
    scc_init = CS.SCCInit(
        zeros(Int, n + 1),
        zeros(Int, n),
        zeros(Int, n),
        zeros(Bool, n),
        zeros(Int, n),
    )
    di_ei = [1, 2, 2, 2, 3, 3, 4, 4, 5, 5, 6, 7, 8, 8]
    di_ej = [2, 2, 3, 5, 7, 4, 3, 8, 6, 1, 7, 6, 4, 7]
    scc_map = CS.scc(di_ei, di_ej, scc_init)
    @test scc_map[1] == scc_map[2] == scc_map[5]
    @test scc_map[3] == scc_map[4] == scc_map[8]
    @test scc_map[6] == scc_map[7]
    @test scc_map[1] != scc_map[3]
    @test scc_map[1] != scc_map[6]
    @test scc_map[3] != scc_map[6]
end

ext: jl


text:
push!(LOAD_PATH, ".")
using ProgressBars
using NPZ
using neural_network

@doc """
Train the neural network on the entire training set.
"""
function run_training_set()
    @info "start reading training data"
    trainingimages::Matrix{Float64} = npzread("dataset/training_images.npy")
    traininglabels::Vector{Int} = parse.(Int, readlines(open("dataset/training_labels.txt")))
    @info "stop reading"

    network::NeuralNetwork = NeuralNetwork([784, 200, 75, 10])

    for (i, image, label) ∈ zip(ProgressBar(1:60000), eachcol(trainingimages), traininglabels)
        train!(network, image, label)
    end

    save_params(network)
end

@time run_training_set()

max_stars_repo_name: ShakyaMajumdar/neural-network

max_issues_repo_name: ShakyaMajumdar/neural-network

max_forks_repo_name: ShakyaMajumdar/neural-network


text:
function hempty{T<:Number}(A::HMat{T})
    C = HMat{T}();
    C.height = A.height;
    C.width = A.width;
    C.trg = A.trg;
    C.src = A.src;
    C.level = A.level;

    C.EPS = A.EPS;
    C.MAXRANK = A.MAXRANK;
    C.MINN = A.MINN;

    if A.blockType == LOWRANK
        C.blockType = LOWRANK;
        C.UMat = zeros(C.height,0);
        C.VMat = zeros(C.width,0);
    elseif A.blockType == DENSE
        C.blockType = DENSE;
        C.DMat = zeros(C.height,C.width);
    elseif A.blockType == HMAT
        C.blockType = HMAT;
        C.childHMat = Array(HMat{T},4,4);
        for i = 1:4
            for j = 1:4
                C.childHMat[i,j] = hempty(A.childHMat[i,j]);
            end
        end
    end
    return C;
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
import Base: broadcast, broadcast!

for (op1,op2) in ((:.+,:+), (:.-,:-), (:.*,:*))
    @eval begin
        @generated function broadcast!{T,N}(::typeof($op1), y::CuArray{T,N}, x1::CuArray{T,N}, x2::CuArray{T,N})
            op = $op2
            f = CuFunction("""
            __global__ void f(Array<$T,$N> y, Array<$T,$N> x1, Array<$T,$N> x2) {
                int idx = blockIdx.x * blockDim.x + threadIdx.x;
                if (idx < y.length()) {
                    int subs[$N];
                    y.idx2sub(idx, subs);
                    y[idx] = x1(subs) $op x2(subs);
                }
            }""")
            quote
                $f(y, x1, x2, dx=length(y))
                y
            end
        end
    end
end


text:
@testset "NO_FIELDS" begin
    @test (@test_deprecated NO_FIELDS) isa NoTangent
end

@testset "extern" begin
    @test extern(@thunk(3)) == 3
    @test extern(@thunk(@thunk(3))) == 3

    @test extern(Tangent{Foo}(x=2.0)) == (;x=2.0)
    @test extern(Tangent{Tuple{Float64,}}(2.0)) == (2.0,)
    @test extern(Tangent{Dict}(Dict(4 => 3))) == Dict(4 => 3)

    # with differentials on the inside
    @test extern(Tangent{Foo}(x=@thunk(0+2.0))) == (;x=2.0)
    @test extern(Tangent{Tuple{Float64,}}(@thunk(0+2.0))) == (2.0,)
    @test extern(Tangent{Dict}(Dict(4 => @thunk(3)))) == Dict(4 => 3)

    z = ZeroTangent()
    @test extern(z) === false
    dne = NoTangent()
    @test_throws Exception extern(dne)
    E = ChainRulesCore.NotImplementedException
    @test_throws E extern(ni)
end


text:
using Hawk
using Base.Test

@constants begin
    a=10
end

@test isconst(:a)


text:
import Calculus:differentiate,jacobian,simplify,SymbolParameter,isminus

function simplify(::SymbolParameter{:-}, args)
    # Remove redundant subtractions: (x-0) == x
    if length(args) == 2 && args[2] == 0
        return args[1]
    end
    # Remove any 0's in a subtraction
    args = map(simplify, filter(x -> x != 0, args))
    if length(args) == 0
        return 0
    # Special Case: simplify(:(x - x)) == 0
    elseif length(args) == 2 && args[1] == args[2]
        return 0
    # Special Case: simplify(:(x - (-y))) == x + y
    elseif length(args) == 2 && isminus(args[2])
        return Expr(:call, :+, args[1], args[2].args[2])
    else
        return Expr(:call, :-, args...)
    end
end

function differentiate(ex::Expr, wrt::Expr)
    if ex.head == :vect || ex.head == :vcat
        return differentiate(SymbolParameter(:vect), ex.args[1:end], wrt)
    elseif ex.head == :ref
        return ex == wrt ? 1 : 0
    elseif ex.head != :call
        error("Unrecognized expression $ex")
    end
    return simplify(differentiate(SymbolParameter(ex.args[1]), ex.args[2:end], wrt))
end

function differentiate(ex::Expr, targets::Vector{Expr})
    n = length(targets)
    exprs = Array(Any, n)
    for i in 1:n
        exprs[i] = differentiate(ex, targets[i])
    end
    return exprs
end

differentiate(ex::Union{Number, Symbol}, wrt::Expr) = 0

function differentiate(::SymbolParameter{:vect}, args, wrt)
    for i = 1:length(args)
        args[i] = differentiate(args[i], wrt)
    end
    return Expr(:vect, args...)
end

function jacobian(ex::Expr, targets::Vector)
    @assert ex.head == :vect || ex.head == :vcat
    exprs = Expr(:vcat)
    for i = 1:length(ex.args)
        push!(exprs.args, Expr(:row))
        for j = 1:length(targets)
            push!(exprs.args[i].args, simplify(differentiate(ex.args[i], targets[j])))
        end
    end
    return exprs
end

differentiate(::SymbolParameter{:max}, args, wrt) = Expr(:if, :($(args[1]) > $(args[2])), differentiate(args[1], wrt), differentiate(args[2], wrt))

differentiate(::SymbolParameter{:min}, args, wrt) = Expr(:if, :($(args[1]) < $(args[2])), differentiate(args[1], wrt), differentiate(args[2], wrt))

differentiate(::SymbolParameter{:abs}, args, wrt) = Expr(:call, :ifelse, :($(args[1]) < 0), -differentiate(args[1], wrt), differentiate(args[1], wrt))

differentiate(::SymbolParameter{:>}, args, wrt) = 0
differentiate(::SymbolParameter{:<}, args, wrt) = 0
differentiate(::SymbolParameter{:≤}, args, wrt) = 0
differentiate(::SymbolParameter{:≥}, args, wrt) = 0

function simplify(ex::Expr)
    if in(ex.head, [:vcat,:row])
        for i = 1:length(ex.args)
            ex.args[i] = simplify(ex.args[i])
        end
        return ex
    end

    if ex.head != :call
        return ex
    end

    if all(a -> isa(a, Number), ex.args[2:end]) && length(ex.args) > 1
        return eval(current_module(), ex)
    end
    new_ex = simplify(SymbolParameter(ex.args[1]), ex.args[2:end])
    while new_ex != ex
        new_ex, ex = simplify(new_ex), new_ex
    end
    return new_ex
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
tp = [Point(0.1234, 0.9876),
      Point(0.0000, 2.0000),
      Point(0.1234, 0.9876),
      Point(0.1234, 0.9876),
      Point(0.1234, 0.9876)]

tq = [Point(0.8765, 0.2345),
      Point(0.0000, 0.0000),
      Point(0.1234, 0.9876),
      Point(0.8765, 0.2345),
      Point(0.1234, 0.9876)]

tr = [2.0, 1.0, 2.0, 0.5, 0.0]

println("Testing circlepoints:")
for i in 1:length(tp)
    p = tp[i]
    q = tq[i]
    r = tr[i]
    (cp, rstatus) = circlepoints(p, q, r)
    println(@sprintf("(%.4f, %.4f), (%.4f, %.4f), %.4f => %s",
                     p.x, p.y, q.x, q.y, r, rstatus))
    for c in cp
        println(@sprintf("    (%.4f, %.4f), %.4f",
                         c.c.x, c.c.y, c.r))
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
Copied from https://github.com/JuliaImages/Images.jl/blob/master/src/algorithms.jl
to ensure that this package doesn't need to depend on Images.jl.
In future the intergral image code should be its own package.
Work in this direction already started here:
https://github.com/JuliaImages/IntegralArrays.jl
=#

"""
```
integral_img = integral_image(img)
```
Returns the integral image of an image. The integral image is calculated by assigning
to each pixel the sum of all pixels above it and to its left, i.e. the rectangle from
(1, 1) to the pixel. An integral image is a data structure which helps in efficient
calculation of sum of pixels in a rectangular subset of an image. See `boxdiff` for more
information.
"""
function integral_image(img::AbstractArray)
    integral_img = Array{accum(eltype(img))}(undef, size(img))
    sd = coords_spatial(img)
    cumsum!(integral_img, img, dims=sd[1])
    for i = 2:length(sd)
        cumsum!(integral_img, integral_img, dims=sd[i])
    end
    integral_img
end

"""
```
sum = boxdiff(integral_image, ytop:ybot, xtop:xbot)
sum = boxdiff(integral_image, CartesianIndex(tl_y, tl_x), CartesianIndex(br_y, br_x))
sum = boxdiff(integral_image, tl_y, tl_x, br_y, br_x)
```
An integral image is a data structure which helps in efficient calculation of sum of pixels in
a rectangular subset of an image. It stores at each pixel the sum of all pixels above it and to
its left. The sum of a window in an image can be directly calculated using four array
references of the integral image, irrespective of the size of the window, given the `yrange` and
`xrange` of the window. Given an integral image -
        A - - - - - - B -
        - * * * * * * * -
        - * * * * * * * -
        - * * * * * * * -
        - * * * * * * * -
        - * * * * * * * -
        C * * * * * * D -
        - - - - - - - - -
The sum of pixels in the area denoted by * is given by S = D + A - B - C.
"""
boxdiff(int_img::AbstractArray{T, 2}, y::UnitRange, x::UnitRange) where {T} = boxdiff(int_img, y.start, x.start, y.stop, x.stop)
boxdiff(int_img::AbstractArray{T, 2}, tl::CartesianIndex, br::CartesianIndex) where {T} = boxdiff(int_img, tl[1], tl[2], br[1], br[2])

function boxdiff(int_img::AbstractArray{T, 2}, tl_y::Integer, tl_x::Integer, br_y::Integer, br_x::Integer) where T
    sum = int_img[br_y, br_x]
    sum -= tl_x > 1 ? int_img[br_y, tl_x - 1] : zero(T)
    sum -= tl_y > 1 ? int_img[tl_y - 1, br_x] : zero(T)
    sum += tl_y > 1 && tl_x > 1 ? int_img[tl_y - 1, tl_x - 1] : zero(T)
    sum
end

accum(::Type{T}) where {T<:Integer} = Int
accum(::Type{Float32})    = Float32
accum(::Type{T}) where {T<:Real} = Float64
accum(::Type{C}) where {C<:Colorant} = base_colorant_type(C){accum(eltype(C))}


text:
g10 = CompleteGraph(10)
@test local_clustering_coefficient(g10) == ones(10)
@test global_clustering_coefficient(g10) == 1
@test local_clustering(g10) == (fill(36, 10), fill(36, 10))
@test triangles(g10) == fill(36, 10)
@test triangles(g10, 1) == 36


text:
"""
 Nufft object for 2d interpolation.
 It contains fft plan and 1d array to pass data to
 nufft2d subroutine from nufft package
"""
struct :: NUFFT2D

    f1d      :: Vector{Complex{Float64}}
    f1       :: Vector{Complex{Float64}}
    f2       :: Vector{Complex{Float64}}
    fcmplx   :: Array{Complex{Float64},2}
    i        :: Vector{Int64}
    j        :: Vector{Int64}
    x        :: Vector{Float64}
    y        :: Vector{Float64}
    epsnufft :: Float64
    n1       :: Int64
    n2       :: Int64
    x1min    :: Float64
    x1max    :: Float64
    x2min    :: Float64
    x2max    :: Float64

    function NUFFT2D(n1, x1min, x1max, n2, x2min, x2max)

        x      = zeros(Float64, n1*n2)
        y      = zeros(Float64, n1*n2)
        i      = zeros(Int64,   n1*n2)
        j      = zeros(Int64,   n1*n2)
	fcmplx = zeros(Complex{Float64}, (n1,n2))
        f1     = zeros(Complex{Float64}, (n1))
        f2     = zeros(Complex{Float64}, (n2))
        f1d    = zeros(Complex{Float64}, (n1*n2))

        new( f1d, f1, f2, fcmplx, i, j, x, y, epsnufft,
	n1, n2, x1min, x1max, x2min, x2max)
    end

end

"""
Compute the fft and prepare data for nufft call
Interpolate value when the fft is already computed.
"""
function interpolate!( nufft, f )

    n1 = nufft.n1
    n2 = nufft.n2

    fft!(fcmplx,f)

    m = n2÷2
    for i=1:n1
      f2=view(fcmplx,i,:)
      for j =
      fcmplx[i,:]=f2([m+1:n2;1:m])
    end
    m = n1/2
    for j=1,n2
        f1=view(fcmplx,:,j)
        fcmplx[:,j]=f1([[(p,p=m+1,n1)], [(p,p=1,m)]])
    end

    xij = (x - x1min) / (x1max - x1min)
    yij = (y - x2min) / (x2max - x2min)

    if ( 0 < xij && xij < 1 && 0 < yij && yij < 1 )

      nufft.i[1] = 1
      nufft.j[1] = 1
      nufft.x[1] = xij * 6π
      nufft.y[1] = yij * 6π

      nufft2d2!(x, y, f1d, 1, nufft.eps, nufft.fcmplx)

    end

end

using FINUFFT
using Random

Random.seed!(1)

nj = 10
nk = 11
ms = 12
mt = 13
mu = 14

tol = 1e-15

# nonuniform data
x = 3π*(1.0 .- 2*rand(nj))
y = 3π*(1.0 .- 2*rand(nj))
c = rand(nj) + 1im*rand(nj)
s = rand(nk)
t = rand(nk)
u = rand(nk)
f = rand(nk) + 1im*rand(nk)

# uniform data
F2D = rand(ms, mt) + 1im*rand(ms,mt)

modevec(m) = -floor(m/2):floor((m-1)/2+1)
k1 = modevec(ms)
k2 = modevec(mt)

#int finufft2d2(BIGINT nj,FLT* xj,FLT *yj,CPX* cj,int iflag,FLT eps,
#             BIGINT ms, BIGINT mt, CPX* fk, nufft_opts opts)
#
# Type-2 2D complex nonuniform FFT.
#
#   cj[j] =  SUM   fk[k1,k2] exp(+/-i (k1 xj[j] + k2 yj[j]))      for j = 0,...,nj-1
#           k1,k2
#   where sum is over -ms/2 <= k1 <= (ms-1)/2, -mt/2 <= k2 <= (mt-1)/2,
#
# Inputs:
#   nj     number of targets (int64, aka BIGINT)
#   xj,yj     x,y locations of targets (each a size-nj FLT array) in [-3pi,3pi]
#   fk     FLT complex array of Fourier transform values (size ms*mt,
#          increasing fast in ms then slow in mt, ie Fortran ordering).
#          Along each dimension the ordering is set by opts.modeord.
#   iflag  if >=0, uses + sign in exponential, otherwise - sign (int)
#   eps    precision requested (>1e-16)
#   ms,mt  numbers of Fourier modes given in x and y (int64)
#          each may be even or odd;
#          in either case the mode range is integers lying in [-m/2, (m-1)/2].
#   opts   struct controlling options (see finufft.h)
# Outputs:
#   cj     size-nj complex FLT array of target values
#          (ie, stored as 2*nj FLTs interleaving Re, Im).
#   returned value - 0 if success, else see ../docs/usage.rst

@testset "FINUFFT 2D2" begin
    out = complex(zeros(nj))
    ref = complex(zeros(nj))
    for j=1:nj
        for ss=1:ms
            for tt=1:mt
                ref[j] += F2D[ss, tt] * exp(1im*(k1[ss]*x[j]+k2[tt]*y[j]))
            end
        end
    end
    nufft2d2!(x, y, out, 1, tol, F2D)
    relerr_2d2 = norm(vec(out)-vec(ref), Inf) / norm(vec(ref), Inf)
    @test relerr_2d2 < 1e-13
    out2 = nufft2d2(x, y, 1, tol, F2D)
    reldiff = norm(vec(out)-vec(out2), Inf) / norm(vec(out), Inf)
    @test reldiff < 1e-14
end


text:
function sector_performance(; client = GLOBAL[], datatype::String = "json", parser = "default")
    @argcheck in(datatype, ["json", "csv"])
    uri = _form_uri_head(client, "SECTOR") * _form_uri_tail(client, nothing, nothing)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end


text:
using MetidaNCA
using Documenter, Weave, PrettyTables, CSV, DataFrames
#using DocumenterLaTeX

include("validation.jl")

#v_out_path = joinpath(dirname(@__FILE__), "src", "validation_report.md")

makedocs(
        modules = [MetidaNCA],
        sitename = "MetidaNCA.jl",
        authors = "Vladimir Arnautov",
        pages = [
            "Home" => "index.md",
            "Examples" => "examples.md",
            "Details" => "details.md",
            "Parameter list" => "parameters.md",
            "API" => "api.md"
            ],
        )

deploydocs(repo = "github.com/PharmCat/MetidaNCA.jl.git", devbranch = "main", forcepush = true
)


text:
# using Revise

using IncrementalInference
using Test

##

@testset "Test priors" begin

##

N=100
graphinits = [false, true]

# TEMP value when not using for loop
graphinit = true

for  graphinit = graphinits

fg = initfg()
fg.solverParams.N = N
fg.solverParams.graphinit = graphinit
fg.solverParams.treeinit = !graphinit

addVariable!(fg, :x0, ContinuousScalar, N=N)
addFactor!(fg, [:x0], Prior(Normal(-1.0, 1.0)))

addVariable!(fg, :x1, ContinuousScalar, N=N)

addVariable!(fg, :x2, ContinuousScalar, N=N)
addFactor!(fg, [:x2], Prior(Normal(+1.0, 1.0)))

addFactor!(fg, [:x0; :x1], LinearRelative(Normal(0.0, 0.01)))
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(0.0, 0.01)))

#solve
tree = solveTree!(fg)
x0_m = getKDEMean(getKDE(getVariable(fg, :x0)))[1]
x1_m = getKDEMean(getKDE(getVariable(fg, :x1)))[1]
x2_m = getKDEMean(getKDE(getVariable(fg, :x2)))[1]

@info ("Testing means = 0 with 2 priors:\ngraphinit=$graphinit\nMeans: x0: $(x0_m), x1: $x1_m, x2: $x2_m")

@test_skip isapprox(x0_m, 0.0, atol = 0.1)
@test_skip isapprox(x1_m, 0.0, atol = 0.1)
@test_skip isapprox(x2_m, 0.0, atol = 0.1)

@warn "priorusetest.jl is testing with large tolerances"
@test isapprox(x0_m, 0.0, atol = 1.0)
@test isapprox(x1_m, 0.0, atol = 1.0)
@test isapprox(x2_m, 0.0, atol = 1.0)

#testing if values are close to one another
testvals = [x0_m, x1_m, x2_m]
meanval = mean(testvals)
@test all(isapprox.(testvals, meanval, atol=0.4))

end

##

for  graphinit = graphinits

fg = initfg()
fg.solverParams.N = N
fg.solverParams.graphinit = graphinit
fg.solverParams.treeinit = !graphinit

addVariable!(fg, :x0, ContinuousScalar, N=N)
addFactor!(fg, [:x0], Prior(Normal(-1.0, 1.0)))

addVariable!(fg, :l0, ContinuousScalar, N=N)
addFactor!(fg, [:l0], Prior(Normal(+1.0, 1.0)))

addVariable!(fg, :l1, ContinuousScalar, N=N)

addFactor!(fg, [:x0; :l0], LinearRelative(Normal(0, 0.01)))
addFactor!(fg, [:x0; :l1], LinearRelative(Normal(0, 0.01)))

addVariable!(fg, :x1, ContinuousScalar, N=N)
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(0, 0.01)))

addVariable!(fg, :x2, ContinuousScalar, N=N)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(0, 0.01)))

addFactor!(fg, [:x2; :l0], LinearRelative(Normal(0, 0.01)))
addFactor!(fg, [:x2; :l1], LinearRelative(Normal(0, 0.01)))

#solve
tree = solveTree!(fg)

x0_m = getKDEMean(getKDE(getVariable(fg, :x0)))[1]
x1_m = getKDEMean(getKDE(getVariable(fg, :x1)))[1]
x2_m = getKDEMean(getKDE(getVariable(fg, :x2)))[1]
l0_m = getKDEMean(getKDE(getVariable(fg, :l0)))[1]
l1_m = getKDEMean(getKDE(getVariable(fg, :l1)))[1]

@info ("Testing means = 0 with 2 priors:\ngraphinit=$graphinit\nMeans: x0: $(x0_m), x1: $x1_m, x2: $x2_m, l0: $l0_m, l1: $l1_m")

@test_skip isapprox(x0_m, 0.0, atol = 0.1)
@test_skip isapprox(x1_m, 0.0, atol = 0.1)
@test_skip isapprox(x2_m, 0.0, atol = 0.1)
@test_skip isapprox(l0_m, 0.0, atol = 0.1)
@test_skip isapprox(l1_m, 0.0, atol = 0.1)

@warn "priorusetest.jl is testing with large tolerances"
@test isapprox(x0_m, 0.0, atol = 1.0)
@test isapprox(x1_m, 0.0, atol = 1.0)
@test isapprox(x2_m, 0.0, atol = 1.0)
@test isapprox(l0_m, 0.0, atol = 1.2)
@test isapprox(l1_m, 0.0, atol = 1.2)

#testing if values are close to one another
@show testvals = [x0_m, x1_m, x2_m, l0_m, l1_m]
@show meanval = mean(testvals)
@test all(isapprox.(testvals, meanval, atol=0.3))

end

##

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the Julia API to convert between the representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/api.jl
# ==============================

# Conversion to DCM
# ------------------------------------------------------------------------------

@testset "Julia conversion API: To DCM" begin
    for T in (Float32, Float64)
        T = Float64

        # Euler angles
        # ======================================================================

        ea = rand(EulerAngles{T})
        D_exp = angle_to_dcm(ea)
        D_api = convert(DCM, ea)
        @test D_api === D_exp
        @test eltype(D_api) === T

        # Euler angle and axis
        # ======================================================================

        av = rand(EulerAngleAxis{T})
        D_exp = angleaxis_to_dcm(av)
        D_api = convert(DCM, av)
        @test D_api === D_exp
        @test eltype(D_api) === T

        # Quaternion
        # ======================================================================

        q = rand(Quaternion{T})
        D_exp = quat_to_dcm(q)
        D_api = convert(DCM, q)
        @test D_api === D_exp
        @test eltype(D_api) === T
    end
end

# Conversion to Euler angle and axis
# ------------------------------------------------------------------------------

@testset "Julia conversion API: To Euler angle and axis" begin
    for T in (Float32, Float64)
        # DCMs
        # ======================================================================

        dcm = rand(DCM{T})
        av_exp = dcm_to_angleaxis(dcm)
        av_api = convert(EulerAngleAxis, dcm)
        @test av_exp === av_api
        @test eltype(av_api) === T

        # Euler angles
        # ======================================================================

        ea = rand(EulerAngles{T})
        av_exp = angle_to_angleaxis(ea)
        av_api = convert(EulerAngleAxis, ea)
        @test av_exp === av_api
        @test eltype(av_api) === T

        # Quaternion
        # ======================================================================

        q = rand(Quaternion{T})
        av_exp = quat_to_angleaxis(q)
        av_api = convert(EulerAngleAxis, q)
        @test av_exp === av_api
        @test eltype(av_api) === T
    end
end

# Conversion to Euler angles
# ------------------------------------------------------------------------------

@testset "Julia conversion API: To Euler angles" begin
    for rot_seq in valid_rot_seqs
        for T in (Float32, Float64)
            # DCM
            # ==================================================================

            dcm = rand(DCM{T})
            ea_exp = dcm_to_angle(dcm, rot_seq)
            ea_api = convert(EulerAngles(rot_seq), dcm)
            @test ea_exp === ea_api
            @test eltype(ea_api) === T

            # Euler angle
            # ==================================================================

            ea = rand(EulerAngles{T})
            ea_exp = angle_to_angle(ea, rot_seq)
            ea_api = convert(EulerAngles(rot_seq), ea)
            @test ea_api === ea_exp
            @test eltype(ea_api) === T

            # Euler angle and axis
            # ==================================================================

            av = rand(EulerAngleAxis{T})
            ea_exp = angleaxis_to_angle(av, rot_seq)
            ea_api = convert(EulerAngles(rot_seq), av)
            @test ea_api === ea_exp
            @test eltype(ea_api) === T

            # Quaternion
            # ==================================================================

            q = rand(Quaternion{T})
            ea_exp = quat_to_angle(q, rot_seq)
            ea_api = convert(EulerAngles(rot_seq), q)
            @test ea_exp === ea_api
            @test eltype(ea_api) === T
        end
    end

    # Conversion without specifying the rotation sequence.
    for T in (Float32, Float64)
        # DCM
        # ======================================================================

        dcm = rand(DCM{T})
        ea_exp = dcm_to_angle(dcm, :ZYX)
        ea_api = convert(EulerAngles, dcm)
        @test ea_exp === ea_api
        @test eltype(ea_api) === T

        # Euler angle
        # ======================================================================

        ea = rand(EulerAngles{T})
        ea_api = convert(EulerAngles, ea)
        @test ea_api === ea
        @test eltype(ea_api) === T

        # In a previous version, there was a bug in which this code was changing
        # the rotation sequence to `:ZYX`.
        vea = EulerAngles[ea]
        @test vea[1] === ea

        # Euler angle and axis
        # ======================================================================

        # Sample a random Euler angle and axis.
        av = rand(EulerAngleAxis{T})
        ea_exp = angleaxis_to_angle(av, :ZYX)
        ea_api = convert(EulerAngles, av)
        @test ea_api === ea_exp
        @test eltype(ea_api) === T

        # Quaternion
        # ======================================================================

        q = rand(Quaternion{T})
        ea_exp = quat_to_angle(q)
        ea_api = convert(EulerAngles, q)
        @test ea_exp === ea_api
        @test eltype(ea_api) === T
    end
end

# Conversion to quaternion
# ------------------------------------------------------------------------------

@testset "Julia conversion API: To Quaternion" begin
    for T in (Float32, Float64)
        # Euler angles
        # ======================================================================

        ea = rand(EulerAngles{T})
        q_exp = angle_to_quat(ea)
        q_api = convert(Quaternion, ea)
        @test q_api === q_exp
        @test eltype(q_api) === T

        # Euler angle and axis
        # ======================================================================

        # Sample a random Euler angle and axis.
        av = rand(EulerAngleAxis{T})
        q_exp = angleaxis_to_quat(av)
        q_api = convert(Quaternion, av)
        @test q_api === q_exp
        @test eltype(q_api) === T

        # DCM
        # ======================================================================

        D = rand(DCM{T})
        q_exp = dcm_to_quat(D)
        q_api = convert(Quaternion, D)
        @test q_api === q_exp
        @test eltype(q_api) === T
    end
end


text:
# Period doubling bifurcation + cobweb for logistic map
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

lo = Systems.logistic(0.4)
rs = (2.8, 3.3)

fig, axs = subplots(1, 2, figsize = (0.66figx, figy))

T = 1000
xs = 0:0.0001:1

function cobweb(t) # transform timeseries t into cobweb (point2D)
    cx = Float64[]; cy = Float64[]
    for i ∈ 1:length(t)-1
        push!(cx, t[i]); push!(cy, t[i])
		push!(cx, t[i]); push!(cy, t[i+1])
    end
    # add line to x axis
    pushfirst!(cx, cx[1])
    pushfirst!(cy, 0)
    return cx, cy
end

for (i, r) in enumerate(rs)
    ax = axs[i]
    ax.set_xlim(0,1)
    ax.set_ylim(0,1)
    ax.set_xticks([0, 1])
    ax.set_xlabel("\$x\$"; labelpad=-30)
	set_parameter!(lo, 1, r)
	f  = lo.f.(xs, Ref([r]), 0)
    f² = lo.f.(f, Ref([r]), 0)
	ax.plot(xs, f, color = "C0")
	ax.plot([0, 1], [0,1]; lw = 2, color = "C3", ls = "--")
	x = trajectory(lo, T)
    cx, cy = cobweb(x)
    if i>1
        ax.plot(xs, f², color = "C1")
        ax.text(0.15, 0.875, "\$f^{(2)}\$"; color = "C1", transform = ax.transAxes, size = 32)
    end
    ax.plot(cx, cy; lw = 1.0, color = "C2")
    ax.text(0.05, 0.875, "\$f\$"; color = "C0", transform = ax.transAxes, size = 32)
end
axs[2].set_yticklabels([])

fig.tight_layout(pad = 0.37)
add_identifiers!(fig)

wsave(plotsdir("4", "period_doubling_bifurcation"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Colors
using FixedPointNumbers

typealias RGB8bit ColorTypes.RGB{FixedPointNumbers.UFixed{UInt8,8}}

type Gradient
   startcolor::RGB8bit
   endcolor::RGB8bit
   nilcolor::RGB8bit
end

type ImageSpec
   width::Int64
   height::Int64
   num_supersamples::Int64
end

type CameraSpec
   gradient::Gradient
   center::Complex
   zoom::Float64
   pow::Float64
   affine::Bool
end

type FractalSpec
   fracfunc_closure::Function
end

max_stars_repo_name: jblondin/Fractal

max_issues_repo_name: jblondin/Fractal

max_forks_repo_name: jblondin/Fractal


text:
using SparseArrays, LinearAlgebra, LuxurySparse

# TODO: add documents
# 1. generate a table from SYM_LIST
# 2. auto attach docs to each gate instance and type

# note!!!
# Here, we don't need gates like swap^α and √X.
# swap^α = rot(SWAP, -α*π)
# √X = rot(X, -α*π)
# since both of them are reflexive.
const ConstGateDefaultType = ComplexF64
const SYM_LIST = [
    (:P0, sparse(ConstGateDefaultType[1 0;0 0])),
    (:P1, sparse(ConstGateDefaultType[0 0;0 1])),
    (:X, PermMatrix([2,1], ConstGateDefaultType[1+0im, 1])),
    (:Y, PermMatrix([2,1], ConstGateDefaultType[-im, im])),
    (:Z, Diagonal(ConstGateDefaultType[1+0im, -1])),
    (:S, Diagonal(ConstGateDefaultType[1, im])),
    (:Sdag, Diagonal(ConstGateDefaultType[1, -im])),
    (:T, Diagonal(ConstGateDefaultType[1, exp(π*im/4)])),
    (:Tdag, Diagonal(ConstGateDefaultType[1, exp(-π*im/4)])),
    (:I2, IMatrix{2, ConstGateDefaultType}()),
    (:H, (elem = 1 / sqrt(2); ConstGateDefaultType[elem elem; elem -elem])),
    (:CNOT, PermMatrix([1, 4, 3, 2], ones(ConstGateDefaultType, 4))),
    (:SWAP, PermMatrix([1, 3, 2, 4], ones(ConstGateDefaultType, 4))),
    (:Toffoli, PermMatrix([1, 2, 3, 8, 5, 6, 7, 4], ones(ConstGateDefaultType, 8))),
    (:Pu, sparse([1], [2], ConstGateDefaultType[1+0im], 2, 2)),
    (:Pd, sparse([2], [1], ConstGateDefaultType[1+0im], 2, 2)),
]

for (NAME, MAT) in SYM_LIST
    GT = Symbol(NAME, "Gate")
    @eval begin
        @const_gate $NAME = $MAT
    end
end

Base.adjoint(::PuGate) = Pd
Base.adjoint(::PdGate) = Pu
Base.adjoint(::SGate) = Sdag
Base.adjoint(::SdagGate) = S
Base.adjoint(::TGate) = Tdag
Base.adjoint(::TdagGate) = T

# Docs
"""
    X
    XGate <: ConstantGate{1}

Pauli X gate. `X` is the instance of `XGate`.
"""
X, XGate

"""
    Y
    YGate  <: ConstantGate{1}

Pauli Y gate. `Y` is the instance of `YGate`.
"""
Y, YGate

"""
    Z
    ZGate  <: ConstantGate{1}

Pauli Z gate. `Z` is the instance of `YGate`.
"""
Z, ZGate


text:
module FrameFun

using FillArrays, GridArrays, StaticArrays
using LinearAlgebra, LowRankApprox
using Printf, RecipesBase, Reexport
using DomainSets
using QuadGK
using CompositeTypes, CompositeTypes.Display

import CompositeTypes: components, component, ncomponents
import DomainSets: factors, nfactors, factor

@reexport using BasisFunctions
import BasisFunctions: coefficienttype, dual

include("extra/FrameFunDomains/FrameFunDomains.jl")
@reexport using .FrameFunDomains

include("platforms/generic/platforms.jl")
include("platforms/generic/parameterpaths.jl")

include("frames/extensionframes.jl")

include("generic/approximationproblems.jl")
include("generic/interface/interface.jl")

# Other platforms
include("platforms/augmentation.jl")
include("platforms/extensionplatform.jl")
include("platforms/weightedsumplatform.jl")
include("platforms/sumplatform.jl")
include("platforms/basesandframes.jl")

include("generic/interface/interface_extension.jl")

export residual, abserror, maxerror, L2error
include("generic/error.jl")

include("applications/DiffEquations.jl")
@reexport using .DiffEquations
# include("applications/WeightedApproximation.jl")
# @reexport using .WeightedApproximation
include("applications/high_dimensional.jl")

include("generic/adaptivity.jl")

include("recipes.jl")

end # module


text:
module preprocessing

export PreprocessingOption, generatePreprocessingOption, preprocessFunction, preprocessPackage!, stdFunctionPreprocess!, preprocessStacktrace, retrieveDocfreq!, DocFreqRemoval, removeCommonFunction!, createDocFreqRemoval

using SHA
using DataStructures
using util

struct Frame
    func::String
    is_crash::Bool # Only for Gnome True if this is the frame where we crashed. For example, in a GDB trace, is_crash true if this is the frame where the signal handler was called.
end

Frame(func::String) = Frame(func, false)

function createFrame(frameJson)
    # Consider frame without
    func_name = get(frameJson, "function", nothing)
    func_name =  isnothing(func_name) ? ""  :  string(func_name)

    is_crash = get(frameJson, "is_crash", "0") == "1"
    if !is_crash && occursin("<signal handler called>", func_name)
        # println("$(func_name) should be is_crash")
        is_crash = true
    end

    return Frame(func_name, is_crash)
end

struct PreprocessingOption
    max_depth::Int32
    filter_recursion::String
    ukn_tkn::String
    unique_ukn_report::Bool
    trim_len::Int32
    rm_sub_stacks::Bool
    select_interesting_stacks::Bool
end

function generatePreprocessingOption(;max_depth=300, filter_recursion="none", ukn_tkn="??", unique_ukn_report=true,
    preprocess_func="function", rm_sub_stacks=true, select_interesting_stacks=false)

    trim_len = 0

    if preprocess_func == "function"
        trim_len = 0
    elseif preprocess_func == "class"
        trim_len = 1
    elseif preprocess_func == "package"
        trim_len = 2
    end

    @info "Trim length: $(trim_len)"

    if select_interesting_stacks
        @warn "select_interesting_stacks only works for Gnome!"
    end

    return PreprocessingOption(max_depth, filter_recursion, ukn_tkn, unique_ukn_report, trim_len, rm_sub_stacks, select_interesting_stacks)
end

function preprocessFunction(func_name::String)::String
    # Strip __GI__ in the begining of function call
    func_name = replace(func_name, r"\(.*\)" => "")
    # Strip _ in the begining of function call
    func_name = replace(func_name, r"^_*GI_+" => "")
    # Remove underscore at the beggining of the function
    func_name = replace(func_name, r"^_+" => "")

    # todo: Add S3M preprocessing: remove numbers and punctuations(except , and .) Do not remove functions that only contain numbers
    # todo: transform to lower character as S3M

    # Java anonymous function and inner function
    func_name = replace(func_name, r"[$][0-9]+" => "\$")
    func_name = replace(func_name, r"invoke[0-9]+" => "invoke")

    # func_name = replace(func_name, r"(?<=[a-zA-Z])[0-9]+$" => "" )
    # func_name = replace(func_name, r"[^a-zA-Z.,_<>$]"=> "" )

    return strip(func_name)
end

function trim_function_name(func_name, trim_len, sep='.')
    if trim_len == 0
        return func_name
    end

    last_idx = 0
    ndots = 0

    for idx in Iterators.reverse(eachindex(func_name))
        if func_name[idx] === sep
            ndots+=1

            if ndots == trim_len
                last_idx = idx - 1
                break
            end
        end
    end

    if ndots != trim_len
        return ""
    end

    return func_name[1:last_idx]
end

function preprocess_stack!(vocab::Vocab, frames::MutableLinkedList{Frame}, ukn_tkn, ukn_id, trim_len, max_depth=300)::Vector{UInt32}
    token_id::UInt32 = 0
    vecLen = length(frames) > max_depth ? max_depth : length(frames)
    function_ids = Vector{Int32}()
    sizehint!(function_ids, length(vecLen))

    for (depth, frame) in enumerate(frames)
        depth > max_depth && break

        if length(frame.func) == 0 || frame.func == ukn_tkn
            token_id = ukn_id
        else
            func_name = trim_function_name(frame.func, trim_len)

            if isempty(func_name)
                continue
            end

            func_name = preprocessFunction(frame.func)
            token_id = setdefault!(vocab, func_name)
        end

        push!(function_ids, token_id)
    end

    if length(function_ids) == 0
        @warn "Stack trace is empty. Fill it with a fake token."
        token_id = setdefault!(vocab, "@@@$(length(vocab.vocab))@@@")
        push!(function_ids, token_id)
    end

    return function_ids
end

function stdFunctionPreprocess!(vocab::Vocab, frames::MutableLinkedList{Frame}, ukn_tkn, ukn_id, max_depth=300)::Vector{UInt32}
    vecLen = length(frames) > max_depth ? max_depth : length(frames)

    if length(frames) == 0
        @warn "Stack trace is empty. Fill it with a fake token."
        function_ids = Vector{UInt32}(1)
        function_ids[1] = setdefault!(vocab, "@@@$(length(vocab.vocab))@@@")
        return function_ids
    end

    function_ids = Vector{UInt32}(undef, vecLen)
    token_id::UInt32 = 0

    for (depth, frame) in enumerate(frames)
        depth > max_depth && break

        if length(frame.func) == 0 || frame.func == ukn_tkn
            token_id = ukn_id
        else
            func_name = preprocessFunction(frame.func)
            token_id = setdefault!(vocab, func_name)
        end

        function_ids[depth] = token_id
    end

    return function_ids
end

# vocab = Dict{String, UInt32}()
# preprocessPackage!(vocab, [Dict("function"=> "__gi__func()"),Dict("function"=> "org.eclipse.package.class.method"),Dict("function"=> 1)])
# vocab = Dict{String, UInt32}("UKN" => 0)
# stdFunctionPreprocess(vocab, [Dict("function"=> "__GI__func()"),Dict("function"=> "UKN"),Dict("function"=> 1)], "UKN", 0)

"""
    rm_sub_stack!(stacktraces::Vector{Vector{UInt32}})

    Remove stacks that are sub sequence of other ones. We consider if st1 is subsequence st2 when
        the prefix of st2 is equal to st1, i.e., st1[1] ==st2[1], st1[2] ==st2[2] .. st1[n] ==st2[n] where n is st1 length.
"""
# function rm_sub_stack(stacktraces::Vector{Vector{UInt32}})::Vector{Vector{UInt32}}
function rm_sub_stack!(stacktraces::Vector{Vector{UInt32}})::Vector{Vector{UInt32}}
    length(stacktraces) == 1 && return stacktraces
    idx_to_rm = Vector{Int64}()

    sorted_stacktraces = sort!(stacktraces, by=v -> length(v))

    for (i, stack1) in enumerate(sorted_stacktraces)
        stack1_len = length(stack1)

        for j = i + 1:length(sorted_stacktraces)
            stack2 = sorted_stacktraces[j]

            stack1_len > length(stack2) && continue

            # Prefix of st2 is equal to st1
            frame_idx = 1

            while frame_idx != stack1_len + 1
                stack2[frame_idx] != stack1[frame_idx] && break
                frame_idx += 1
            end

            if frame_idx - 1 === stack1_len
                push!(idx_to_rm, i)
                break
            end
        end
    end

    # if length(idx_to_rm) > 0
    #     println(idx_to_rm)
    #     println(stacktraces)
    # end

    length(idx_to_rm) > 0 && deleteat!(sorted_stacktraces, idx_to_rm)

    return sorted_stacktraces
end

function rmDuplicateStacks(stacktraces::Vector{Vector{UInt32}})::Vector{Vector{UInt32}}
    hash2stack = Dict{String,Vector{UInt32}}([(bytes2hex(sha2_512(string(stack))), stack) for (stack_idx,  stack) in enumerate(stacktraces)])

    if length(hash2stack) == length(stacktraces)
        return stacktraces
    end

    return [stack for stack in values(hash2stack)]
end

# println(rmDuplicateStacks([[0x00001,0x00002,0x00003],[0x00002,0x00003,0x00004,0x00005],[0x00003,0x00003,0x00003,0x00003]]))
# println(rmDuplicateStacks([[0x00001,0x00002,0x00003],[0x00001,0x00002,0x00003],[0x00003,0x00003,0x00003,0x00003]]))

function validateStacktraceFrames!(mainVector, stack, report_id, ukn_tkn,  is_nested)
    new_stacktrace = MutableLinkedList{Frame}()
    last_depth = -1

    for frameJson in stack["frames"]
        depth = typeof(frameJson["depth"]) == Int64 ? frameJson["depth"] : parse(Int64, frameJson["depth"])
        pos_diff = depth - last_depth - 1

        frame = createFrame(frameJson)

        if pos_diff > 0
            # Fix stack trace
            for i in 1:pos_diff
                # It is missing frames from last_depth until depth
                # missing_depth = last_depth + i
                push!(new_stacktrace, Frame(ukn_tkn))
        end

            # @warn "report stacktrace $(report_id) is missing frames."
        elseif pos_diff < 0
            # There are multiple stacks in the stack variable
            if length(new_stacktrace) > 0
                push!(mainVector, new_stacktrace)
            end

            last_depth = depth
            new_stacktrace = MutableLinkedList{Frame}()

            for i in 0:depth - 1
                push!(new_stacktrace, Frame(ukn_tkn))
            end
            # logging.getLogger().error("More than one stack in report {}.".format(report_id, stacktrace['frames']))
        end

        last_depth = depth
        append!(new_stacktrace, frame)
    end

    if length(new_stacktrace) == 0 && !is_nested
        @warn "preprocessed_frame is empty for report $(report_id)"
        return
    end

    push!(mainVector, new_stacktrace)
end

function validateStacktrace!(main_frames, stack, report_id, ukn_tkn)
    validateStacktraceFrames!(main_frames, stack, report_id, ukn_tkn, false)

    nested = get(stack, "nested", nothing)

    if !isnothing(nested) && length(stack["nested"]) > 0
        validateStacktraceFrames!(main_frames, stack["nested"][1], report_id, ukn_tkn, true)
    end
end

function removeRecursiveCalls(frames::Vector{UInt32}, ukn_function::UInt32, recursion_removal::String)::Vector{UInt32}
    if recursion_removal == "none"
        return frames
    end

    clean_stack = MutableLinkedList{UInt32}()

    if recursion_removal == "brodie"
        # Quickly Finding Known Software Problems via Automated Symptom Matching - Mark Brodie 2006
        previous_function = frames[1]
        push!(clean_stack, frames[1])

        for (idx, fr) in enumerate(frames)
            idx == 1 && continue

            if fr == ukn_function || fr != previous_function
                push!(clean_stack, fr)
            end

            previous_function = fr
        end
    elseif recursion_removal == "modani"
        # Automatically Identifying Known Software Problems - Natwar Modani 2007
        idx = 1

        while idx < length(frames) + 1
            fr = frames[idx]
            end_idx = -9999999

            push!(clean_stack, fr)

            if fr != ukn_function
                for bk_idx in (length(frames)):-1:idx + 1
                    bk_func = frames[bk_idx]

                    if fr == bk_func
                        end_idx = bk_idx
                        break
                    end

                end
            end

            if end_idx > 0
                idx = end_idx + 1
            else
                idx += 1
            end
        end
    else
        error("Invalid argument value for recursion_removal: $(recursion_removal)")
    end

    if length(clean_stack) != length(frames)
        return collect(clean_stack)
    else
        return frames
    end
end

function extract_interesting_stacks(stacktraceLinkedList::MutableLinkedList{MutableLinkedList{Frame}})::MutableLinkedList{MutableLinkedList{Frame}}
    # Following: https://bazaar.launchpad.net/~bgo-maintainers/bugzilla-traceparser/3.4/view/head:/lib/TraceParser/Trace.pm line 256
    length(stacktraceLinkedList) == 1 && return stacktraceLinkedList

    interesting_stacks = MutableLinkedList{MutableLinkedList{Frame}}()

    for stack in stacktraceLinkedList
        for frame in stack
            if frame.is_crash
                push!(interesting_stacks, stack)
                break
            end
        end
    end

    length(interesting_stacks)  > 0 && return interesting_stacks

    # Search for threads that have a function with
    # "signal" or "segv" in the name.

    # If a trace lacks <signal handler called>, we determine the
    # "interesting thread" by looking for a thread that has
    # functions that match this regex.
    possible_crash_functions = r"signal|segv|sighandler"i
    for stack in stacktraceLinkedList
            for frame in stack
                if occursin(possible_crash_functions, frame.func)
                    push!(interesting_stacks, stack)
                    break
            end
        end
    end

    length(interesting_stacks)  > 0 && return interesting_stacks

    # If we still don't have a thread, return every first thread whose
    # last function isn't some form of wait or one of the ignored
    # functions.

    # Or if that fails, by a thread whose last function *doesn't* match
    # this regex.
    wait_function =  r"wait|sleep|poll"i

    # However, some wait functions are interesting--for example,
    # if we're waiting on a lock, that's interesting during
    # deadlock traces.
    interesting_wait_function = r"lock"i;
    ignore_functions = Set(["kernel_vsyscall", "libc_start_main",  "raise", "abort", "poll", "??"])

    for stack in stacktraceLinkedList
        first_func = stack[1].func
        if ((!occursin(wait_function, first_func) || occursin(interesting_wait_function, first_func)) && first_func ∉ ignore_functions)
            push!(interesting_stacks, stack)
        end
    end

    length(interesting_stacks)  > 0 && return interesting_stacks

    return stacktraceLinkedList
end

function preprocessStacktrace(report_id, stacktracesJson, vocab::Vocab,  opt::PreprocessingOption)
    ukn_tkn = opt.ukn_tkn

    if opt.unique_ukn_report
        # Create token for unknown function name. Give a unique id to the unknown values in stacktrace
        # Following ABRT, we only compare the function names and consider two unknown function(??) as different
        ukn_tkn = "$(ukn_tkn)$(report_id)"

        haskey(vocab.vocab, ukn_tkn) &&  error("Token for unknown function name in the report $(report_id) already exists.")
    end

    ukn_id = setdefault!(vocab, ukn_tkn, true)

    stacktraceLinkedList = MutableLinkedList{MutableLinkedList{Frame}}()

    if stacktracesJson isa Vector
        for stack in stacktracesJson
            validateStacktrace!(stacktraceLinkedList, stack, report_id, ukn_tkn)
        end
    else
        validateStacktrace!(stacktraceLinkedList, stacktracesJson, report_id, ukn_tkn)
    end

    if opt.select_interesting_stacks
        stacktraceLinkedList = extract_interesting_stacks(stacktraceLinkedList)
    end

    finalStacktraces = collect(stacktraceLinkedList)

    preprocessedStacktraces = Vector{Vector{UInt32}}(undef, length(finalStacktraces))

    for (idx, stack) in enumerate(finalStacktraces)
        preprocessedStacktraces[idx] = preprocess_stack!(vocab, stack, ukn_tkn, ukn_id, opt.trim_len, opt.max_depth)
        preprocessedStacktraces[idx] = removeRecursiveCalls(preprocessedStacktraces[idx], ukn_id, opt.filter_recursion)
    end

    if opt.rm_sub_stacks
        preprocessedStacktraces = rm_sub_stack!(preprocessedStacktraces)
        # preprocessedStacktraces = rmDuplicateStacks(preprocessedStacktraces)
    end

    length(preprocessedStacktraces) == 0 && error("$(report_id) contains 0 zero stacktraces".format(report_id))

    return preprocessedStacktraces
end

abstract type CommonFunctionRemoval end

struct DocFreqRemoval <: CommonFunctionRemoval
    k::Float64 # threshold
    trim::Bool
end

function createDocFreqRemoval(name, k)
    if name == "none"
        return DocFreqRemoval(-1.0, false)
    elseif name == "threshold"
        return DocFreqRemoval(k / 100.00, false)
    elseif name == "threshold_trim"
        return DocFreqRemoval(k / 100.00, true)
    end
end

isToRemove(funcId, removalStrategy, docFreq) = ((docFreq.doc_freq[funcId] > removalStrategy.k * docFreq.nDocs) && (!docFreq.static_df_ukn || funcId ∉ docFreq.vocab.ukn_set))

function removeCommonFunction!(report::Report, removalStrategy::DocFreqRemoval, docFreq::DocFreq)
    if removalStrategy.k <= 0.0
        return
    end

    idxs_to_remove_by_stack = Vector{BitVector}(undef, length(report.stacks))
    n_empty_stacks = 0
    empty_vectors = falses(length(report.stacks))

    for (stack_idx, stacktrace) in enumerate(report.stacks)
        idxs_to_remove =  get_indices_to_remove(stacktrace, removalStrategy, docFreq)
        idxs_to_remove_by_stack[stack_idx] = idxs_to_remove

        if sum(idxs_to_remove) == length(stacktrace)
            n_empty_stacks += 1
            empty_vectors[stack_idx] = true
        end
    end

    if n_empty_stacks == length(report.stacks)
    # If all stacks will be empty, so do not remove the common functions
        return
    end

    for stack_idx = 1:length(report.stacks)
        if empty_vectors[stack_idx]
            continue
        end

        deleteat!(report.stacks[stack_idx], idxs_to_remove_by_stack[stack_idx])
    end

    deleteat!(report.stacks, empty_vectors)
end

function get_indices_to_remove(stack::Vector{UInt32}, removalStrategy::DocFreqRemoval, docFreq::DocFreq)
    idxToRemove = falses(length(stack))

    if removalStrategy.trim
        for (idx, funcId) in enumerate(stack)
            (idx > 1 && !idxToRemove[idx - 1]) && break

            idxToRemove[idx] = isToRemove(funcId, removalStrategy, docFreq)
        end

        for (idx, funcId) in Iterators.reverse(enumerate(stack))
            (idx < length(stack) && !idxToRemove[idx + 1]) && break

            idxToRemove[idx] = isToRemove(funcId, removalStrategy, docFreq)
        end
    else
        for (idx, funcId) in enumerate(stack)
           idxToRemove[idx] = isToRemove(funcId, removalStrategy, docFreq)
        end
    end

    return idxToRemove
    # deleteat!(stack, idxToRemove)
end

end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
# Detailed models for performing mapping to Asap3/Asap4

# Constructors for basic components of the array
include("styles.jl")
include("components.jl")

# Array layouts
include("asap2.jl")
include("asap3.jl")
include("asap4.jl")
include("build_arch.jl")

# Specialized Array Layouds
include("asap3_vpr.jl")


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => Hyperrectangle(low=[fill(0.001, 2); zeros(576)],
                            high=[fill(0.0015, 2); zeros(576)]),
    "info" => Dict(
        "name" => "mna1",
        "dim" => 578,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
module BamIO
export refid, position, l_read_name, mq, n_cigar_op, flag, l_seq, next_refid, next_position, tlen, read_name, cigar, seq, qualities, eof, nextread, reference_name
export is_paired, is_proper_paired, is_unmapped, mate_unmapped, is_reverse, mate_reverse, is_read1, is_read2, not_primary, fail_qc, is_duplicate, supplementary_alignment

include("bam.jl")

# package code goes here

end # module


text:
@testset "Dual Type" begin
    model = JuMP.Model()
    x = @variable(model, [1:2])
    @constraint(model, sum(x) == 1)
    h = hrep(model)
    p = polyhedron(h)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, LPHRep{Float64}, Polyhedra.Hull{Float64, Vector{Float64}, Int}}
    h = hrep(zeros(Int, 2, 2), zeros(Int, 2))
    p = polyhedron(h)
    @test p isa Polyhedra.DefaultPolyhedron{Rational{BigInt}, MixedMatHRep{Rational{BigInt}, Matrix{Rational{BigInt}}}, MixedMatVRep{Rational{BigInt}, Matrix{Rational{BigInt}}}}
    h = hrep(spzeros(2, 2), zeros(2))
    p = polyhedron(h)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, SparseMatrixCSC{Float64, Int}}, MixedMatVRep{Float64, Matrix{Float64}}}
    v = vrep(zeros(2, 3))
    p = polyhedron(v)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}
    v = vrep(spzeros(2, 3))
    p = polyhedron(v)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, SparseMatrixCSC{Float64, Int}}}
end

struct BadPoly{T} <: Polyhedra.Polyhedron{T}
end

@testset "Unimplemented methods" begin
    p = BadPoly{Int}()
    @test_throws ErrorException detectvlinearity!(p)
    @test_throws ErrorException detecthlinearity!(p)
    @test_throws ErrorException eliminate(p, [1], FourierMotzkin())
    @test_throws ErrorException eliminate(p, [1], BlockElimination())
    @test_throws ErrorException hrep(p)
    @test_throws ErrorException vrep(p)
    @test_throws ErrorException loadpolyhedron!(p, "a", "b")
    @test_throws ErrorException loadpolyhedron!(p, "a", "ine")
    @test_throws ErrorException loadpolyhedron!(p, "a", "ext")
end

@testset "Polyhedra.DefaultPolyhedron constructor with nothing" begin
    vr = convexhull([-1, 0], [0, -1]) + conichull([1, 1], [-1, 1])
    hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
    p = Polyhedra.DefaultPolyhedron{Int, typeof(hr), typeof(vr)}(hr, nothing, lp_solver)
    @test hrep(p) === hr
    @test !vrepiscomputed(p)
    @test p.solver === lp_solver
    p = Polyhedra.DefaultPolyhedron{Int, typeof(hr), typeof(vr)}(nothing, vr, lp_solver)
    @test !hrepiscomputed(p)
    @test vrep(p) === vr
    @test p.solver === lp_solver
end


text:
using NeuralVerification
const NV = NeuralVerification

"""
Samples a point in the hyperrectangle [los, his] uniformly at random.
"""
function sample(los::Vector{Float64}, his::Vector{Float64})
    @assert length(los) == length(his)
    n = length(los)
    return los .+ rand(n) .* (his - los)
end

"""
Computes result and intermediate values of nnet(input).
"""
function compute_values(nnet, input; collect_f=last)
    x = input
    elem = collect_f((x, x))
    collected = [elem]

    for layer in nnet.layers
        x̂ = NV.affine_map(layer, x)
        x = layer.activation(x̂)
        push!(collected, collect_f((x̂, x)))
    end

    return collected
end

"""
Checks, if given bounds [lbs, ubs] on all neurons in nnet are violated by
randomly sampling n points from the input_set.
"""
function check_bounds(nnet, input_set::Hyperrectangle, lbs, ubs; n=100, ϵ_tol=1e-10, print_success=true)
    counter_examples = []
    for i in 1:n
        x = sample(low(input_set), high(input_set))
        values = compute_values(nnet, x, collect_f=first)
        for (i, (lb, ub)) in enumerate(zip(lbs, ubs))
            # lbs, ubs don't count input layer
            lb_correct = all(values[i+1] .>= lb .- ϵ_tol)
            ub_correct = all(values[i+1] .<= ub .+ ϵ_tol)

            if ~lb_correct || ~ub_correct
                max_lb_error = maximum(abs.(min.(0, values[i+1] - lb)))
                max_ub_error = maximum(max.(0, values[i+1] - ub))
                println("layer ", i, " wrong for ", x, ", largest deviations: ", [max_lb_error, max_ub_error])
                push!(counter_examples, x)
            end
        end
    end

    if length(counter_examples) == 0 && print_success
        println("test passed")
    elseif length(counter_examples) > 0
        println("FAIL !!!")
    end

    return counter_examples
end

"""
Checks bounds randomly by generating input set with specified width, calculating
bounds for it according to the generate_bounds function and then sampling from the
input set.
"""
function check_random_bounds(nnet, width, generate_bounds; n=100, n_rand=100, ϵ_tol=1e-10, print_individual_tests=true)
    @assert generate_bounds != nothing
    dim = size(nnet.layers[1].weights, 2)
    counter_examples = []

    for i in 1:n_rand
        low = randn(dim)
        high = low .+ width
        input_set = Hyperrectangle(low=low, high=high)

        lbs, ubs = generate_bounds(input_set)

        input_counters = check_bounds(nnet, input_set, lbs, ubs, n=n, ϵ_tol=ϵ_tol, print_success=print_individual_tests)
        if length(input_counters) > 0
            push!(counter_examples, input_counters)
        end
    end

    if length(counter_examples) == 0
        println("All tests passed")
    else
        println("FAIL !!!")
    end

    return counter_examples
end

"""
Creates network with nice weights for hand calculation.
(Nice meaning small integers)

args:
layer_dims (Int64[]) - vector of layer dimensions (including the input layer!)

kwargs:
possible_vals (list) - list to draw weights from, defaults to -1.:1.
"""
function create_nice_random_network(layer_dims::Vector{Int64}; possible_vals=nothing)
    possible_vals = isnothing(possible_vals) ? (-1.:1.) : possible_vals

    layers = []
    n_layers = length(layer_dims) - 1 # because we don't create input layer
    d_last = layer_dims[1]
    for (i, d) in enumerate(layer_dims[2:end])
        weights = rand(possible_vals, d, d_last)
        bias = rand(possible_vals, d)

        if i < n_layers
            push!(layers, NV.Layer(weights, bias, NV.ReLU()))
        else
            push!(layers, NV.Layer(weights, bias, NV.Id()))
        end

        d_last = d
    end

    return Network(layers)
end

"""
Creates network with gaussian distributed weights.

layer_dims (Int64[]) - vector of layer dimensions (including the input layer!)
"""
function create_random_network(layer_dims::Vector{Int64})
    layers = []
    n_layers = length(layer_dims) - 1 # because we don't create input layer
    d_last = layer_dims[1]
    for (i, d) in enumerate(layer_dims[2:end])
        weights = randn(d, d_last)
        bias = randn(d)

        if i < n_layers
            push!(layers, NV.Layer(weights, bias, NV.ReLU()))
        else
            push!(layers, NV.Layer(weights, bias, NV.Id()))
        end

        d_last = d
    end

    return Network(layers)
end

function generate_nice_random_input_set(nnet, width; possible_vals=nothing)
    possible_vals = isnothing(possible_vals) ? (-1.:1.) : possible_vals

    dim = size(nnet.layers[1].weights, 2)
    low = rand(possible_vals, dim)
    high = low .+ width
    return Hyperrectangle(low=low, high=high)
end

function generate_counterexample_nn(width, generate_bounds; layer_range=nothing, layer_width_range=nothing,
                                    n_nns=100, n_sample=100, n_in=2, n_out=1)

    layer_range = isnothing(layer_range) ? (1:10) : layer_range
    layer_width_range = isnothing(layer_width_range) ? (2:2) : layer_width_range

    counter_nets = []
    counters = []
    input_sets = []

    for hl in layer_range
        # number of hidden layers with 2 neurons +1 for input layer +1 for output layer
        layers = rand(layer_width_range, hl+2)
        layers[1] = n_in
        layers[end] = n_out

        for n_nn in 1:n_nns
            net = create_nice_random_network(layers)
            input_set = generate_nice_random_input_set(net, 0.1)
            lbs, ubs = generate_bounds(net, input_set)

            counters_current = check_bounds(net, input_set, lbs, ubs, n=n_sample, print_success=false)

            if length(counters_current) > 0
                println("Found counterexample with ", length(layers)-1, " layers!")
                push!(counter_nets, net)
                push!(counters, counters_current)
                push!(input_sets, input_set)
            end
        end
    end

    return counter_nets, counters, input_sets
end


text:
include("../../imports.jl")  # Load on main node.

# NOTE: Change this.
# SIMNAME = "$(basename(@__DIR__))-test"
SIMNAME = "$(basename(@__DIR__))"
SCRATCH_DIR = "/scratchdata/$(ENV["USER"])/cytof/results/repfam"
RESULTS_DIR_PREFIX = joinpath(SCRATCH_DIR, SIMNAME)
AWS_BUCKET_PREFIX = "s3://cytof-repfam/$(SIMNAME)"

# NOTE: read configs
include("settings.jl")

for setting in settings
  setting[:results_dir] = "$(RESULTS_DIR_PREFIX)/$(setting[:outdir_suffix])"
  setting[:aws_bucket] = "$(AWS_BUCKET_PREFIX)/$(setting[:outdir_suffix])"
end

println("Sourcing files ..."); flush(stdout)
include("simfn.jl")

function sim(setting)
  results_dir = setting[:results_dir]
  aws_bucket = setting[:aws_bucket]
  mkpath(results_dir)

  println("Running $(results_dir)")
  flush(stdout)

  CytofRepFAM.Model.redirect_stdout_to_file("$(results_dir)/log.txt") do
    simfn(setting)
  end

  # Send to S3.
  CytofRepFAM.Model.s3sync(from=results_dir, to=aws_bucket,
                           tags=`--exclude '*.nfs'`)

  return
end

# NOTE:
# - f::Function: A function which takes one argument of type Dict{Any}
# - settings::Vector{Dict}): A vector of settings
println("Starting jobs ..."); flush(stdout)
_ = sim(settings[1])
println("DONE with all runs!")


text:
using Test
using Printf: @sprintf
using LinearAlgebra: dot
using StaticArrays: SVector
using QuantumLattices.Essentials.Terms
using QuantumLattices.Essentials.Spatials: Point, PID, Bond, Bonds, Lattice, decompose, acrossbonds, zerothbonds
using QuantumLattices.Essentials.DegreesOfFreedom: IDFConfig, Table, IID, Index, Internal, FilteredAttributes, OID, Operator, Operators, coordabsent, coordpresent, Boundary
using QuantumLattices.Interfaces: rank, update!, expand!, id, kind, reset!
using QuantumLattices.Prerequisites: Float, decimaltostr
using QuantumLattices.Prerequisites.CompositeStructures: NamedContainer
using QuantumLattices.Mathematics.AlgebraOverFields: ID, SimpleID, idtype
import QuantumLattices.Interfaces: dimension, expand
import QuantumLattices.Essentials.DegreesOfFreedom: isHermitian, otype
import QuantumLattices.Essentials.Terms: couplingcenter, couplingcenters, abbr

struct TID <: IID nambu::Int end
Base.adjoint(sl::TID) = TID(3-sl.nambu)

struct TIndex{S} <: Index{PID{S}, TID}
    scope::S
    site::Int
    nambu::Int
end
Base.fieldnames(::Type{<:TIndex}) = (:scope, :site, :nambu)
Base.union(::Type{P}, ::Type{I}) where {P<:PID, I<:TID} = TIndex{fieldtype(P, :scope)}

struct TFock <: Internal{TID} end
dimension(f::TFock) = 2
Base.getindex(f::TFock, i::Int) = TID(i)

struct TCID <: SimpleID
    center::Int
    nambu::Int
end
Base.fieldnames(::Type{<:TCID}) = (:center, :nambu)

struct TCoupling{V<:Number, I<:ID{TCID}} <: Coupling{V, I}
    value::V
    id::I
end
Base.repr(tc::TCoupling) = @sprintf "%s ph(%s)@(%s)" decimaltostr(tc.value) join(tc.id.nambus, ", ") join(tc.id.centers, ", ")
function expand(tc::TCoupling, pids::NTuple{R, PID}, focks::NTuple{R, TFock}, kind::Union{Val{K}, Nothing}=nothing) where {R, K}
    nambus = tc.id.nambus
    pids = NTuple{rank(tc), eltype(pids)}(pids[tc.id.centers[i]] for i = 1:rank(tc))
    return ((tc.value, NTuple{rank(tc), TIndex{fieldtype(pids|>eltype, :scope)}}(TIndex(pids[i].scope, pids[i].site, nambus[i]) for i = 1:rank(tc))),)
end

struct TOperator{V<:Number, I<:ID{OID}} <: Operator{V, I}
    value::V
    id::I
end

abbr(::Type{<:Term{ST, :TermMu} where ST}) = :tmu
isHermitian(::Type{<:Term{ST, :TermMu} where ST}) = true
otype(T::Type{<:Term{ST, :TermMu} where ST}, I::Type{<:OID}) = TOperator{T|>valtype, ID{I, T|>rank}}

abbr(::Type{<:Term{ST, :TermHopping} where ST}) = :thp
isHermitian(::Type{<:Term{ST, :TermHopping} where ST}) = false
otype(T::Type{<:Term{ST, :TermHopping} where ST}, I::Type{<:OID}) = TOperator{T|>valtype, ID{I, T|>rank}}

@testset "Subscript" begin
    sub = @subscript (x1, 4, 4, x2) with x1 < x2
    @test sub == deepcopy(sub)
    @test isequal(sub, deepcopy(sub))
    @test rank(sub) == rank(typeof(sub)) == 2
    @test dimension(sub) == dimension(typeof(sub)) == 4
    @test string(sub) == "(x1, 4, 4, x2) with $(sub.identifier)"
    @test sub(Val('M'), 1, 2) == (1, 4, 4, 2)
    @test sub(Val('C'), 1, 2) == true
    @test sub(Val('C'), 2, 1) == false

    sub = @subscript (x1, x2, x1, x2)
    @test sub(Val('M'), 1, 2) == (1, 2, 1, 2)
    @test sub(Val('C'), 1, 2) == true

    sub = Subscript{4}()
    @test rank(sub) == 1 && dimension(sub) == 4
    @test string(sub) == "(*, *, *, *)"
    @test sub(Val('M'), 2) == (2, 2, 2, 2)
    @test sub(Val('C'), 2) == true

    sub = Subscript((1, 2, 2, 1))
    @test rank(sub) == 0 && dimension(sub) == 4
    @test string(sub) == "(1, 2, 2, 1)"
    @test sub(Val('M')) == (1, 2, 2, 1)
    @test sub(Val('C')) == true
end

@testset "Subscripts" begin
    sub1 = @subscript (x1, 2) with x1 < 2
    sub2 = @subscript (y1, y2, y1, y2) with y1 < y2
    subs = Subscripts(sub1, sub2)
    @test rank(subs) == rank(typeof(subs)) == 3
    @test rank(subs, 1) == rank(typeof(subs), 1) == 1
    @test rank(subs, 2) == rank(typeof(subs), 2) == 2
    @test dimension(subs) == dimension(typeof(subs)) == 6
    @test dimension(subs, 1) == dimension(typeof(subs), 1) == 2
    @test dimension(subs, 2) == dimension(typeof(subs), 2) == 4
    @test subs(Val('M'), (1, 2, 3)) == (1, 2, 2, 3, 2, 3)
    @test subs(Val('C'), (1, 2, 3)) == true
    @test subs(Val('C'), (2, 2, 3)) == false
    @test subs(Val('C'), (1, 3, 2)) == false
    @test subs(Val('C'), (2, 3, 2)) == false

    sub3 = Subscript((6, 6))
    @test subs == sub1 * sub2
    @test subs * sub3 == Subscripts(sub1, sub2, sub3)
    @test sub3 * subs == Subscripts(sub3, sub1, sub2)
    @test subs * Subscripts(sub3) == Subscripts(sub1, sub2, sub3)
    @test expand(subs * sub3, (3, 3, 2, 3, 2, 3, 8, 9))|>collect == [(1, 2, 1, 2, 1, 2, 6, 6), (1, 2, 1, 3, 1, 3, 6, 6), (1, 2, 2, 3, 2, 3, 6, 6)]
end

@testset "centerrelated" begin
    @test couplingcenter(Coupling, 2, 2, Val(1)) == 1
    @test couplingcenter(Coupling, 1, 2, Val(1)) == 1
    @test couplingcenters(Coupling, ('*', '*'), Val(1)) == (1, 1)
    @test couplingcenters(Coupling, (1, 2, 3), Val(3)) == (1, 2, 3)
    @test couplingcenters("(1, 2, 3, 4)") == (1, 2, 3, 4)
end

@testset "couplings" begin
    @test @couplings(TCoupling(1.0, ID(TCID(1, 1), TCID(2, 1)))) == Couplings(TCoupling(1.0, ID(TCID(1, 1), TCID(2, 1))))
    @test @couplings(TCoupling(1.0, ID(TCID(1, 1), TCID(2, 1)))+TCoupling(1.0, ID(TCID(1, 2), TCID(2, 2)))) == TCoupling(1.0, ID(TCID(1, 1), TCID(2, 1)))+TCoupling(1.0, ID(TCID(1, 2), TCID(2, 2)))
end

@testset "TermFunction" begin
    ta = TermAmplitude()
    @test ta() == 1

    ta = TermAmplitude(x->x+3.0)
    @test ta(1.0) == 4.0

    tcs = TCoupling(1.0, ID(TCID(1, 1), TCID(2, 1))) + TCoupling(2.0, ID(TCID(1, 2), TCID(2, 2)))
    termcouplings = TermCouplings(tcs)
    @test termcouplings == deepcopy(TermCouplings(tcs))
    @test isequal(termcouplings, deepcopy(TermCouplings(tcs)))
    @test termcouplings() == tcs

    tcs1 = TCoupling(1.0, ID(TCID(1, 1), TCID(2, 1))) + TCoupling(1.0, ID(TCID(1, 2), TCID(2, 2)))
    tcs2 = TCoupling(1.0, ID(TCID(1, 2), TCID(2, 1))) + TCoupling(1.0, ID(TCID(2, 1), TCID(1, 2)))
    termcouplings = TermCouplings(i->(tcs1, tcs2)[(i-1)%2+1])
    @test termcouplings(1) == tcs1
    @test termcouplings(2) == tcs2

    termmodulate = TermModulate(:t)
    @test termmodulate(t=1) == 1
    @test termmodulate(mu=1) == nothing
    @test termmodulate() == nothing

    termmodulate = TermModulate(:t, t->t*2.0)
    @test termmodulate(2) == 4
end

@testset "Term" begin
    point = Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0))
    config = IDFConfig{TFock}(pid->TFock(), [PID(1, 1)])
    term = Term{'F', :TermMu, 2}(:mu, 1.5, 0, couplings=TCoupling(1.0, ID(TCID(1, 2), TCID(1, 1))), amplitude=(bond->3), modulate=false)
    @test term|>statistics == term|>typeof|>statistics == 'F'
    @test term|>kind == term|>typeof|>kind == :TermMu
    @test term|>id == term|>typeof|>id == :mu
    @test term|>valtype == term|>typeof|>valtype == Float
    @test term|>rank == term|>typeof|>rank == 2
    @test term|>abbr == term|>typeof|>abbr == :tmu
    @test term|>isHermitian == term|>typeof|>isHermitian == true
    @test term == deepcopy(term)
    @test isequal(term, deepcopy(term))
    @test string(term) == "TermMu{2F}(id=mu, value=1.5, bondkind=0, factor=1.0)"
    @test repr(term, point, config) == "tmu: 4.5 ph(2, 1)@(1, 1)"
    @test +term == term
    @test -term == term*(-1) == replace(term, factor=-term.factor)
    @test 2*term == term*2 == replace(term, factor=term.factor*2)
    @test term/2 == term*(1/2) == replace(term, factor=term.factor/2)

    p1 = Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0))
    p2 = Point(PID(1, 2), (0.0, 0.0), (0.0, 0.0))
    config = IDFConfig{TFock}(pid->TFock(), [PID(1, 1), PID(1, 2)])
    tcs1 = Couplings(TCoupling(1.0, ID(TCID(1, 2), TCID(1, 2))))
    tcs2 = Couplings(TCoupling(1.0, ID(TCID(1, 1), TCID(1, 1))))
    term = Term{'F', :TermMu, 2}(:mu, 1.5, 0, couplings=bond->(tcs1, tcs2)[bond.pid.site%2+1], amplitude=bond->3, modulate=true)
    @test repr(term, p1, config) == "tmu: 4.5 ph(1, 1)@(1, 1)"
    @test repr(term, p2, config) == "tmu: 4.5 ph(2, 2)@(1, 1)"
    @test one(term) == replace(term, value=1.0)
    @test zero(term) == replace(term, value=0.0)
    @test term.modulate(mu=4.0) == 4.0
    @test term.modulate(t=1.0) == nothing
    @test update!(term, mu=4.25) == replace(term, value=4.25)
    @test term.value == 4.25
end

@testset "expand" begin
    point = Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0))
    config = IDFConfig{TFock}(pid->TFock(), [PID(1, 1)])
    table = Table(config, filter(attr->(attr ≠ :nambu), FilteredAttributes(TIndex)))
    term = Term{'F', :TermMu, 2}(:mu, 1.5, 0, couplings=TCoupling(1.0, ID(TCID(1, 2), TCID(1, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(TOperator(+2.25, (TIndex(1, 1, 2), TIndex(1, 1, 1)), seqs=(1, 1)))
    @test expand(term, point, config, table, true, coordabsent) == operators
    @test expand(term, point, config, table, false, coordabsent) == operators*2

    bond = Bond(1, Point(PID('b', 2), (1.5, 1.5), (1.0, 1.0)), Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0)))
    config = IDFConfig{TFock}(pid->TFock(), [PID('a', 1), PID('b', 2)])
    table = Table(config, filter(attr->(attr ≠ :nambu), FilteredAttributes(TIndex)))
    term = Term{'F', :TermHopping, 2}(:t, 1.5, 1, couplings=TCoupling(1.0, ID(TCID(1, 2), TCID(2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(TOperator(4.5, (TIndex('a', 1, 2), TIndex('b', 2, 1)), rcoords=(SVector(0.5, 0.5), SVector(1.5, 1.5)), icoords=(SVector(0.0, 0.0), SVector(1.0, 1.0)), seqs=(1, 2)))
    @test expand(term, bond, config, table, true) == operators
    @test expand(term, bond, config, table, false) == operators+operators'

    lattice = Lattice("Tuanzi", [Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0))], vectors=[[1.0, 0.0]], neighbors=1)
    bonds = Bonds(lattice)
    config = IDFConfig{TFock}(pid->TFock(), lattice.pids)
    table = Table(config, filter(attr->(attr ≠ :nambu), FilteredAttributes(TIndex)))
    term = Term{'F', :TermMu, 2}(:mu, 1.5, 0, couplings=TCoupling(1.0, ID(TCID(1, 2), TCID(1, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(TOperator(+2.25, (TIndex(1, 1, 2), TIndex(1, 1, 1)), seqs=(1, 1)))
    @test expand(term, bonds, config, table, true, coordabsent) == operators
    @test expand(term, bonds, config, table, false, coordabsent) == operators*2
end

@testset "Parameters" begin
    ps1 = Parameters{(:t1, :t2, :U)}(1.0im, 1.0, 2.0)
    ps2 = Parameters{(:t1, :U)}(1.0im, 2.0)
    ps3 = Parameters{(:t1, :U)}(1.0im, 2.1)
    @test match(ps1, ps2) == true
    @test match(ps1, ps3) == false
end

@testset "Generator" begin
    lattice = Lattice("Tuanzi", [Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0)), Point(PID(1, 2), (0.5, 0.0), (0.0, 0.0))], vectors=[[1.0, 0.0]], neighbors=1)
    bonds = Bonds(lattice)
    config = IDFConfig{TFock}(pid->TFock(), lattice.pids)
    table = Table(config, filter(attr->(attr ≠ :nambu), FilteredAttributes(TIndex)))
    boundary = Boundary()
    t = Term{'F', :TermHopping, 2}(:t, 2.0, 1, couplings=TCoupling(1.0, ID(TCID(1, 2), TCID(2, 1))))
    μ = Term{'F', :TermMu, 2}(:μ, 1.0, 0, couplings=TCoupling(1.0, ID(TCID(1, 2), TCID(1, 1))), modulate=true)
    tops1 = expand(t, filter(acrossbonds, bonds, Val(:exclude)), config, table, true, coordpresent)
    tops2 = expand(one(t), filter(acrossbonds, bonds, Val(:include)), config, table, true, coordpresent)
    μops = expand(one(μ), filter(zerothbonds, bonds, Val(:include)), config, table, true, coordpresent)

    optp = TOperator{Float, ID{OID{TIndex{Int}, SVector{2, Float}, SVector{2, Float}, Int}, 2}}
    genops = GenOperators(tops1, NamedContainer{(:μ,)}((μops,)), NamedContainer{(:t, :μ)}((tops2, Operators{optp|>idtype, optp}())))
    @test genops == deepcopy(genops) && isequal(genops, deepcopy(genops))
    @test genops == GenOperators((t, μ), bonds, config, table, true, coordpresent)
    @test genops|>eltype == genops|>typeof|>eltype == optp
    @test genops|>idtype == genops|>typeof|>idtype == optp|>idtype
    @test expand!(Operators{idtype(optp), optp}(), genops, boundary, t=2.0, μ=1.5) == tops1+tops2*2.0+μops*1.5
    @test empty!(deepcopy(genops)) == empty(genops) == GenOperators(empty(μops), NamedContainer{(:μ,)}((empty(μops),)), NamedContainer{(:t, :μ)}((empty(μops), empty(μops))))
    @test reset!(deepcopy(genops), (t, μ), bonds, config, table, true, coordpresent) == genops

    gen = Generator((t, μ), bonds, config, table, true, boundary)
    @test gen == deepcopy(gen) && isequal(gen, deepcopy(gen))
    @test Parameters(gen) == Parameters{(:t, :μ)}(2.0, 1.0)
    @test expand!(Operators{idtype(optp), optp}(), gen) == expand(gen) == tops1+tops2*2.0+μops
    @test expand(gen, :t) == tops1+tops2*2.0
    @test expand(gen, :μ) == μops
    @test expand(gen, 1)+expand(gen, 2)+expand(gen, 3)+expand(gen, 4) == expand(gen)
    @test expand(gen, :μ, 1)+expand(gen, :μ, 2) == μops
    @test expand(gen, :t, 3) == tops1
    @test expand(gen, :t, 4) == tops2*2.0
    @test empty!(deepcopy(gen)) == Generator{coordpresent}((t, μ), empty(bonds), empty(config), empty(table), true, boundary) == empty(gen)
    @test reset!(empty(gen), lattice) == gen
    @test update!(gen, μ=1.5)|>expand == tops1+tops2*2.0+μops*1.5
end


text:
# Autogenerated wrapper script for libatomic_ops_jll for i686-w64-mingw32
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "PATH"
LIBPATH_default = ""

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ';')
    global LIBPATH = join(vcat(LIBPATH_list, [Sys.BINDIR, joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ';')

end  # __init__()


text:
@test Benitor.fibonacci(0) == 0
@test Benitor.fibonacci(1) == 1
@test Benitor.fibonacci(2) == 1
@test Benitor.fibonacci(4) == 3


text:
tests = [
    "uniloss",
    "multiloss",
    "prediction",
    "risks",
    "regularizers"
]

for t in tests
    tfile = string(t, ".jl")
    println("  * $tfile ...")
    include(tfile)
end


text:
# Something it seems I need - from https://github.com/JuliaLang/julia/issues/32567
import Base: redirect_stdout
function redirect_both(f::Function, stdout_buf, stderr_buf)
    old_stdout = stdout
    old_stderr = stderr

    stdout_rd,stdout_wr = redirect_stdout()
    stderr_rd,stderr_wr = redirect_stderr()

    try
        ret = f()
        return ret
    finally
        # This is ridiculous - readavailable will block without any output and
        # there seems to be no way to tell whether there is any output from the
        # pipes. Is it hidden in libuv somewhere?
        print(stdout, "\0")
        print(stderr, "\0")
        flush(stdout_wr)
        flush(stderr_wr)
        Libc.flush_cstdio()
        flush(stdout)
        flush(stderr)
        close(stdout_wr)
        close(stderr_wr)

        # # Throw away the extra null added.
        write(stdout_buf, readavailable(stdout_rd)[begin:end-1])
        write(stderr_buf, readavailable(stderr_rd)[begin:end-1])

        redirect_stdout(old_stdout)
        redirect_stderr(old_stderr)
    end
end

function PrettyTime(secs, sigfigs=2)
    time = secs*1e6

    factors = [1000, 1000, 60, 60, 24, Inf]
    cumfactors = cumprod(factors)
    names = ["μs", "ms", "secs", "mins", "hrs", "days"]
    ind = findfirst(time .< cumfactors)

    out = time
    if ind > 1
        out /= cumfactors[ind-1]
    end
    out = round(out, sigdigits=sigfigs)

    return string(out) * " " * names[ind]
end


text:
export StandardDiffusion

mutable struct StandardDiffusion <: DiffusionComponent
    "Scale for linear term."
    α::Float
    "Scale for quadratic term."
    β::Float
    "Denominator term for avoiding division by zero."
    ε::Float
    "Squared characteristic sound speed."
    cs²::Float

    function StandardDiffusion(; α = 1.0, β = 2.0, ε = 1e-3)
        @check α ≥ 0
        @check β ≥ 0
        @check ε ≥ 0
        new(α, β, ε, 0.0)
    end
end

function updatediffusion!(
    standard_diffusion::StandardDiffusion,
    energy_component::EnergyComponent,
    eos::EquationOfState,
)
    standard_diffusion.cs² =
        soundspeed²_from(eos, mean(get_specific_energies(energy_component)))
end

function compute_Δa_diffusion(
    r::Number,
    vᵣ::Number,
    hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
    ρᵢ_ρⱼ::Tuple{Number,Number},
    ddr_W_r_hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
    mass_component::MassComponent,
    standard_diffusion::StandardDiffusion,
)
    m = get_particle_mass(mass_component)
    α = standard_diffusion.α
    β = standard_diffusion.β
    ε = standard_diffusion.ε
    cs² = standard_diffusion.cs²

    hᵢⱼ = mean(hᵢ_hⱼ)
    ρᵢⱼ = mean(ρᵢ_ρⱼ)
    ddr_W_r_hᵢⱼ = mean(ddr_W_r_hᵢ_hⱼ)
    μᵢⱼ = vᵣ * r * hᵢⱼ / (r^2 + ε * hᵢⱼ^2)
    Πᵢⱼ = μᵢⱼ < 0 ? (β * μᵢⱼ - α * cs²) * μᵢⱼ / ρᵢⱼ : 0.0
    m * Πᵢⱼ * ddr_W_r_hᵢⱼ
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
import Cosmology.age
using Roots

"""
    age(h::AbstractGadgetHeader)

Computes the age of the universe given the properties from `AbstractGadgetHeader`.
"""
function age(h::AbstractGadgetHeader, units::Bool=true)

    c = cosmology(h)
    t = age(c, h.z)

    if units
        return t
    else
        return t |> ustrip
    end
end

"""
    redshift(c::AbstractCosmology, t_Gyrs::Real)

Computes the redshift at a given age of the universe in Gyrs for a cosmology `t`.
"""
function redshift(t_Gyrs::Real, c::Cosmology.AbstractCosmology)

    age_helper(z) = ustrip(age(c, z)) - t_Gyrs

    return find_zero(age_helper, [ Inf , 0.0] )
end

"""
    redshift(t::Unitful.AbstractQuantity, h::AbstractGadgetHeader)

See [`redshift`](@ref).
"""
redshift(t::Unitful.AbstractQuantity, h::AbstractGadgetHeader) = redshift( ustrip( t |> u"Gyr" ), cosmology(h))

"""
    redshift(t::Real, h::AbstractGadgetHeader)

See [`redshift`](@ref).
"""
redshift(t::Real, h::AbstractGadgetHeader) = redshift( t, cosmology(h))

"""
    redshift(t::Unitful.AbstractQuantity, h::AbstractGadgetHeader)

See [`redshift`](@ref).
"""
redshift(t::Unitful.AbstractQuantity, c::Cosmology.AbstractCosmology) = redshift( ustrip( t |> u"Gyr" ), c )


text:
## Sequence space defintions

# A Fun for SequenceSpace can be an iterator
Base.start(::Fun{SequenceSpace}) = 1
Base.next(f::Fun{SequenceSpace},st) = f[st],st+1
Base.done(f::Fun{SequenceSpace},st) = false # infinite length

getindex(f::Fun{SequenceSpace},k::Integer) =
    k ≤ ncoefficients(f) ? f.coefficients[k] : zero(eltype(f))
getindex(f::Fun{SequenceSpace},K) = eltype(f)[f[k] for k in K]

dotu(f::Fun{SequenceSpace},g::Fun{SequenceSpace}) =
    mindotu(f.coefficients,g.coefficients)
dotu(f::Fun{SequenceSpace},g::AbstractVector) =
    mindotu(f.coefficients,g)
dotu(f::AbstractVector,g::Fun{SequenceSpace}) =
    mindotu(f,g.coefficients)

Base.norm(f::Fun{SequenceSpace}) = norm(f.coefficients)
Base.norm(f::Fun{SequenceSpace},k::Int) = norm(f.coefficients,k)
Base.norm(f::Fun{SequenceSpace},k::Number) = norm(f.coefficients,k)

## Constant space defintions

Fun(c::Number)=Fun([c],ConstantSpace())
Fun(c::Number,d::ConstantSpace)=Fun([c],d)

dimension(::ConstantSpace) = 1

#TODO: Change
setdomain{CS<:AnyDomain}(f::Fun{CS},d::Domain) = Number(f)*ones(d)

canonicalspace(C::ConstantSpace)=C
spacescompatible(a::ConstantSpace,b::ConstantSpace)=domainscompatible(a,b)

Base.ones(S::ConstantSpace)=Fun(ones(1),S)
Base.ones(S::Union{AnyDomain,UnsetSpace})=ones(ConstantSpace())
Base.zeros(S::Union{AnyDomain,UnsetSpace})=zeros(ConstantSpace())
evaluate(f::AbstractVector,::ConstantSpace,x...)=f[1]
evaluate(f::AbstractVector,::ConstantSpace,x::Array)=f[1]*ones(x)

evaluate(f::AbstractVector,::ZeroSpace,x...)=zero(eltype(f))
evaluate(f::AbstractVector,::ZeroSpace,x::Array)=zeros(x)

Base.convert{CS<:ConstantSpace,T<:Number}(::Type{T},f::Fun{CS}) =
    convert(T,f.coefficients[1])

# promoting numbers to Fun
# override promote_rule if the space type can represent constants
Base.promote_rule{CS<:ConstantSpace,T<:Number}(::Type{Fun{CS}},::Type{T}) = Fun{CS,T}
Base.promote_rule{CS<:ConstantSpace,T<:Number,V}(::Type{Fun{CS,V}},::Type{T}) =
    Fun{CS,promote_type(T,V)}
Base.promote_rule{T<:Number,IF<:Fun}(::Type{IF},::Type{T}) = Fun

# When the union of A and B is a ConstantSpace, then it contains a one
conversion_rule(A::ConstantSpace,B::UnsetSpace)=NoSpace()
conversion_rule(A::ConstantSpace,B::Space)=(union_rule(A,B)==B||union_rule(B,A)==B)?A:NoSpace()

conversion_rule(A::ZeroSpace,B::Space) = A
maxspace_rule(A::ZeroSpace,B::Space) = B

Conversion(A::ZeroSpace,B::ZeroSpace) = ConversionWrapper(ZeroOperator(A,B))
Conversion(A::ZeroSpace,B::Space) = ConversionWrapper(ZeroOperator(A,B))

union_rule(A::ConstantSpace,B::Space)=ConstantSpace(domain(B))⊕B

## Special Multiplication and Conversion for constantspace

#  TODO: this is a special work around but really we want it to be blocks
Conversion{T,D}(a::ConstantSpace,b::Space{T,D,2})=ConcreteConversion{typeof(a),typeof(b),
        promote_type(op_eltype_realdomain(a),eltype(op_eltype_realdomain(b)))}(a,b)

Conversion(a::ConstantSpace,b::Space)=ConcreteConversion(a,b)
bandinds{CS<:ConstantSpace,S<:Space}(C::ConcreteConversion{CS,S})=1-ncoefficients(ones(rangespace(C))),0
function getindex{CS<:ConstantSpace,S<:Space,T}(C::ConcreteConversion{CS,S,T},k::Integer,j::Integer)
    if j != 1
        throw(BoundsError())
    end
    on=ones(rangespace(C))
    k ≤ ncoefficients(on)?T(on.coefficients[k]):zero(T)
end

# this is identity operator, but we don't use MultiplicationWrapper to avoid
# ambiguity errors

defaultMultiplication{CS<:ConstantSpace}(f::Fun{CS},b::ConstantSpace) =
    ConcreteMultiplication(f,b)
defaultMultiplication{CS<:ConstantSpace}(f::Fun{CS},b::Space) =
    ConcreteMultiplication(f,b)
defaultMultiplication(f::Fun,b::ConstantSpace) = ConcreteMultiplication(f,b)

bandinds{CS1<:ConstantSpace,CS2<:ConstantSpace,T}(D::ConcreteMultiplication{CS1,CS2,T}) =
    0,0
getindex{CS1<:ConstantSpace,CS2<:ConstantSpace,T}(D::ConcreteMultiplication{CS1,CS2,T},k::Integer,j::Integer) =
    k==j==1?T(D.f.coefficients[1]):one(T)

rangespace{CS1<:ConstantSpace,CS2<:ConstantSpace,T}(D::ConcreteMultiplication{CS1,CS2,T}) =
    D.space

rangespace{F<:ConstantSpace,T}(D::ConcreteMultiplication{F,UnsetSpace,T}) =
    UnsetSpace()
bandinds{F<:ConstantSpace,T}(D::ConcreteMultiplication{F,UnsetSpace,T}) =
    error("No range space attached to Multiplication")
getindex{F<:ConstantSpace,T}(D::ConcreteMultiplication{F,UnsetSpace,T},k::Integer,j::Integer) =
    error("No range space attached to Multiplication")

bandinds{CS<:ConstantSpace,F<:Space,T}(D::ConcreteMultiplication{CS,F,T}) = 0,0
getindex{CS<:ConstantSpace,F<:Space,T}(D::ConcreteMultiplication{CS,F,T},k::Integer,j::Integer) =
    k==j?T(D.f):zero(T)
rangespace{CS<:ConstantSpace,F<:Space,T}(D::ConcreteMultiplication{CS,F,T}) = D.space

bandinds{CS<:ConstantSpace,F<:Space,T}(D::ConcreteMultiplication{F,CS,T}) = 1-ncoefficients(D.f),0
function getindex{CS<:ConstantSpace,F<:Space,T}(D::ConcreteMultiplication{F,CS,T},k::Integer,j::Integer)
    Op = Multiplication(D.f,space(D.f))
    k≤ncoefficients(D.f) && j==1?T(Op[k,1]):zero(T)
end
rangespace{CS<:ConstantSpace,F<:Space,T}(D::ConcreteMultiplication{F,CS,T}) =
    rangespace(Multiplication(D.f,space(D.f)))

# functionals always map to Constant space
function promoterangespace(P::Operator,A::ConstantSpace,cur::ConstantSpace)
    @assert isafunctional(P)
    domain(A)==domain(cur)?P:SpaceOperator(P,domainspace(P),A)
end

for op = (:*,:.*,:./,:/)
    @eval $op{CS<:ConstantSpace}(f::Fun,c::Fun{CS}) = f*convert(Number,c)
end

## Multivariate case
union_rule(a::TensorSpace,b::ConstantSpace{AnyDomain})=TensorSpace(map(sp->union(sp,b),a.spaces))
## Special spaces

function Base.convert{TS<:TensorSpace,T<:Number}(::Type{T},f::Fun{TS})
    if all(sp->isa(sp,ConstantSpace),space(f).spaces)
        convert(T,f.coefficients[1])
    else
        error("Cannot convert $f to type $T")
    end
end

Base.convert{CS1<:ConstantSpace,CS2<:ConstantSpace,T<:Number,TT,d}(::Type{T},f::Fun{TensorSpace{Tuple{CS1,CS2},TT,d}}) =
    convert(T,f.coefficients[1])

isconstspace(sp::TensorSpace) = all(isconstspace,sp.spaces)


text:
module Diagnostics

abstract type Format end
abstract type Lint end
abstract type Action end

mutable struct Diagnostic{C}
    loc::UnitRange
    actions::Vector{Action}
    message::String
end
Diagnostic(r::UnitRange) = Diagnostic(r, [], "")

@enum(ErrorCodes,
ParseFailure)

@enum(FormatCodes,
Useelseif,
Indents,
CamelCase,
LowerCase,
MissingWS,
ExtraWS,
CommaWS)

@enum(LintCodes,
DuplicateArgumentName,
ArgumentFunctionNameConflict,
SlurpingPosition,
KWPosition,
ImportInFunction,
DuplicateArgument,
LetNonAssignment,
RangeNonAssignment,
CondAssignment,
DeadCode,
DictParaMisSpec,
DictGenAssignment,
MisnamedConstructor,
LoopOverSingle,
AssignsToFuncName,
PossibleTypo,

Deprecation,
functionDeprecation,
typeDeprecation,
immutableDeprecation,
abstractDeprecation,
bitstypeDeprecation,
typealiasDeprecation,
parameterisedDeprecation)

struct Deletion <: Action
    range::UnitRange
end

struct AddWS <: Action
    range::UnitRange
    length::Int
end

struct TextEdit <: Action
    range::UnitRange
    text::String
end

end

# Formatting

islbracket(t::Token) = t.kind == Tokens.LPAREN ||
                        t.kind == Tokens.LBRACE ||
                        t.kind == Tokens.LSQUARE

isrbracket(t::Token) = t.kind == Tokens.RPAREN ||
                        t.kind == Tokens.RBRACE ||
                        t.kind == Tokens.RSQUARE

function format_op(ps, prec)
    !ps.formatcheck && return
    loc = ps.t.startbyte:ps.t.endbyte + 1
    if (prec == ColonOp || prec == PowerOp || prec == DeclarationOp || prec == DotOp) && ps.t.kind != Tokens.ANON_FUNC && ps.t.kind != Tokens.PRIME
        if ps.lws.kind != EmptyWS && ps.ws.kind != EmptyWS
            diag = Diagnostic{Diagnostics.ExtraWS}(loc, Diagnostics.Action[Diagnostics.Deletion(ps.ws.startbyte:ps.nt.startbyte),Diagnostics.Deletion(ps.lws.startbyte:ps.t.startbyte)], "Unexpected white space around operator")
            push!(ps.diagnostics, diag)
        elseif ps.lws.kind != EmptyWS
            diag = Diagnostic{Diagnostics.ExtraWS}(loc, Diagnostics.Action[Diagnostics.Deletion(ps.lws.startbyte:ps.t.startbyte)], "Unexpected white space around operator")
            push!(ps.diagnostics, diag)
        elseif ps.ws.kind != EmptyWS
            diag = Diagnostic{Diagnostics.ExtraWS}(loc, Diagnostics.Action[Diagnostics.Deletion(ps.ws.startbyte:ps.nt.startbyte)], "Unexpected white space around operator")
            push!(ps.diagnostics, diag)
        end
    elseif ps.t.kind == Tokens.PRIME
        if ps.lws.kind != EmptyWS
            diag = Diagnostic{Diagnostics.ExtraWS}(loc, Diagnostics.Action[Diagnostics.Deletion(ps.ws.startbyte:ps.nt.startbyte),Diagnostics.Deletion(ps.lws.startbyte:ps.t.startbyte)], "Unexpected white space around operator")
            push!(ps.diagnostics, diag)
        end
    elseif ps.t.kind == Tokens.ISSUBTYPE || ps.t.kind == Tokens.DDDOT
    else
        if ps.lws.kind == EmptyWS && ps.ws.kind == EmptyWS
            diag = Diagnostic{Diagnostics.MissingWS}(loc, Diagnostics.Action[Diagnostics.AddWS(ps.nt.startbyte:ps.nt.startbyte, 1),Diagnostics.AddWS(ps.t.startbyte:ps.t.startbyte, 1)], "Missing white space around operator")
            push!(ps.diagnostics, diag)
        elseif ps.lws.kind == EmptyWS
            diag = Diagnostic{Diagnostics.MissingWS}(loc, Diagnostics.Action[Diagnostics.AddWS(ps.t.startbyte:ps.t.startbyte, 1)], "Missing white space around operator")
            push!(ps.diagnostics, diag)
        elseif ps.ws.kind == EmptyWS
            diag = Diagnostic{Diagnostics.MissingWS}(loc, Diagnostics.Action[Diagnostics.AddWS(ps.nt.startbyte:ps.nt.startbyte, 1)], "Missing white space around operator")
            push!(ps.diagnostics, diag)
        end
    end
end

function format_comma(ps)
    !ps.formatcheck && return
    if ps.lws.kind != EmptyWS && !(islbracket(ps.lt))
        push!(ps.diagnostics, Diagnostic{Diagnostics.ExtraWS}(ps.t.startbyte + (0:1), [Diagnostics.Deletion(ps.lws.startbyte:ps.t.startbyte)], "Unexpected white space preceding comma"))
    end
    if ps.ws.kind == EmptyWS && !(isrbracket(ps.nt))
        push!(ps.diagnostics, Diagnostic{Diagnostics.MissingWS}(ps.t.startbyte + (0:1), [Diagnostics.AddWS(ps.nt.startbyte:ps.nt.startbyte, 1)], "Missing white space following comma"))
    end
end

function format_lbracket(ps)
    !ps.formatcheck && return
    if ps.ws.kind == WS
        push!(ps.diagnostics, Diagnostic{Diagnostics.ExtraWS}(ps.t.startbyte + (0:1), [Diagnostics.Deletion(ps.ws.startbyte:ps.nt.startbyte)], "Unexpected white space following $(ps.t.val)"))
    end
end

function format_rbracket(ps)
    !ps.formatcheck && return
    if ps.lws.kind == WS
        push!(ps.diagnostics, Diagnostic{Diagnostics.ExtraWS}(ps.t.startbyte + (0:1), [Diagnostics.Deletion(ps.lws.startbyte:ps.t.startbyte)], "Unexpected white space preceding $(ps.t.val)"))
    end
end

function format_indent(ps, start_col)
    !ps.formatcheck && return
    if (start_col > 0 && ps.nt.startpos[2] != start_col)
        dindent = start_col - ps.nt.startpos[2]
        if dindent > 0
            push!(ps.diagnostics, Diagnostic{Diagnostics.MissingWS}(ps.nt.startbyte + (0:dindent), [Diagnostics.AddWS(ps.nt.startbyte:ps.nt.startbyte, dindent)], ""))
        else
            push!(ps.diagnostics, Diagnostic{Diagnostics.ExtraWS}(ps.nt.startbyte + (dindent + 1:0), [Diagnostics.Deletion(ps.nt.startbyte + (dindent:0))], ""))
        end
    end
end

function format_kw(ps)
    !ps.formatcheck && return
    if ps.ws.kind == WS
        if length(ps.ws.val) > 1
            push!(ps.diagnostics, Diagnostic{Diagnostics.ExtraWS}(ps.t.startbyte:ps.nt.startbyte, [Diagnostics.Deletion(ps.ws.startbyte + 1:ps.nt.startbyte)], ""))
        end
    end
end

function format_no_rws(ps)
    !ps.formatcheck && return
    if ps.ws.kind != EmptyWS
        push!(ps.diagnostics, Diagnostic{Diagnostics.ExtraWS}(ps.t.startbyte:ps.nt.startbyte, [Diagnostics.Deletion(ps.ws.startbyte:ps.nt.startbyte)], ""))
    end
end

function format_typename(ps, sig)
    !ps.formatcheck && return
#     start_loc = ps.nt.startbyte - sig.span
#     id = get_id(sig)
#     sig isa EXPR && return
#     val = string(id.val)
#     # Abitrary limit of 3 for uppercase acronym
#     if islower(first(val)) || (length(val) > 3 && all(isupper, val))
#         push!(ps.diagnostics, Diagnostic{Diagnostics.CamelCase}(start_loc + (1:sizeof(val))))
#     end
end

function format_funcname(ps, id, offset)
    !ps.formatcheck && return
    # start_loc = ps.nt.startbyte - offset
    # !(id isa Symbol) && return
    # val = string(id)
    # if !islower(val) #!all(islower(c) || isdigit(c) || c == '!' for c in val)
    #     push!(ps.diagnostics, Diagnostic{Diagnostics.LowerCase}(start_loc + (1:sizeof(val))))
    # end
end