text:
"""
    unroll(moving_average::Vector{Float64}, window::Int64; initial_conditions::U = nothing, assert_natural::Bool = false) where { U <: Union{ Tuple{Vararg{Union{Int64,Float64}}},Nothing} }

Retrive original time series (i.e. unroll) from its moving average `moving_average`.
# Arguments
- `moving_average::Vector{Float64}`: the time series representing the moving average to unroll ;
- `window:::Int64`: the width of the moving average ;
- `initial_conditions::U = nothing`: the initial values of the original time series to be recovered. It may be a `Tuple` of `window-1` float or integer values, or `nothing` if initial conditions are unknown ;
- `assert_natural::Bool = false` default boolean argument. If true, then the pipeline will try to recover a time series of natural numbers only. More then one acceptable time series (where "acceptable" means that it reproduces `moving_average`) may be found and returned .

NB: If ```isnothing(initial_conditions) && !assert_natural``` , then only an approximate method may be used, see this [StackExchange post](https://stats.stackexchange.com/a/68002).
"""
function unroll(moving_average::Vector{Float64}, window::Int64; initial_conditions::U = nothing, assert_natural::Bool = false) where { U <: Union{ Tuple{Vararg{Union{Int64,Float64}}},Nothing} }

    reconstructed_time_series = assert_natural ? Vector{Int64}[] : Vector{Float64}[]

    if isnothing(initial_conditions)
        if assert_natural
            reconstructed_time_series = unroll_iterative(moving_average, window)
        else
            push!(reconstructed_time_series, unroll_linear_approximation(moving_average, window))
        end
    elseif length(initial_conditions) == window-1
        push!(reconstructed_time_series, unroll_recursive(moving_average, window,initial_conditions))
    else
        error("`initial_conditions` type must be either Nothing or a NTuple{window-1,Union{Float64,Int64}}")
    end

end

"""
    unroll_iterative(moving_average::Vector{Float64}, window::Int64)

Unroll `moving_average` (interpreting it as a moving average whose window width is `window`), returning the original time series assuming it is composed of only natural numbers.

The methodology is as follows:
1. Consider the minimum of `moving_average`, that we will name `minimum_average` ;
2. Produce all possible sets of `n₋ + n₊ + 1` naturals that could have `minimum_average` as mean. This is perfomed by obtaining all permutations of all the partititions of `minimum_average*(n₋ + n₊ + 1)` via Combinatorics.jl, and we will refer to each of the resulting array as a "possibility". These are irganized in an array of iterators `possibilities`;
3. For each for each `possibility` and for each `element` in `moving_average[(minimum_index +1):end]` compute the natural `x` (zero included) to be pushed to `possibility`'s end so that `sum(possibility[(i+1):(i+1 + n₋ + n₊ - 1)]) == minimum_average*(n₋ + n₊ + 1)`. if there is such `x`, push it to the `possibility`'s end and go to the next possibility, else remove the `possibility`. When this loop finishes, perform the same loops backward for each `element in reverse(moving_average[1:(minimum_index - 1)])`, this time pushing the `x`s to the `possibility`'s top. This allows for obtaining the set of all possible time series ;
4. Return the remaining possibilities.
"""
function unroll_iterative(moving_average::Vector{Float64}, window::Int64) #  n₋::Int64, n₊::Int64, window = n₋ + n₊ + 1

    # compute linear approximation
    # mean_vec::Vector{Float64} = repeat([1/ window ], window ) #n₋ + n₊ + 1
    # mean_mat::Matrix{Float64} = hcat([pad_left_right(mean_vec, i, length(moving_average) - i - 1) for i in 0:(length(moving_average) - 1)]...)'
    # linear_approximation::Vector{Float64} = pinv(mean_mat) * moving_average #mean_mat \ moving_average

    # find minimum element of the moving average, and compute the numerator of the such average
    minimum_index::Int64 = argmin(moving_average)
    minimum_window_total_cases::Int64 = round(Int64,moving_average[minimum_index] * window) #n₋ + n₊ + 1
    #println("minimum_window_total_cases = ", minimum_window_total_cases)

    # if minimum_window_total_cases == 0, set collected_partitions manually since Combinatorics would return an "undefined reference"
    collected_partitions = Vector{Int64}[]
    if minimum_window_total_cases == 0
        collected_partitions = [[0]]
    # Else, compute it using Combinatorics.jl
    else
        collected_partitions = collect(partitions(minimum_window_total_cases))
    end

    # Get all permutations of partitions of such numerator. organize them as the column of a matrix.
    possibilities::Vector{Combinatorics.Permutations{Vector{Int64}}} = permutations.([pad!(partition,window) for partition in collected_partitions if length(partition) <= window ])

    # Initialize arrays that will contain the index of the columns (possibilities) to keep (`to_be_kept`) and the new row to be added after each iteration (see for loop below)
    to_be_kept = Set{Vector{Int64}}()

    # Pre-compute the numerators of all the averages
    numerators::Vector{Int64} = round.(Ref(Int64),moving_average .* window)

    #l::Int64 = length(possibilities)
    # loop over possibilities, and store the ones that reproduce the moving average
    for perm_it in possibilities
        # if k%1000 == 0
        #     println("doing $k \\ $l")
        # end
        for possibility in perm_it
            valid::Bool = true
            # go forward
            for (i,numerator) in enumerate(@view(numerators[(minimum_index +1):end]))
                diff::Int64 = numerator - sum(@view(possibility[(i+1):(i+1 + window -1 - 1)]))
                if diff < 0
                    valid = false
                    break
                else
                    push!(possibility,diff)
                end
            end

            # go backward
            if valid
                for numerator in reverse(@view(numerators[1:(minimum_index - 1)]))
                diff::Int64 = numerator - sum( @view(possibility[1:(window-1)]))
                if diff < 0
                    valid = false
                    break
                else
                    pushfirst!(possibility,diff)
                end
            end

            if valid
                push!(to_be_kept,possibility)
            end

        end

        end
    end
    # println("unroll_rolling_mean_of_natural_series_iterators. Returning...")
    # println("length(to_be_kept) = ", length(to_be_kept))
    # return accepted possibilities
    return collect(to_be_kept)

end

"""
    unroll_linear_approximation(moving_average::Vector{Float64}, window::Int64)

Compute the linear approximation of the time series whose moving average of window `window` is `moving_average`. For details, please refer to https://stats.stackexchange.com/a/68002 .
"""
function unroll_linear_approximation(moving_average::Vector{Float64}, window::Int64)

 # compute linear approximation
 mean_vec::Vector{Float64} = repeat([1/ window ], window )
 mean_mat::Matrix{Float64} = hcat([pad_left_right(mean_vec, i, length(moving_average) - i - 1) for i in 0:(length(moving_average) - 1)]...)'
 linear_approximation::Vector{Float64} = pinv(mean_mat) * moving_average # equivalently, mean_mat \ moving_average

 return linear_approximation

end

"""
    unroll_recursive(moving_average::Vector{Float64}, window::Int64 , initial_conditions::Tuple{Vararg{Union{Int64,Float64}}} )

Return the raw time series of which `moving_average` is the moving average, where `moving_average[i]` is the moving average of the `[(i-n₋):(i+n₊)]` slice of the raw time series to be returned. Assume the raw time series is composed of only natural numbers.

The methodology is as follow:
1. Initialize `deaveraged = collect(initial)` ;
2. for each `i ∈ eachindex(moving_average) `, set `deaveraged[i+window-1] = round(Int64,window*moving_average[i] - sum(@view(deaveraged[i:(i+window-1-1)])))` ;
3. Return `deaveraged` .
"""
function unroll_recursive(moving_average::Vector{Float64}, window::Int64 , initial_conditions::Tuple{Vararg{Union{Int64,Float64}}}) #NTuple{n₋+n₊,Int64} n₋::Int64, n₊::Int64

    # check that initial conditions are of the correct size
    if length(initial_conditions) != window - 1
        error("initial_conditions must have length equal to (window-1) = ", window-1)
    end
     # initialize raw time series using provided `initial_conditions`
    deaveraged::Vector{Union{Int64,Float64}} = collect(initial_conditions)
    # pre-resize `deaveraged` to accomodate all the reconstructed time series elements
    resize!(deaveraged, length(initial_conditions) + length(moving_average))

    # loop over moving_average
    for i ∈ eachindex(moving_average) #in 1:(length(moving_average))
        deaveraged[i+window-1] = round(Int64,window*moving_average[i] - sum(@view(deaveraged[i:(i+window-1-1)])))
    end

    # return deaveraged time series
    return deaveraged

end

#######################################

# function unroll_iterative(moving_average::Vector{Float64}, n₋::Int64, n₊::Int64) #  n₋::Int64, n₊::Int64, window = n₋ + n₊ + 1

#     # compute linear approximation
#     mean_vec::Vector{Float64} = repeat([1/(n₋ + n₊ + 1 )], n₋ + n₊ + 1 ) #n₋ + n₊ + 1
#     mean_mat::Matrix{Float64} = hcat([pad_left_right(mean_vec, i, length(moving_average) - i - 1) for i in 0:(length(moving_average) - 1)]...)'
#     linear_approximation::Vector{Float64} = pinv(mean_mat) * moving_average #mean_mat \ moving_average

#     # find minimum element of the moving average, and compute the numerator of the such average
#     minimum_index::Int64 = argmin(moving_average)
#     minimum_window_total_cases::Int64 = round(Int64,moving_average[minimum_index] * (n₋ + n₊ + 1 ))
#     println("minimum_window_total_cases = ", minimum_window_total_cases)

#     # compute the numeratir of the frst average
#     # first_window_total_cases::Int64 = round(Int64,moving_average[1] * (n₋ + n₊ + 1 ))
#     # if first_window_total_cases == 0, set collected_partitions manually since Combinatorics would return a "indefined reference"
#     #println("first_window_total_cases  = $first_window_total_cases ")
#     collected_partitions = Int64[]
#     if minimum_window_total_cases == 0
#         collected_partitions = [[0]]
#     # else, compute it using Combinatorics.jl
#     else
#         collected_partitions = collect(partitions(minimum_window_total_cases))
#     end

#     # get all permutations of partitions of such numerator. organize them as the column of a matrix.
#     possibilities::Vector{Combinatorics.Permutations{Vector{Int64}}} = permutations.([pad!(partition, n₋ + n₊ + 1) for partition in collected_partitions if length(partition) <= (n₋ + n₊ + 1) ])
#     #for partition in second_collected_partitions if length(partition) <= (n₋ + n₊ + 1) ])
#     # initialize arrays that will contain the index of the columns (possibilities) to keep (`to_be_kept`) and the new row to be added after each iteration (see for loop below)
#     to_be_kept = Set{Vector{Int64}}()

#     numerators::Vector{Int64} = round.(Ref(Int64),moving_average .* (n₋ + n₊ + 1))

#     #l::Int64 = length(possibilities)
#     for (k,perm_it) in enumerate(possibilities)
#         # if k%1000 == 0
#         #     println("doing $k \\ $l")
#         # end
#         for possibility in perm_it
#             valid::Bool = true
#             # go forward
#             for (i,numerator) in enumerate(@view(numerators[(minimum_index +1):end])) # enumerate(moving_average[(minimum_index +1):end]) #(n₋ + n₊ + 1)
#                 diff::Int64 =  numerator  - sum( @view(possibility[(i+1):(i+1 + n₋ + n₊ - 1)]) ) # round(Int64,average*(n₋ + n₊ + 1 ))
#                 if diff < 0
#                     valid = false
#                     break
#                 else
#                     push!(possibility,diff)
#                 end
#             end

#             # go backward
#             if valid
#                 for numerator in reverse(@view(numerators[1:(minimum_index - 1)])) #reverse(moving_average[1:(minimum_index - 1)])
#                 diff::Int64 = numerator - sum( @view(possibility[1:(n₋ + n₊)]) ) # round(Int64,average*(n₋ + n₊ + 1 ))
#                 if diff < 0
#                     valid = false
#                     break
#                 else
#                     pushfirst!(possibility,diff)
#                 end
#             end

#             if valid
#                 push!(to_be_kept,possibility) #[1:(n₋ + n₊ )] #@view(
#             end

#         end

#         end
#     end

#     println("unroll_rolling_mean_of_natural_series_iterators. Returning...")

#     # end

#     println("length(to_be_kept) = ", length(to_be_kept))

#     # For every initial conditions, which are all slices `possibility[1:(n₋ + n₊)]` for the remaining possibilities, apply `deaverage(moving_average, n₋, n₊, Tuple(possibility[1:(n₋ + n₊ )])  )` and keep the resulting reconstructed raw series iff it contains no negative numbers. Return all possible reconstructed time series, togheter with a linear approximation.
#     #return [pred_poss for pred_poss in deaverage.(Ref(moving_average), Ref(n₋), Ref(n₊), Tuple.(collect(to_be_kept))  ) if all(pred_poss .>= 0)] , linear_approximation

#     return collect(to_be_kept), linear_approximation

# end

###############################################################################


text:
using BunchKaufmanUtils, StaticArrays, LinearAlgebra
using Test

@testset "BunchKaufmanUtils.jl" begin
    θ = pi/6
    V = [cos(θ) -sin(θ); sin(θ) cos(θ)]
    D = Diagonal([-1, 1])
    A = V*D*V'
    F = bunchkaufman(A)
    @test !ispossemidef(F)
    v = [1,2]
    @test pseudosolve(F, v) ≈ F \ v ≈ A \ v
    @test F*v ≈ A*v
    AS = SMatrix{2,2}(A)
    vs = SVector{2}(v)
    FS = @inferred(bunchkaufman(AS))
    @test size(FS) == (2,2)
    @test size(FS, 1) == 2
    @test size(FS, 3) == 1
    @test @inferred(pseudosolve(FS, vs)) ≈ @inferred(FS \ vs) ≈ @inferred(AS \ vs)

    A = [1 1; 1 1]
    v = [2, 2]
    F = bunchkaufman(A; check=false)
    x = pseudosolve(F, v)
    @test A*x ≈ v

    # Trivial all-zeros case
    F = bunchkaufman(reshape([0], (1, 1)), check=false)
    @test pseudosolve(F, [0]) == [0]
end

# Tests modified from the stdlib Bunch-Kaufman tests
using Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted

@testset "stdlib tests" begin
    n = 10

    # Split n into 2 parts for tests needing two matrices
    n1 = div(n, 2)
    n2 = 2*n1

    Random.seed!(1234321)

    areal = randn(n,n)/2
    aimg  = randn(n,n)/2
    a2real = randn(n,n)/2
    a2img  = randn(n,n)/2
    breal = randn(n,2)/2
    bimg  = randn(n,2)/2
    vreal, vimg = randn(n), randn(n)

    @testset "$eltya argument A" for eltya in (Float32, Float64, Int) # , ComplexF32, ComplexF64)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
        vT = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(vreal, vimg) : vreal)
        asym = transpose(a) + a                  # symmetric indefinite
        aher = a' + a                  # Hermitian indefinite
        apd  = a' * a                  # Positive-definite
        for (a, a2, aher, apd) in ((a, a2, aher, apd),
                                   (view(a, 1:n, 1:n),
                                    view(a2, 1:n, 1:n),
                                    view(aher, 1:n, 1:n),
                                    view(apd , 1:n, 1:n)))
            ε = εa = eps(abs(float(one(eltya))))

            @testset "$uplo Bunch-Kaufman factor of indefinite matrix" for uplo in (:U,) # (:U, :L)
                bc1 = bunchkaufman(Hermitian(aher, :U))
                @test bc1 \ aher ≈ Matrix(I, n, n)
                @test pseudosolve(bc1, aher) ≈ Matrix(I, n, n)
                @test bc1*a ≈ aher*a
                @test bc1*vT ≈ aher*vT
                @testset for rook in (false, true)
                    @test pseudosolve(bunchkaufman(Symmetric(transpose(a) + a, uplo), rook), transpose(a) + a) ≈ Matrix(I, n, n)
                end
            end

            @testset "$eltyb argument B" for eltyb in (Float32, Float64, Int) # ,ComplexF32, ComplexF64)
                b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
                for b in (b, view(b, 1:n, 1:2))
                    εb = eps(abs(float(one(eltyb))))
                    ε = max(εa,εb)

                    @testset "$uplo Bunch-Kaufman factor of indefinite matrix" for uplo in (:U,) #(:L, :U)
                        bc1 = bunchkaufman(Hermitian(aher, uplo))
                        @test aher*(bc1\b) ≈ b atol=1000ε
                        @test aher*(pseudosolve(bc1, b)) ≈ b atol=1000ε
                    end

                    @testset "$uplo Bunch-Kaufman factors of a pos-def matrix" for uplo in (:U,) #(:U, :L)
                        @testset "rook pivoting: $rook" for rook in (false, true)
                            bc2 = bunchkaufman(Hermitian(apd, uplo), rook)
                            @test pseudosolve(bc2, apd) ≈ Matrix(I, n, n)
                            @test apd*(pseudosolve(bc2, b)) ≈ b rtol=eps(cond(apd))
                        end
                    end
                end
            end
        end
    end

    @testset "test example due to @timholy in PR 15354" begin
        A = rand(6,5); A = complex(A'*A) # to avoid calling the real-lhs-complex-rhs method
        F = bunchkaufman(A);
        v6 = rand(ComplexF64, 6)
        v5 = view(v6, 1:5)
        @test pseudosolve(F, v5) ≈ F\v6[1:5]
    end
end

ext: jl


text:
############################################################
# joPDE ####################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export joPDE, joPDEexception, joPDEadjoint, joPDEadjointException, subsample

############################################################

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections

# Forward
struct joPDE{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
	info::Info
	model::ModelAll
	geometry::Geometry
	options::Options
	fop::Function              # forward
	fop_T::Nullable{Function}  # transpose
end

mutable struct joPDEexception <: Exception
	msg :: String
end

# Adjoint
struct joPDEadjoint{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
	name::String
	m::Integer
	n::Integer
	info::Info
	model::ModelAll
	geometry::Geometry
	options::Options
	fop::Function              # forward
	fop_T::Nullable{Function}  # transpose
end

mutable struct joPDEadjointException <: Exception
	msg :: String
end

############################################################
## outer constructors

function joPDE(name::String,info::Info,model::ModelAll, geometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
	(DDT == Float32 && RDT == Float32) || throw(joPDEexception("Domain and range types not supported"))
	if typeof(geometry) == GeometryOOC
		m = sum(geometry.nsamples)
	else
		m = 0
		for j=1:info.nsrc m += length(geometry.xloc[j])*geometry.nt[j] end
	end
	n = info.n * sum(info.nt)
	if info.nsrc > 1
		srcnum = 1:info.nsrc
	else
		srcnum = 1
	end
	return F = joPDE{Float32,Float32}(name, m, n, info, model, geometry, options,
							  src -> time_modeling(model, src.geometry, src.data, geometry, [], [], srcnum, 'F', 1, options),
							  rec -> time_modeling(model, geometry, [], rec.geometry, rec.data, [], srcnum, 'F', -1, options)
							  )
end

function joPDEadjoint(name::String,info::Info,model::ModelAll, geometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
	(DDT == Float32 && RDT == Float32) || throw(joPDEadjointException("Domain and range types not supported"))
	if typeof(geometry) == GeometryOOC
		m = sum(geometry.nsamples)
	else
		m = 0
		for j=1:info.nsrc m += length(geometry.xloc[j])*geometry.nt[j] end
	end
	n = info.n * sum(info.nt)
	if info.nsrc > 1
		srcnum = 1:info.nsrc
	else
		srcnum = 1
	end
	return F = joPDEadjoint{Float32,Float32}(name, m, n, info, model, geometry, options,
									 rec -> time_modeling(model, geometry, [], rec.geometry, rec.data, [], srcnum, 'F', -1, options),
									 src -> time_modeling(model, src.geometry, src.data, geometry, [], [], srcnum, 'F', 1, options),
								     )
end

############################################################
## overloaded Base functions

# conj(joPDE)
conj{DDT,RDT}(A::joPDE{DDT,RDT}) =
	joPDE{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
	get(A.fop),
		A.fop_T
		)

# transpose(jo)
transpose{DDT,RDT}(A::joPDE{DDT,RDT}) =
	joPDEadjoint{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
		get(A.fop_T),
		A.fop
		)

# ctranspose(jo)
ctranspose{DDT,RDT}(A::joPDE{DDT,RDT}) =
	joPDEadjoint{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
		get(A.fop_T),
		A.fop
		)

# conj(jo)
conj{DDT,RDT}(A::joPDEadjoint{DDT,RDT}) =
    joPDEadjoint{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
        get(A.fop),
        A.fop_T
        )

# transpose(jo)
transpose{DDT,RDT}(A::joPDEadjoint{DDT,RDT}) =
	joPDE{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
		get(A.fop_T),
		A.fop
		)

# ctranspose(jo)
ctranspose{DDT,RDT}(A::joPDEadjoint{DDT,RDT}) =
	joPDE{DDT,RDT}(A.name,A.n,A.m,A.info,A.model,A.geometry,A.options,
		get(A.fop_T),
		A.fop
		)

############################################################
## overloaded Base *(...jo...)

# *(joPDE,joRHS)
function *{ADDT,ARDT,vDT}(A::joPDE{ADDT,ARDT},v::joRHS{vDT})
	A.n == size(v,1) || throw(joPDEexception("shape mismatch"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(joPDE,joRHS):",A.name,typeof(A),vDT]," / "))
	V = A.fop(v)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joPDE,joData):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

function *{ADDT,ARDT,vDT}(A::joPDEadjoint{ADDT,ARDT},v::joRHS{vDT})
	A.n == size(v,1) || throw(joPDEadjointException("shape mismatch"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(joPDE,joRHS):",A.name,typeof(A),vDT]," / "))
	V = A.fop(v)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joPDE,joData):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(joPDE,joProjection)
function *{ARDT,BDDT,CDT}(A::joPDE{CDT,ARDT},B::joProjection{BDDT,CDT})
	A.n == size(B,1) || throw(joPDEexception("shape mismatch"))
	return joModeling(A.info,A.model,B.geometry,A.geometry;options=A.options,DDT=CDT,RDT=ARDT)
end

function *{ARDT,BDDT,CDT}(A::joPDEadjoint{CDT,ARDT},B::joProjection{BDDT,CDT})
	A.n == size(B,1) || throw(joPDEadjointException("shape mismatch"))
	return joModeling(A.info,A.model,A.geometry,B.geometry,options=A.options,DDT=CDT,RDT=ARDT)'
end

# *(num,joPDE)
function *{ADDT,ARDT}(a::Number,A::joPDE{ADDT,ARDT})
	return joPDE{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1 -> jo_convert(ARDT,a*A.fop(v1),false),
		v2 -> jo_convert(ADDT,a*A.fop_T(v2),false)
		)
end

function *{ADDT,ARDT}(a::Number,A::joPDEadjoint{ADDT,ARDT})
	return joPDEadjoint{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1 -> jo_convert(ARDT,a*A.fop(v1),false),
		v2 -> jo_convert(ADDT,a*A.fop_T(v2),false)
		)
end

############################################################
## overloaded Basees +(...joPDE...), -(...joPDE...)

# +(joPDE,num)
function +{ADDT,ARDT}(A::joPDE{ADDT,ARDT},b::Number)
	return joPDE{ADDT,ARDT}("("*A.name*"+N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1 -> A.fop(v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
		v2 -> A.fop_T(v2)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
		)
end

function +{ADDT,ARDT}(A::joPDEadjoint{ADDT,ARDT},b::Number)
	return joPDEadjoint{ADDT,ARDT}("("*A.name*"+N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1->A.fop(v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
		v2->A.fop_T(v2)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
		)
end

# -(joPDE,num)
function -{ADDT,ARDT}(A::joPDE{ADDT,ARDT},b::Number)
	return joPDE{ADDT,ARDT}("("*A.name*"-N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1 -> A.fop(v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
		v2 -> A.fop_T(v2)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
		)
end

function -{ADDT,ARDT}(A::joPDEadjoint{ADDT,ARDT},b::Number)
	return joPDEadjoint{ADDT,ARDT}("("*A.name*"-N)",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1->A.fop(v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
		v2->A.fop_T(v2)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v2
		)
end

# -(joPDE)
-{DDT,RDT}(A::joPDE{DDT,RDT}) =
	joPDE{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1->-A.fop(v1),
		v2->-get(A.fop_T)(v2)
		)

-{DDT,RDT}(A::joPDEadjoint{DDT,RDT}) =
	joPDEadjoint{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.geometry,A.options,
		v1->-A.fop(v1),
		v2->-get(A.fop_T)(v2)
		)

############################################################
## Additional overloaded functions

# Subsample Modeling operator
function subsample{ADDT,ARDT}(F::joPDE{ADDT,ARDT}, srcnum)
	geometry = subsample(F.geometry,srcnum)		# Geometry of subsampled data container
	info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
	return joPDE(F.name, info, F.model, geometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

function subsample{ADDT,ARDT}(F::joPDEadjoint{ADDT,ARDT}, srcnum)
	geometry = subsample(F.geometry,srcnum)		# Geometry of subsampled data container
	info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
	return joPDEadjoint(F.name, info, F.model, geometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

getindex(F::joPDE,a) = subsample(F,a)
getindex(F::joPDEadjoint,a) = subsample(F,a)

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
# Reference: DOI: 10.1088/0143-0807/31/1/004
function getOverlapCore((R₁, α₁, d₁, ijk₁)::Tuple{Vector{Float64},
                                                  Float64, Float64, Vector{Int}},
                        (R₂, α₂, d₂, ijk₂)::Tuple{Vector{Float64},
                                                  Float64, Float64, Vector{Int}})
    α = α₁ + α₂
    res = d₁ * d₂ * (π/α)^1.5 * exp(-α₁ * α₂ / α * sum(abs2, R₁-R₂))
    for i in 1:3
        l₁, l₂ = ijk₁[i], ijk₂[i]
        res *= (-1.0)^(l₁) * factorial(l₁) * factorial(l₂) / α^(l₁+l₂)
        temp = 0
        for i₁ in 0:(l₁÷2), i₂ in 0:(l₂÷2)
            Ω = l₁ + l₂ - 2*(i₁ + i₂)
            for o in 0:(Ω÷2)
                temp += (-1)^o * factorial(Ω) *
                        α₁^(l₂ - i₁ - 2i₂ - o) *
                        α₂^(l₁ - i₂ - 2i₁ - o) *
                         α^(2 * (i₁ + i₂) + o) *
                        (R₁[i] - R₂[i])^(Ω-2o) /
                        ( 4^(i₁ + i₂ + o) *
                          factorial(i₁) *
                          factorial(i₂) *
                          factorial(o ) *
                          factorial(l₁ - 2i₁) *
                          factorial(l₂ - 2i₂) *
                          factorial(Ω  - 2o ) )
            end
        end
        res *= temp
    end
    res
end

@inline function getOverlap(bf1::FloatingGTBasisFuncs{<:Any, GN1, 1},
                            bf2::FloatingGTBasisFuncs{<:Any, GN2, 1}) where {GN1, GN2}
    R₁, R₂ = (bf1, bf2) .|> centerCoordOf
    ijk₁, ijk₂ = get.(Ref(ijkOrbitalList), getindex.(getfield.((bf1, bf2), :ijk), 1), false)
    res = 0
    for gf1 in bf1.gauss::NTuple{GN1, GaussFunc}, gf2 in bf2.gauss::NTuple{GN2, GaussFunc}
        α₁ = gf1.xpn()
        α₂ = gf2.xpn()
        d₁ = gf1.con()
        d₂ = gf2.con()
        bf1.normalizeGTO && (d₁ *= getNijkα(ijk₁..., α₁))
        bf1.normalizeGTO && (d₁ *= getNijkα(ijk₂..., α₂))
        res += getOverlapCore((R₁, α₁, d₁, ijk₁), (R₂, α₂, d₂, ijk₂))
    end
    res
end

getOverlap(b1::CompositeGTBasisFuncs{<:Any, 1}, b2::CompositeGTBasisFuncs{<:Any, 1}) =
[getOverlap(bf1, bf2) for bf1 in unpackBasisFuncs(b1), bf2 in unpackBasisFuncs(b2)] |> sum

getOverlap(b1::CompositeGTBasisFuncs, b2::CompositeGTBasisFuncs) =
[getOverlap(bf1, bf2) for bf1 in b1, bf2 in b2]

function getNijkα(i, j, k, α)
    (2α/π)^0.75 * ( (8α)^(i+j+k) * factorial(i) * factorial(j) * factorial(k) /
                (factorial(2i) * factorial(2j) * factorial(2k)) )^0.5
end


text:
const rewrite_rules = [
    :string_genus,
    :yang_baxter_star,
    :yang_baxter_triangle,
    :charge_rm_v,
    :charge_rm_f,
    :z_fusion,
    :x_fusion,
    :perm_rz,
    :identity,
    :genus_fusion,
    :swap_genus,
]

function rewrite!(tait::Tait, m::Match{:string_genus})
    rem_vertex!(tait, m.vertices[1])
    return tait
end

function rewrite!(tait::Tait, m::Match{:yang_baxter_star})
    hes = m.half_edges
    v0 = m.vertices[1]
    vs = [dst(tait, he) for he in hes]
    p1, p2, p3 = (quon_param(tait, he) for he in hes)
    q1, q2, q3 = yang_baxter_param_inv(p1, p2, p3)
    add_edge!(tait, vs[1], vs[2], face(tait, hes[1]), q1)
    add_edge!(tait, vs[1], vs[3], face(tait, hes[3]), q2)
    add_edge!(tait, vs[2], vs[3], face(tait, hes[2]), q3)
    rem_vertex!(tait, v0; update = true)
    return tait
end

function rewrite!(tait::Tait, m::Match{:yang_baxter_triangle})
    hes = m.half_edges
    vs = [src(tait, he) for he in hes]
    f = face(tait, hes[1])
    for ids in [(1, 2, 3), (2, 3, 1), (3, 1, 2)]
        p1, p2, p3 = (quon_param(tait, hes[ids[1]]), quon_param(tait, hes[ids[3]]), quon_param(tait, hes[ids[2]]))
        if !is_singular_yang_baxter(p1, p2, p3)
            q1, q2, q3 = yang_baxter_param(p1, p2, p3)
            new_v = add_vertex!(tait, f)
            add_edge!(tait, vs[ids[1]], new_v, face(tait, hes[ids[1]]), q1)
            add_edge!(tait, vs[ids[2]], new_v, face(tait, hes[ids[2]]), q2)
            add_edge!(tait, vs[ids[3]], new_v, face(tait, hes[ids[3]]), q3)
            for he in hes
                rem_edge!(tait, he; update = true)
            end
            break
        end
    end
    return tait
end

function rewrite!(tait::Tait{QuonParam{T}}, m::Match{:charge_rm_v}) where T
    hes = m.half_edges
    v = m.vertices[1]
    p0 = QuonParam{T}(zero(T), true)
    if (σ_inv(tait, hes[1]) == hes[end]) || (length(hes) == 1)
        for he in hes
            tait.quon_params[he] = p0
            tait.quon_params[twin(tait, he)] = p0
            rewrite!(tait, Match{:identity}([], [he]))
        end
        return tait
    end

    hes_identity = Int[]
    for i = 1:(length(hes)-1)
        he1 = hes[i]
        he2 = hes[i+1]
        v1 = dst(tait, he1)
        v2 = dst(tait, he2)
        f = face(tait, he1)
        push!(hes_identity, add_edge!(tait, v1, v2, f, p0)[1])
    end
    for i = 1:(length(hes))
        rem_edge!(tait, hes[i]; update = true)
        rewrite!(tait, Match{:identity}([], [hes_identity[i]]))
    end
    return tait
end

function rewrite!(tait::Tait{QuonParam{T}}, m::Match{:charge_rm_f}) where T
    hes = m.half_edges
    f = face(tait, hes[1])
    if next(tait, hes[end]) != hes[1]
        add_edge!(tait, src(tait, hes[1]), dst(tait, hes[end]), f, QuonParam{T}(Pi, false))
    end
    for he in hes
        rem_edge!(tait, he; update = true)
    end
    return tait
end

function rewrite!(tait::Tait{QuonParam{T}}, m::Match{:z_fusion}) where T
    he1, he2 = m.half_edges
    twin1 = twin(tait, he1)
    twin2 = twin(tait, he2)
    p1 = quon_param(tait, he1)
    p2 = quon_param(tait, he2)
    p = add_orthorgonal(p1, p2)
    p0 = QuonParam{T}(zero(p2.param), false)
    tait.quon_params[he1] = p
    tait.quon_params[twin1] = p
    tait.quon_params[he2] = p0
    tait.quon_params[twin2] = p0
    rewrite!(tait, Match{:identity}([], [he2]))
    return tait
end

function rewrite!(tait::Tait{QuonParam{T}}, m::Match{:x_fusion}) where T
    he1, he2 = m.half_edges
    twin1 = twin(tait, he1)
    twin2 = twin(tait, he2)
    p1 = quon_param(tait, he1)
    p2 = quon_param(tait, he2)
    p = add_parallel(p1, p2)
    p0 = QuonParam{T}(zero(p2.param), true)
    tait.quon_params[he1] = p
    tait.quon_params[twin1] = p
    tait.quon_params[he2] = p0
    tait.quon_params[twin2] = p0
    rewrite!(tait, Match{:identity}([], [he2]))
    return tait
end

function rewrite!(tait::Tait, m::Match{:perm_rz})
    v_ct, v_g = m.vertices
    f = face(tait, m.half_edges[2])
    old_he = m.half_edges[1]
    add_edge!(tait, v_ct, v_g, f, quon_param(tait, old_he))
    rem_edge!(tait, old_he; update = true)
    return tait
end

function rewrite!(tait::Tait, m::Match{:identity})
    he = m.half_edges[1]
    p = quon_param(tait, he)
    is_zero(p.param) || (p = change_direction(p))
    if is_parallel(p)
        contract_edge!(tait, he)
    else
        rem_edge!(tait, he; update = true)
    end
    return tait
end

function rewrite!(tait::Tait{QuonParam{T}}, m::Match{:genus_fusion}) where T
    he_g1, _ = m.half_edges
    g1, g2 = m.vertices
    f = face(tait, he_g1)
    new_he1, _ = add_edge!(tait, g1, g2, f, QuonParam{T}(zero(T), true))
    rewrite!(tait, Match{:identity}([], [new_he1]))
end

function rewrite!(tait::Tait, m::Match{:swap_genus})
    g = m.vertices[1]
    delete!(tait.genuses, g)
    return tait
end


text:
#InspectDR: Support for pan/zoom
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const PAN_STEPRATIO = 0.25 #Percentage of current extents
const ZOOM_STEPRATIO = 2.0 #How much to zoom in/out with mousewheel + keybindings

#==Drawing functions
===============================================================================#
function selectionbox_draw(ctx::CairoContext, selbb::BoundingBox, graphbb::BoundingBox, hallowed::Bool, vallowed::Bool)
	xmin = selbb.xmin; xmax = selbb.xmax
	ymin = selbb.ymin; ymax = selbb.ymax
	if !hallowed
		xmin = graphbb.xmin
		xmax = graphbb.xmax
	end
	if !vallowed
		ymin = graphbb.ymin
		ymax = graphbb.ymax
	end

	Cairo.save(ctx) #-----
	Cairo.set_source(ctx, COLOR_BLACK)
	setlinestyle(ctx, :dash, 1.0)
	Cairo.rectangle(ctx, BoundingBox(xmin, xmax, ymin, ymax))
	Cairo.stroke(ctx)
	Cairo.restore(ctx) #-----

	nothing
end

function selectionbox_draw(ctx::CairoContext, w::PlotWidget)
	if !w.sel.enabled; return; end
	selectionbox_draw(ctx, w.sel.bb, w.graphbb, w.hallowed, w.vallowed)
end

#==Switching to h/v lock (or releasing both)
===============================================================================#
function locdir_h(pwidget::PlotWidget)
	pwidget.hallowed = true
	pwidget.vallowed = false
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end
function locdir_v(pwidget::PlotWidget)
	pwidget.hallowed = false
	pwidget.vallowed = true
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end
function locdir_any(pwidget::PlotWidget)
	pwidget.hallowed = true
	pwidget.vallowed = true
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end

#==Scale/position widget control
===============================================================================#
function scalectrl_enabled(pwidget::PlotWidget)
	return Gtk.GAccessor.sensitive(pwidget.w_xscale)
end
function scalectrl_enabled(pwidget::PlotWidget, v::Bool)
	Gtk.GAccessor.sensitive(pwidget.w_xscale, v)
	Gtk.GAccessor.sensitive(pwidget.w_xpos, v)
end
#Apply current scale/position scrollbar values to plot extents:
function scalectrl_apply(pwidget::PlotWidget)
	xscale = getproperty(pwidget.xscale, :value, Int)
	xpos = getproperty(pwidget.xpos, :value, Float64)

	#Use transformed coordinate system:
	ext_full = rescale(getextents_full(pwidget.src), pwidget.src.axes)
	span = ext_full.xmax - ext_full.xmin
	center = (ext_full.xmax + ext_full.xmin) / 2
	vspan = span/xscale #Visible span
	xmin = center + span*xpos - vspan/2
	xmax = xmin + vspan

	#Update extents & redraw
	ext_new = PExtents2D(xmin, xmax, DNaN, DNaN)
	ext = merge(getextents_xfrm(pwidget.src), ext_new)
	setextents_xfrm(pwidget.src, ext)
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end

#==Basic zoom control
===============================================================================#
#Zoom to bounding box (in device coordinates):
function zoom(pwidget::PlotWidget, bb::BoundingBox)
	p1 = Point2D(bb.xmin, bb.ymin)
	p2 = Point2D(bb.xmax, bb.ymax)

	ext = getextents_xfrm(pwidget.src)
	xf = Transform2D(ext, pwidget.graphbb)
	p1 = ptmap_rev(xf, p1)
	p2 = ptmap_rev(xf, p2)

	setextents_xfrm(pwidget.src, PExtents2D(
		min(p1.x, p2.x), max(p1.x, p2.x),
		min(p1.y, p2.y), max(p1.y, p2.y)),
		pwidget.hallowed, pwidget.vallowed
	)
	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end

#Zoom in/out @ point (pt in plot coordinates)
function zoom(pwidget::PlotWidget, ext::PExtents2D, pt::Point2D, ratio::Float64)
	xspan = ext.xmax - ext.xmin
	yspan = ext.ymax - ext.ymin
	Δx = pt.x - ext.xmin
	Δy = pt.y - ext.ymin
	xmin = pt.x-ratio*Δx
	ymin = pt.y-ratio*Δy

	setextents_xfrm(pwidget.src, PExtents2D(
		xmin, xmin + ratio*xspan,
		ymin, ymin + ratio*yspan
	))
	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end

#Zoom in/out, centered on current extents
function zoom(pwidget::PlotWidget, ratio::Float64)
	ext = getextents_xfrm(pwidget.src)
	pt = Point2D((ext.xmin+ext.xmax)/2, (ext.ymin+ext.ymax)/2)
	zoom(pwidget, ext, pt, ratio)
end
zoom_out(pwidget::PlotWidget, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, stepratio)
zoom_in(pwidget::PlotWidget, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, 1/stepratio)

#Zoom in/out around specified device coordinates:
function zoom(pwidget::PlotWidget, x::Float64, y::Float64, ratio::Float64)
	pt = Point2D(x, y)
	ext = getextents_xfrm(pwidget.src)
	xf = Transform2D(ext, pwidget.graphbb)
	pt = ptmap_rev(xf, pt)
	zoom(pwidget, ext, pt, ratio)
end
zoom_out(pwidget::PlotWidget, x::Float64, y::Float64, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, x, y, stepratio)
zoom_in(pwidget::PlotWidget, x::Float64, y::Float64, stepratio::Float64=ZOOM_STEPRATIO) =
	zoom(pwidget, x, y, 1/stepratio)

function zoom_xfull(pwidget::PlotWidget)
	#TODO
end
function zoom_full(pwidget::PlotWidget, hallowed::Bool=true, vallowed::Bool=true)
	setextents(pwidget.src, PExtents2D(), hallowed, vallowed) #Reset desired extents
	if hallowed
		scalectrl_enabled(pwidget, false) #Suppress updates from setproperty!
		setproperty!(pwidget.xscale, :value, Int(1))
		setproperty!(pwidget.xpos, :value, Float64(0))
		scalectrl_enabled(pwidget, true)
		scalectrl_apply(pwidget)
	else
		render(pwidget)
		Gtk.draw(pwidget.canvas)
	end
end
zoom_hfull(pwidget::PlotWidget) = zoom_full(pwidget, true, false)
zoom_vfull(pwidget::PlotWidget) = zoom_full(pwidget, false, true)

#==Box-zoom control
===============================================================================#
function boxzoom_setstart(pwidget::PlotWidget, x::Float64, y::Float64)
	locdir_any(pwidget)
	pwidget.sel.enabled = true
	pwidget.sel.bb = BoundingBox(x, x, y, y)
end
function boxzoom_cancel(pwidget::PlotWidget)
	pwidget.sel.enabled = false
	Gtk.draw(pwidget.canvas)
end
function boxzoom_complete(pwidget::PlotWidget, x::Float64, y::Float64)
	pwidget.sel.enabled = false
	bb = pwidget.sel.bb
	pwidget.sel.bb = BoundingBox(bb.xmin, x, bb.ymin, y)
	zoom(pwidget, pwidget.sel.bb)
end
#Set end point of boxzoom area:
function boxzoom_setend(pwidget::PlotWidget, x::Float64, y::Float64)
	bb = pwidget.sel.bb
	pwidget.sel.bb = BoundingBox(bb.xmin, x, bb.ymin, y)
	Gtk.draw(pwidget.canvas)
end

#==Basic pan control
===============================================================================#
function pan_xratio(pwidget::PlotWidget, panstepratio::Float64)
	ext = getextents_xfrm(pwidget.src)
	panstep = panstepratio*(ext.xmax-ext.xmin)
	setextents_xfrm(pwidget.src, PExtents2D(
		ext.xmin+panstep, ext.xmax+panstep,
		ext.ymin, ext.ymax
	))
	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end
function pan_yratio(pwidget::PlotWidget, panstepratio::Float64)
	ext = getextents_xfrm(pwidget.src)
	panstep = panstepratio*(ext.ymax-ext.ymin)
	setextents_xfrm(pwidget.src, PExtents2D(
		ext.xmin, ext.xmax,
		ext.ymin+panstep, ext.ymax+panstep
	))
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end

pan_left(pwidget::PlotWidget) = pan_xratio(pwidget, -PAN_STEPRATIO)
pan_right(pwidget::PlotWidget) = pan_xratio(pwidget, PAN_STEPRATIO)
pan_up(pwidget::PlotWidget) = pan_yratio(pwidget, PAN_STEPRATIO)
pan_down(pwidget::PlotWidget) = pan_yratio(pwidget, -PAN_STEPRATIO)

#==Mouse-pan control
===============================================================================#
#Δy/Δy: in device coordinates
function mousepan_delta(pwidget::PlotWidget, ext::PExtents2D, Δx::Float64, Δy::Float64)
	#Convert to plot coordinates:
	xf = Transform2D(ext, pwidget.graphbb)
	Δvec = vecmap_rev(xf, Point2D(-Δx, -Δy))

	setextents_xfrm(pwidget.src, ext) #Restore original extents before overwriting
	setextents_xfrm(pwidget.src, PExtents2D(
		ext.xmin+Δvec.x, ext.xmax+Δvec.x,
		ext.ymin+Δvec.y, ext.ymax+Δvec.y),
		pwidget.hallowed, pwidget.vallowed
	)
	scalectrl_enabled(pwidget, false) #Scroll bar control no longer valid
	render(pwidget)
	Gtk.draw(pwidget.canvas)
end
function mousepan_setstart(pwidget::PlotWidget, x::Float64, y::Float64)
	locdir_any(pwidget)
	pwidget.sel.bb = BoundingBox(x, x, y, y) #Tracks start/end pos
	pwidget.sel.ext_start = getextents_xfrm(pwidget.src)
end
function mousepan_cancel(pwidget::PlotWidget)
	mousepan_delta(pwidget, pwidget.sel.ext_start, 0.0, 0.0)
end
function mousepan_complete(pwidget::PlotWidget, x::Float64, y::Float64)
	#Already panned.
end
#Set new point of mousepan operation:
function mousepan_move(pwidget::PlotWidget, x::Float64, y::Float64)
	bb = pwidget.sel.bb
	bb = BoundingBox(bb.xmin, x, bb.ymin, y)
	pwidget.sel.bb = bb
	Δx = bb.xmax-bb.xmin; Δy = bb.ymax-bb.ymin
	mousepan_delta(pwidget, pwidget.sel.ext_start, Δx, Δy)
end

#Last line


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

module SimdPerf
import Perftests: @perf, meta

# Compute y += a*x using @simd for vectors x and y
function simd_axpy( a, x, y )
    # LLVM's auto-vectorizer typically vectorizes this loop even without @simd
    @simd for i=1:length(x)
        @inbounds y[i] += a*x[i]
    end
end

# Compute a sum reduction using @simd for vector x, with certain starting and ending indices
function sum_reduce(x, istart, iend)
    s = zero(eltype(x))
    @simd for i = istart:iend
        @inbounds s += x[i]
    end
    s
end
function flog_sum_reduce( m, x )
    s = zero(eltype(x))
    for j=1:m
        # Try different starting and ending indices.
        s += sum_reduce(x,j,length(x)-(j-1))
    end
    return s
end

# Inner product of x and y
function simd_dot( x, y )
    s = zero(eltype(x))
    @simd for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end

# This guy, we will keep in his own file since he's nontrivial
include("seismic_fdtd.jl")

for t in [Float32,Float64]
    n = 1000
    x = rand(t,n)
    y = rand(t,n)
    a = convert(t,0.5)

    # Test axpy performance
    @perf simd_axpy(a,x,y) meta("axpy", string(t), "SIMD BLAS axpy for type $t")

    # Test sum reduction performance
    @perf flog_sum_reduce(1000,x) meta("sum", string(t), "SIMD sum reduction over array of type $t")

    # Test dot performance
    @perf simd_dot(x, y) meta("dot", string(t), "SIMD inner product for type $t")

    # Test 2D seismic simulation
    m = 200
    n = 200
    A = fill(convert(t,0.2),m,n)
    B = fill(convert(t,0.25),m,n)
    U = rand(t,m,n) .- convert(t,.5)
    Vx = zeros(t,m,n)
    Vy = zeros(t,m,n)
    @perf flog_fdtd(10,U,Vx,Vy,A,B) meta("seismic_fdtd", string(t), "2D finite-difference seismic simulation for $t")
end

end # module


text:
using RDatasets, Distributions, Random, LinearAlgebra
Random.seed!(0)

K = 3
df = dataset("cluster", "xclara")
n,_ = size(df)
dataPoints = [convert(Array{Float64,1},df[i,:]) for i in 1:n]
shuffle!(dataPoints)

xMin,xMax = minimum(first.(dataPoints)),maximum(first.(dataPoints))
yMin,yMax = minimum(last.(dataPoints)),maximum(last.(dataPoints))

means = [[rand(Uniform(xMin,xMax)),rand(Uniform(yMin,yMax))]  for _ in 1:K]
labels = rand(1:K,n)
prevMeans = -means

while norm(prevMeans - means) > 0.001
    prevMeans = means
    labels = [argmin([norm(means[i]-x) for i in 1:K]) for x in dataPoints]
    means = [sum(dataPoints[labels .== i])/sum(labels .==i) for i in 1:K]
end

countResult = [sum(labels .== i) for i in 1:K]
println("Counts of clusters (manual implementation): ", countResult)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "PrGP" begin

    v1 = [1, 0, 0]
    v2 = [0, 1, 0]
    t  = [3, 3, 4]
    @test optimal_coeff(v1, v2, t) ≈ [3, 3]

    t  = [3, 3, 100101010]
    @test optimal_coeff(v1, v2, t) ≈ [3, 3]

    v1 = [1, 0, 1]
    v2 = [0, 1, 0]
    t  = [3, 3, 4]
    @test optimal_coeff(v1, v2, t) ≈ [3.5, 3.0]

    v1 = [1, 1, 1]
    v2 = [-1, 1, 1]
    @test optimal_coeff(v1, v2, t) ≈ [3.25, 0.25]

    v1 = v2 = [1,1,1]
    t = [10, 10, 10,]
    @test optimal_coeff(v1, v2, t) ≈ [5, 5]

    cfg = Config(population_size=100, num_generations=10, p_crossover=0.5, p_mutation=0.5, mutation_step=0.01, min_depth=1, max_depth=6, tournament_size=4, pset=PrimitiveSet(), create_rand_tree=generate_grow, create_population=ramped_half_and_half, fitness=rmse)
    dataset = Dataset(zeros(1, 1), [3, 3, 4])

    e1 = Entry(1, 1, JuliaGSGP.none, JuliaGSGP.dummy_node(), (-1, -1, -1), (-1, -1, -1), [1, 0, 1], [1, 0, 1], 10, 10)
    e2 = Entry(1, 1, JuliaGSGP.none, JuliaGSGP.dummy_node(), (-1, -1, -1), (-1, -1, -1), [0, 1, 0], [0, 1, 0], 10, 10)
    e3 = projection_crossover_entry(e1, e2, dataset, dataset, cfg)

    @test e3.train_semantics ≈ [3.5, 3., 3.5]

end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
"""
This set of functions queries if a given coordinate/grid is within a region defined either
by a set of region [N,S,E,W] bounds or by longitude/latitude vectors.

Other additional features include:
    * Transformation of longitude coordinates from [-180,180] to [0,360] and vice versa

"""

## Is Point in GeoRegion?

"""
    ispointinregion(plon, plat, rlon, rlat; throw=true) -> Bool

Check if a point with coordinates (`plon`,`plat`) is found within a region defined by the longitude and latitude vectors `rlon` and `rlat`

Arguments:
* `plon::Real` : Longitude of the point in question.
* `plat::Real` : Latitude of the point in question.
* `rlon::Vector{<:Real}` : Longitude vector spanning the region.  Points should be evenly spaced.
* `rlat::Vector{<:Real}` : Latitude vector spanning the region.  Points should be evenly spaced.

Keyword Arguments:
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function ispointinregion(
    plon::Real, plat::Real,
    rlon::Vector{<:Real}, rlat::Vector{<:Real};
    throw::Bool=true
)

    rW = mod(minimum(rlon),360); rS = minimum(rlat);
    rE = mod(maximum(rlon),360); rN = maximum(rlat);
    tlon = abs((rlon[2]-rlon[1])*3/4);
    tlat = abs((rlat[2]-rlat[1])*3/4);

    return checkpoint(plon,plat,rN,rS,rE,rW,tlon,tlat,throw)

end

"""
    ispointinregion(plon, plat, regionbounds; tlon=0, tlat=0, throw=true) -> Bool

Check if a point with coordinates (`plon`,`plat`) is found within a region defined by the bounds `regionbounds` in [`N`,`S`,`E`,`W`] format.

Arguments:
* `plon::Real` : Longitude of the point in question.
* `plat::Real` : Latitude of the point in question.
* `regionbounds::Vector{<:Real}` : Vector defining the [North,South,East,West] bounds of the region.

Keyword Arguments:
* `tlon::Real` : Threshold for longitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the longitude bounds
* `tlat::Real` : Threshold for latitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the latitude bounds
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function ispointinregion(
    plon::Real, plat::Real,
    regionbounds::Vector{<:Real};
    tlon::Real=0, tlat::Real=0,
    throw::Bool=true
)

    rN,rS,rE,rW = regionbounds; rW = mod(rW,360); rE = mod(rE,360);

    return checkpoint(plon,plat,rN,rS,rE,rW,tlon,tlat,throw)

end

"""
    ispointinregion(pcoord, rlon, rlat; throw=true) -> Bool

Check if a point with coordinates `pcoord = [plon,plat]` is found within a region defined by the longitude and latitude vectors `rlon` and `rlat`.

Arguments:
* `pcoord::Vector{<:Real}` : Vector defining coordinates of point in `[plon,plat]`
* `rlon::Vector{<:Real}` : Longitude vector spanning the region.  Points should be evenly spaced.
* `rlat::Vector{<:Real}` : Latitude vector spanning the region.  Points should be evenly spaced.

Keyword Arguments:
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function ispointinregion(
    pcoord::Vector{<:Real},
    rlon::Vector{<:Real}, rlat::Vector{<:Real};
    throw::Bool=true
)

    plon,plat = pcoord;
    rW = mod(minimum(rlon),360); rS = minimum(rlat);
    rE = mod(maximum(rlon),360); rN = maximum(rlat);
    tlon = abs((rlon[2]-rlon[1])*3/4);
    tlat = abs((rlat[2]-rlat[1])*3/4);

    return checkpoint(plon,plat,rN,rS,rE,rW,tlon,tlat,throw)

end

"""
    ispointinregion(pcoord, regionbounds; tlon=0, tlat=0, throw=true) -> Bool

Check if a point with coordinates `pcoord = [plon,plat]` is found within a region defined by the bounds `regionbounds` in [`N`,`S`,`E`,`W`] format.

Arguments:
* `pcoord::Vector{<:Real}` : Vector defining coordinates of point in `[plon,plat]`
* `regionbounds::Vector{<:Real}` : Vector defining the [North,South,East,West] bounds of the region.

Keyword Arguments:
* `tlon::Real` : Threshold for longitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the longitude bounds
* `tlat::Real` : Threshold for latitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the latitude bounds
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function ispointinregion(
    pcoord::Vector{<:Real},
    regionbounds::Vector{<:Real};
    tlon::Real=0, tlat::Real=0,
    throw::Bool=true
)

    plon,plat = pcoord; rN,rS,rE,rW = regionbounds;
    rW = mod(rW,360); rE = mod(rE,360);

    return checkpoint(plon,plat,rN,rS,rE,rW,tlon,tlat,throw)

end

"""
    checkpoint(plon, plat, rN, rS, rE, rW, tlon, tlat, throw) -> Bool

Check if a point with longitude and latitude coordinates `plon` and `plat` is found within a region defined by the bounds [`rN`,`rS`,`rE`,`rW`] format, with `tlon` and `tlat` being the threshold/margin of error for the longitude and latitude respectively.

Arguments:
* `plon::Real` : Longitude of the point in question.
* `plat::Real` : Latitude of the point in question.
* `rN::Real` : North bound of region
* `rS::Real` : South bound of region
* `rE::Real` : East bound of region
* `rW::Real` : West bound of region
* `tlon::Real` : Threshold for longitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the longitude bounds
* `tlat::Real` : Threshold for latitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the latitude bounds
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function checkpoint(
    plon::Real, plat::Real,
    rN::Real, rS::Real, rE::Real, rW::Real,
    tlon::Real, tlat::Real,
    throw::Bool
)

    if (rW < rE && (plon < (rW-tlon) || plon > (rE+tlon))) ||
        (rW > rE && (plon < (rW-tlon) && plon > (rE+tlon))) ||
        plat < (rS-tlat) || plat > (rN+tlat)

        if throw
             error("Requested coordinates [$(plon),$(plat)] are not within the specified region boundaries.")
        else return false
        end

    else; return true

    end

end

## Is Grid in GeoRegion?

"""
    isgridinregion(gridbounds, regionbounds; tlon=0, tlat=0, throw=true) -> Bool

Check if a grid defined by `gridbounds = [gN,gS,gE,gW]` is found within a predefined `GeoRegion`.

Arguments:
* `gridbounds::Vector{<:Real}` : Vector defining the [North,South,East,West] bounds of the grid in question.
* `regionbounds::Vector{<:Real}` : Vector defining the [North,South,East,West] bounds of the region.

Keyword Arguments:
* `tlon::Real` : Threshold for longitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the longitude bounds
* `tlat::Real` : Threshold for latitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the latitude bounds
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function isgridinregion(
    gridbounds::Vector{<:Real},
    regionbounds::Vector{<:Real};
    tlon::Real=0, tlat::Real=0,
    throw::Bool=true
)

    N,S,E,W = gridbounds; rN,rS,rE,rW = regionbounds;
    E = mod(E,360); rE = mod(rE,360); W = mod(W,360); rW = mod(rW,360);

    return checkgrid(N,S,E,W,rN,rS,rE,rW,tlon,tlat,throw)

end

"""
    isgridinregion(gridbounds, rlon, rlat; throw=true) -> Bool

Check if a grid defined by `gridbounds = [gN,gS,gE,gW]` is found within a predefined `GeoRegion`.

Arguments:
* `gridbounds::Vector{<:Real}` : Vector defining the [North,South,East,West] bounds of the grid in question.
* `rlon::Vector{<:Real}` : Longitude vector spanning the region.  Points should be evenly spaced.
* `rlat::Vector{<:Real}` : Latitude vector spanning the region.  Points should be evenly spaced.

Keyword Arguments:
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function isgridinregion(
    gridbounds::Vector{<:Real},
    rlon::Vector{<:Real},rlat::Vector{<:Real};
    throw::Bool=true
)

    N,S,E,W = gridbounds; E = mod(E,360); W = mod(W,360);
    rE = mod(maximum(rlon),360); rN = maximum(rlat);
    rW = mod(minimum(rlon),360); rS = minimum(rlat);
    tlon = abs((rlon[2]-rlon[1])*3/4);
    tlat = abs((rlat[2]-rlat[1])*3/4);

    return checkgrid(N,S,E,W,rN,rS,rE,rW,tlon,tlat,throw)

end

"""
    checkgrid(gN, gS, gE, gW, rN, rS, rE, rW, tlon, tlat, throw) -> Bool

Check if a grid with bounds [`gN`,`gS`,`gE`,`gW`] is found within a region defined by the bounds [`rN`,`rS`,`rE`,`rW`] format, with `tlon` and `tlat` being the threshold/margin of error for the longitude and latitude respectively.

Arguments:
* `gN::Real` : North bound of grid
* `gS::Real` : South bound of grid
* `gE::Real` : East bound of grid
* `gW::Real` : West bound of grid
* `rN::Real` : North bound of region
* `rS::Real` : South bound of region
* `rE::Real` : East bound of region
* `rW::Real` : West bound of region
* `tlon::Real` : Threshold for longitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the longitude bounds
* `tlat::Real` : Threshold for latitude bounds in °.  If set to 0, means that there is no leniency for the point to fall outside the latitude bounds
* `throw::Bool` : If `throw=true`, then if (`plon`,`plat`) is not within the region, an error is thrown and the program stops running.
"""
function checkgrid(
    gN::Real, gS::Real, gE::Real, gW::Real,
    rN::Real, rS::Real, rE::Real, rW::Real,
    tlon::Real, tlat::Real,
    throw::Bool
)

    if (rW < rE && (gE < (rW-tlon) || gE > (rE+tlon))) ||
        (rW > rE && (gE < (rW-tlon) && gE > (rE+tlon))) ||
        (rW < rE && (gW < (rW-tlon) || gW > (rE+tlon))) ||
        (rW > rE && (gW < (rW-tlon) && gW > (rE+tlon))) ||
        gN < (rS-tlat) || gN > (rN+tlat) ||
        gS < (rS-tlat) || gS > (rN+tlat)

        if throw
             error("Requested grid points [$gN,$gS,$gE,$gW] are not within the specified region boundaries.")
        else return false
        end

    else; return true

    end

end

## Tranformation of Coordinates

from180to0360(lon::Real) = lon + 360
from0360to180(lon::Real) = lon - 360


text:
# Autogenerated wrapper script for lazygit_jll for aarch64-linux-musl
export lazygit

JLLWrappers.@generate_wrapper_header("lazygit")
JLLWrappers.@declare_executable_product(lazygit)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        lazygit,
        "bin/lazygit",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function main_validate(args)
    ln2("="^70)
    ln2("main_validate: BEGIN")

    # Résolution de l'action
    ln1("main_validate: fichier d'instance : $(args[:infile])")
    inst = Instance(args[:infile])

    ln1("main_validate: fichier de la solution : $(args[:solfile])")
    sol = Solution(inst, args[:solfile])

    ln1("main_validate: examen de la solution (par get_viol_description)")
    nbviols, violtxt = get_viol_description(sol)

    if nbviols == 0
        msg = "Solution correcte de coût : $(sol.cost)"
        println(to_sc(msg, :GREEN))
        ln3(to_s(sol))
    else
        msg = "Solution incorrecte : il y a $nbviols erreurs !"
        println(to_sc(msg, :RED))
        ln1(violtxt)
        ln3(to_s(sol))
    end

    ln1("main_validate: END")
    ln2("="^70)
end

main_validate(Args.args)


text:
module BitConverter

export bytes, to_big, to_int

"""
    bytes(x::Integer; len::Integer, little_endian::Bool)
    -> Vector{len, UInt8}

Convert an Integer `x` to a Vector{UInt8}
Options (not available for `x::BigInt`):
- `len` to define a minimum Vector lenght in bytes, result will show no leading
zero by default.
- set `little_endian` to `true` for a result in little endian byte order, result
in big endian order by default.

    julia> bytes(32974)
    2-element Array{UInt8,1}:
     0x80
     0xce

    julia> bytes(32974, len=4)
    4-element Array{UInt8,1}:
     0x00
     0x00
     0x80
     0xce

    julia> bytes(32974, little_endian=true)
    2-element Array{UInt8,1}:
     0xce
     0x80

"""
function bytes(x::Integer; len::Integer=0, little_endian::Bool=false)
    result = reinterpret(UInt8, [hton(x)])
    i = findfirst(x -> x != 0x00, result)
    if len != 0
        i = length(result) - len + 1
    end
    result = result[i:end]
    if little_endian
        reverse!(result)
    end
    return result
end

function bytes(x::BigInt)
    n_bytes_with_zeros = x.size * sizeof(Sys.WORD_SIZE)
    uint8_ptr = convert(Ptr{UInt8}, x.d)
    n_bytes_without_zeros = 1

    if ENDIAN_BOM == 0x04030201
        # the minimum should be 1, else the result array will be of
        # length 0
        for i in n_bytes_with_zeros:-1:1
            if unsafe_load(uint8_ptr, i) != 0x00
                n_bytes_without_zeros = i
                break
            end
        end

        result = Array{UInt8}(undef, n_bytes_without_zeros)

        for i in 1:n_bytes_without_zeros
            @inbounds result[n_bytes_without_zeros + 1 - i] = unsafe_load(uint8_ptr, i)
        end
    else
        for i in 1:n_bytes_with_zeros
            if unsafe_load(uint8_ptr, i) != 0x00
                n_bytes_without_zeros = i
                break
            end
        end

        result = Array{UInt8}(undef, n_bytes_without_zeros)

        for i in 1:n_bytes_without_zeros
            @inbounds result[i] = unsafe_load(uint8_ptr, i)
        end
    end
    return result
end

"""
    Int(x::Vector{UInt8}; little_endian::Bool)
    -> Integer

Convert a Vector{UInt8} an Integer, a BigInt above 64 bits
Optionally set `little_endian` to `true` if the input as such byte order,
input is considered big endian by default.

    julia> Int([0x01, 0x00])
    256

    julia> Int([0x01, 0x00], little_endian=true)
    1

"""
function to_int(x::Vector{UInt8}; little_endian::Bool=false)
    if length(x) > 8
        to_big(x)
    else
        missing_zeros = div(Sys.WORD_SIZE, 8) -  length(x)
        if missing_zeros > 0
            if little_endian
                for i in 1:missing_zeros
                    push!(x,0x00)
                end
            else
                for i in 1:missing_zeros
                    pushfirst!(x,0x00)
                end
            end
        end
        if ENDIAN_BOM == 0x04030201 && little_endian
        elseif ENDIAN_BOM == 0x04030201 || little_endian
            reverse!(x)
        end
        return reinterpret(Int, x)[1]
    end
end

"""
    to_big(x::Vector{UInt8}) -> BigInt

Convert a Vector{UInt8} of any lenght to a BigInt.
Considers the input a big endian.

    julia> to_big([0x01, 0x00])
    256

"""
function to_big(x::Vector{UInt8})
    hex = bytes2hex(x)
    return parse(BigInt, hex, base=16)
end

end # module

@deprecate Int(x::Vector{UInt8}; little_endian::Bool) to_int(x; little_endian=little_endian)
@deprecate big(x::Vector{UInt8}) to_big(x)


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
root = 0
isroot = rank == root

for T in Base.uniontypes(MPI.MPIDatatype)

    A = isroot ? ArrayType{T}(1:size) : ArrayType{T}(undef, 1)

    # Non Allocating version
    B = ArrayType{T}(undef, 1)
    MPI.Scatter!(A, B, root, comm)
    @test Array(B)[1] == T(rank+1)

    # In place version
    B = isroot ? copy(A) : ArrayType{T}(undef, 1)
    if root == MPI.Comm_rank(comm)
        MPI.Scatter!(UBuffer(B,1), MPI.IN_PLACE, root, comm)
    else
        MPI.Scatter!(nothing, B, root, comm)
    end
    @test Array(B)[1] == T(rank+1)

    # Test throwing
    if isroot
        B = ArrayType{T}(undef, 0)
        @test_throws DivideError MPI.Scatter!(A, B, root, comm)
        B = ArrayType{T}(undef, 8)
        @test_throws AssertionError MPI.Scatter!(A, B, root, comm)
    end
end

MPI.Finalize()
@test MPI.Finalized()


text:
file = "inputs/models/plant_coffee.yml"

@testset "read_model()" begin
    model = read_model(file)

    @test all([haskey(model,i) for i in ["Metamer","Leaf"]])

    model = model["Leaf"]
    # @test typeof(model)
    @test typeof(model) <: LeafModels
    @test typeof(model.stomatal_conductance) == Medlyn{Float64}
    @test typeof(model.interception) == Translucent{Float64}
    @test typeof(model.photosynthesis) == Fvcb{Float64}

    @test model.stomatal_conductance.g0 == -0.03
    @test model.stomatal_conductance.g1 == 12.0

    @test model.interception.transparency == 0.0
    @test model.interception.optical_properties == σ(0.15, 0.9)

    @test model.photosynthesis.VcMaxRef == 200.0 # Given in the file
    @test model.photosynthesis.O₂ == 210.0 # Use default value
end;


text:
using Base.Test
using AutoEnvs

function test_basics()
    # ctor
    filepath = Pkg.dir("NGSIM", "data", "trajdata_i80_trajectories-0400-0415.txt")
    n_veh = 2
    params = Dict(
        "trajectory_filepaths"=>[filepath],
        "H"=>50,
        "n_veh"=>n_veh
    )
    env = MultiagentNGSIMEnv(params)

    # reset, step
    x = reset(env)
    @test size(x) == (2, 66)
    a = zeros(n_veh, 2)
    nx, r, terminal, infos = step(env, a)
    @test x != nx
    @test terminal == [false, false]
    nnx, r, terminal, infos = step(env, a)
    @test nx != nnx

    # obs spec
    shape, spacetype, infos = observation_space_spec(env)
    @test spacetype == "Box"
    @test in("high", keys(infos))
    @test in("low", keys(infos))

    # does accel reflect applied value
    fns = obs_names(env)
    acc_idx = [i for (i,n) in enumerate(fns) if "accel" == n][1]
    tur_idx = [i for (i,n) in enumerate(fns) if "turn_rate_global" == n][1]
    a = zeros(n_veh, 2)
    a[:, 2] = 1
    for _ in 1:10
        nx, _, _, _ = step(env, a)
    end

    @test all(abs.(nx[:, acc_idx]) .<= 1e-1)
    @test all(abs.(nx[:, tur_idx] - 1) .<= 1e-1)

    # test infos
    reset(env)
    _, _, _, infos = step(env, ones(n_veh, 2))
    @test all(infos["rmse_pos"] .!= 0.)

    # test multiple episodes
    reset(env, [true for _ in 1:n_veh])
    for _ in 1:50
        nx, _, terminal, _ = step(env, a)
    end
    @test all(terminal)
    _, _, terminal, _ = step(env, a)
    @test all(.!terminal)
    for _ in 1:49
        _, _, terminal, _ = step(env, a)
    end
    @test all(terminal)

end

function test_render()
    srand(2)
    n_veh = 50
    filepath = Pkg.dir("NGSIM", "data", "trajdata_i101_trajectories-0805am-0820am.txt")
    params = Dict(
        "trajectory_filepaths"=>[filepath],
        "H"=>50,
        "primesteps"=>50,
        "n_veh"=>n_veh,
        "remove_ngsim_veh"=>true
    )
    env = MultiagentNGSIMEnv(params)

    reset(env)
    a = ones(n_veh, 2)
    for _ in 1:50
        render(env)
        x, r, terminal, _ = step(env, a)
    end
end

@time test_basics()
# manual tests
# @time test_render()


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Random
using ElasticArrays, ArraysOfArrays

@testset "array_shape" begin
    @test_throws ArgumentError @inferred(ValueShapes._valshapeoftype(Vector{Int}))

    @test @inferred(valshape(rand(3))) == ArrayShape{Float64,1}((3,))
    @test @inferred(valshape(rand(3, 4, 5))) == ArrayShape{Float64,3}((3, 4, 5))

    @inferred(ValueShapes.nonabstract_eltype(ArrayShape{Complex,3}((3, 4, 5)))) == Complex{Float64}

    @test @inferred(totalndof(ArrayShape{Float64,1}((3,)))) == 3
    @test @inferred(totalndof(ArrayShape{Complex,3}((3, 4, 5)))) == 120

    @test valshape(@inferred(Array(undef, ArrayShape{Complex,3}((3, 4, 5))))) == ArrayShape{Complex{Float64},3}((3, 4, 5))
    @test typeof(@inferred(Array(undef, ArrayShape{Complex,3}((3, 4, 5))))) == Array{Complex{Float64},3}
    @test size(@inferred(Array(undef, ArrayShape{Complex,3}((3, 4, 5))))) == (3, 4, 5)
    @test typeof(@inferred(Array{Float32}(undef, ArrayShape{Real,1}((3,))))) == Array{Float32,1}

    @test valshape(@inferred(ElasticArray(undef, ArrayShape{Complex,3}((3, 4, 5))))) == ArrayShape{Complex{Float64},3}((3, 4, 5))
    @test typeof(@inferred(ElasticArray(undef, ArrayShape{Complex,3}((3, 4, 5))))) == ElasticArray{Complex{Float64},3,2}
    @test size(@inferred(ElasticArray(undef, ArrayShape{Complex,3}((3, 4, 5))))) == (3, 4, 5)
    @test typeof(@inferred(ElasticArray{Float32}(undef, ArrayShape{Real,2}((3,4))))) == ElasticArray{Float32,2,1}

    shape = ArrayShape{Real}(2,3)
    A = @inferred(shape(undef))
    @test typeof(A) == Array{Float64,2}
    @test size(A) == (2, 3)
    @test @inferred(valshape(A)) == ArrayShape{Float64}(2,3)
    @test @inferred(shape([1, 2, 3, 4, 5, 6])) == [1 3 5; 2 4 6]
    @test_throws ArgumentError shape([1, 2, 3, 4, 5, 6, 7])

    @test eltype(@inferred(Vector{Float32}(undef, shape))) == Float32
    @test eltype(eltype(@inferred(VectorOfSimilarVectors{Float32}(shape)))) == Float32

    @test valshape(shape.(push!(@inferred(VectorOfSimilarVectors(shape)), @inferred(Vector(undef, shape))))[1]) == valshape(shape(undef))
end


text:
############################################################################
#
#  Airframe.jl
#
#  Purpose:    Single file that contains all functions needed to simulate a
#              6-DOF Airframe.
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################

# Initialize Position
function InitPos( time, latitude, longitude, altitude )

    ## Transformation Matrices ECEF, ECI, NED
    Tr_ei = TR_EI(time);
    Tr_ie = Tr_ei';
    Tr_ne = TR_NE(latitude, longitude);

    ## Position Vectors
    rm_ecef = LLH2ECEF(latitude, longitude, altitude);
    rm_eci = Tr_ie*rm_ecef;

    println(altitude);
    ## Populating Position Structure
    position = WGS_Pos(time, latitude, longitude, altitude, rm_ecef, rm_eci, Tr_ei, Tr_ne);
    return position;
end

############################################################################
# Initialize Velocity
function InitVel( vmag::Float64, gamma::Float64, azimuth::Float64, position::WGS_Pos )

    time = position.time;
    Tr_ei = position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = position.Tr_ne;
    Tr_en = Tr_ne';

    wb_ei_eci = [0.0;0.0;WEII3];

    vm_ned = [
        vmag*cos(gamma)*cos(azimuth)
        vmag*cos(gamma)*sin(azimuth)
        -vmag*sin(gamma)
    ];

    vm_ecef = Tr_en*vm_ned;
    vm_eci = Tr_ie*vm_ecef + cross(wb_ei_eci, position.rm_eci);

    velocity = Velocity_States( time, vmag, gamma, azimuth, vm_ned, vm_ecef, vm_eci );
    return velocity;
end

############################################################################
# Initialize Euler
function InitEuler( Roll::Float64, Pitch::Float64, Heading::Float64 )
    Tr_bn = TR_BN( Roll, Pitch, Heading);

    euler = Euler_Angles(Roll, Pitch, Heading, Tr_bn);

    return euler;
end

############################################################################
# Initialize Rates
function InitAngleRates(wb_bn_body::Array{Float64,1}, position::WGS_Pos, euler::Euler_Angles)

    ## Transformation Matrices
    Tr_ei = position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = position.Tr_ne;
    Tr_en = Tr_ne';
    Tr_bn = euler.Tr_bn;
    Tr_nb = Tr_bn';
    Tr_bi = Tr_bn*Tr_ne*Tr_ei;
    Tr_ib = Tr_bi';

    ## Earth Rotation Rate
    wb_ne_body = [0.0;0.0;0.0];
    wb_ne_eci = Tr_ib*wb_ne_body;
    wb_ei_eci = [0.0;0.0;WEII3];
    wb_ei_ecef = Tr_ei*wb_ei_eci;
    wb_bi_body = wb_bn_body + wb_ne_body + Tr_bi*wb_ei_eci;

    ## Transport Rate
    angle_rates = Angle_Rates( wb_bi_body, wb_bn_body, wb_ei_eci, wb_ei_ecef, wb_ne_eci )
    return angle_rates
end

############################################################################
# Initialize
function Airframe_Initialize( tsim::Float64,
                     time::Float64, dt::Float64,
                     position::WGS_Pos, velocity::Velocity_States,
                     euler::Euler_Angles, angle_rates::Angle_Rates, mass_props::MassProps)

    ## Time Computations
    dt_eci = time - tsim;

    ## Transformation Matrices
    Tr_ei = position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = position.Tr_ne;
    Tr_en = Tr_ne';
    Tr_bn = euler.Tr_bn;
    Tr_nb = Tr_bn';
    Tr_bi = Tr_bn*Tr_ne*Tr_ei;
    Tr_ib = Tr_bi';

    quat = QuatInit( Tr_bi );
    quat_bi = Quaternions(quat[1], quat[2], quat[3], quat[4] );

    zz = [0.0;0.0;0.0];

    xo = [ position.rm_eci; velocity.vm_eci; angle_rates.wb_bi_body; quat];

    ## Quad Rotor Intialization
    airframe_init = Aero6DOF(time, dt, position, velocity, euler, angle_rates, zz,  mass_props, quat_bi, xo);

    return airframe_init
end

############################################################################
# Update
function Airframe_Update(airframe::Aero6DOF,
                ftot_body::Array{Float64,1},
                mtot_body::Array{Float64,1},
                )
    dt = airframe.dt;
    airframe_out = airframe;

    ## Time
    time = airframe.time
    wb_ei_eci = [0.0;0.0;WEII3];
    wb_ei_ecef = wb_ei_eci;

    ## Transformation Matrices
    Tr_ei = airframe.position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = airframe.position.Tr_ne;
    Tr_en = Tr_ne';
    Tr_bn = airframe.euler.Tr_bn;
    Tr_nb = Tr_bn';
    Tr_bi = Tr_bn*Tr_ne*Tr_ei;
    Tr_ib = Tr_bi';

    ## Mass Properties
    mass = airframe.massprops.mass;
    II = airframe.massprops.II;
    cg_vec_body = airframe.massprops.cg;

    ## Quadrotor Velocity
    vel = airframe.velocity;
    vm_eci = airframe.velocity.vm_eci;
    vm_ecef = Tr_ei*vm_eci - cross(wb_ei_ecef, airframe.position.rm_ecef);
    vm_ned = Tr_ne*vm_ecef;
    vm_body = Tr_bn*vm_ned;

    ## Quadrotor Position
    pos = airframe.position;

    ## Attitude
    wb_body = airframe.angle_rates;
    wb_bi_body = wb_body.wb_bi_body;

    ## Quaternions
    quat_bi = airframe.quat_bi;
    quat_bi_vec = [quat_bi.q0; quat_bi.q1; quat_bi.q2; quat_bi.q3];

    ## State Vector
    xo = [ airframe.position.rm_eci;
           airframe.velocity.vm_eci;
           airframe.angle_rates.wb_bi_body;
           quat_bi_vec;
           ];

    ## Gravity
    grav_ecef = WGS84_GRAVITY(airframe.position.rm_ecef);
    grav_eci = Tr_ie*grav_ecef - cross( wb_ei_eci, cross(wb_ei_eci, airframe.position.rm_eci) );
    grav_ecef = Tr_ei*grav_eci;
    grav_body = Tr_bi*grav_eci;
    grav_ned = Tr_nb*grav_body;
    fgrav_body = mass*grav_body;

    ## Total Forces in Body Frame
    acc_meas_body = ftot_body/mass;
    ftot_eci = Tr_ib*ftot_body;
    airframe_out.accel_meas = acc_meas_body;

    ## Computing Derivatives
    dx_pos = vm_eci;
    dx_vel = Tr_ib*(1.0/mass*( ftot_body + fgrav_body - cross( wb_bi_body, mass*vm_body ) ) );
    dx_wb  = II^-1*( mtot_body - cross( wb_bi_body, II*wb_bi_body ) );
    dx_quat = QuatDerivative(quat_bi_vec, wb_bi_body)

    ## Derivative Vector
    dx = [dx_pos; dx_vel; dx_wb; dx_quat];

    ## Euler Integration
    xo = xo + dt*dx;
    time = time + dt;

    ## Updating Parameters
    airframe_out.xo = xo;
    airframe_out.time = time;

    ## Updating States
    rm_eci = xo[1:3];
    vm_eci = xo[4:6];
    wb_bi_body = xo[7:9];
    quat_bi_body = xo[10:13];

    ## Updating Quaternions
    quat_bi_body = quat_bi_body/norm(quat_bi_body);
    airframe_out.quat_bi = Quaternions(quat_bi_body[1],
                                   quat_bi_body[2],
                                   quat_bi_body[3],
                                   quat_bi_body[4] );

    ## Updating Transformation Matrix ECI to ECEF;
    Tr_ei = TR_EI(time);
    Tr_ie = Tr_ei';

    ## Updating Positon
    rm_ecef = Tr_ei*rm_eci;
    LLH = ECEF2LLH(rm_ecef);
    latitude = LLH[1];
    longitude = LLH[2];
    altitude = LLH[3];
    pos = InitPos( time, latitude, longitude, altitude );
    airframe_out.position = pos;

    ## Updating Transformations
    quat_bi_body = quat_bi_body/norm(quat_bi_body);
    Tr_bi = QuatToDCM(quat_bi_body);
    Tr_ib = Tr_bi';
    Tr_ne = pos.Tr_ne;
    Tr_en = Tr_ne';
    Tr_bn = Tr_bi*Tr_ie*Tr_en;
    Tr_nb = Tr_bn';
    Tr_eb = Tr_en*Tr_nb;
    Tr_be = Tr_eb';

    ## Updating Euler
    eulervec = ComputeEuler( Tr_bn );
    roll = eulervec[1];
    pitch = eulervec[2];
    heading = eulervec[3];
    airframe_out.euler = InitEuler( roll, pitch, heading );

    ## Updating Body Rates
    wb_bi_ecef = Tr_eb*wb_bi_body;
    wb_ei_ecef = wb_ei_eci;
    wb_ei_body = Tr_be*wb_ei_ecef;
    wb_bn_body = wb_bi_body - wb_ei_body;
    airframe_out.angle_rates = InitAngleRates(wb_bn_body, airframe_out.position, airframe_out.euler);

    ## Updating Velocity
    vm_ecef = Tr_ei*vm_eci - cross(wb_ei_ecef, pos.rm_ecef);
    vm_ned = pos.Tr_ne*vm_ecef;
    vmag = norm(vm_ned);
    gamma = atan( -vm_ned[3], sqrt(vm_ned[1]^2 + vm_ned[2]^2));
    azimuth = atan( vm_ned[2], vm_ned[1] );
    vel = InitVel( vmag, gamma, azimuth, pos );
    airframe_out.velocity = vel;

    return airframe_out;
end


text:
using AbstractPlotting

using Test

include("quaternions.jl")

# if get(ENV, "IS_TRAVIS_CI", "false") == "false"
#   exit(0);
# end

## begin CI-only testing
# if get(ENV, "IS_TRAVIS_CI", "false") == "false"
#   exit(0);
# end

 ## begin CI-only testing=

using MakieGallery

 # Load the database
database = MakieGallery.load_database()

## Exceptions are made on basis of title and not index,
 # because index may change as MakieGallery changes.

# All these require FFMpeg and need to save,
# and are therefore ignored
ffmpeg_exs  = [
            "Animation", "Lots of Heatmaps","Animated Scatter",
            "Chess Game", "Record Video",  "Animated surface and wireframe",
            "Moire", "pong", "pulsing marker", "Travelling wave",
            "Type recipe for molecule simulation"
]

# All these require GLMakie and so are ignored
glmakie_exs = [
            "Textured Mesh", "Load Mesh", "Wireframe of a Mesh",
            "FEM mesh 3D", "Normals of a Cat", "Line GIF"
]

# Requires GDAL (a GL package) so ignored
gdal_exs = [
        "WorldClim visualization"
]

 # Requires GLMakie and ModernGL, so ignored
moderngl_exs = [
            "Explicit frame rendering"
]

# use Stepper plots, which save, so ignored
save_exs = [
          "Axis theming", "Labels", "Color Legend",
          "Stepper demo"
]

# hopefullly fixed by next tag of MakieGallery (already in AbstractPlotting
color_exs = ["colormaps"]

# curl fails with this for some reason, so it has been ignored.
curl_exs = ["Earth & Ships"]

# combine all exceptions into a single Set
exc_str = union(ffmpeg_exs, glmakie_exs, gdal_exs, moderngl_exs, save_exs, color_exs, curl_exs)

# iterate over database
for i in 1:length(database)

   # skip if the title is in the list of exceptions
  if database[i].title ∈ exc_str

     print("Skipping " * database[i].title * "\n(removed from tests explicitly)\n")

     continue  # skip the entry

   end

   print("Running " * database[i].title * "\n(index $i)\n")

   MakieGallery.eval_example(database[i]);  # evaluate the entry

   end


text:
## Thanks to randy3k for pointing this out,
## `RCall` needs to be precompiled with the current R.
## <https://github.com/Non-Contradiction/JuliaCall/issues/9>
## as well as coming up with the solution

CurrentRhome = normpath(ARGS[1])

## println(Rhome)

ENV["R_HOME"] = CurrentRhome

const julia07 = VERSION > v"0.6.5"

if julia07
    using Pkg
else
    macro pkg_str(x) begin nothing end end
end

function installed(name)
    @static if julia07
        get(Pkg.installed(), name, nothing)
    else
        Pkg.installed(name)
    end
end

if installed("Suppressor") == nothing
    Pkg.add("Suppressor")
end;

using Suppressor

if installed("RCall") == nothing
    Pkg.add("RCall")
end;

try
    using RCall
catch e
    Pkg.build("RCall")
end

if normpath(RCall.Rhome) != CurrentRhome
    if installed("RCall") >= v"0.10.2"
        Pkg.build("RCall")
    else
        Base.compilecache("RCall")
    end
end

max_stars_repo_name: pthariensflame/JuliaCall

max_issues_repo_name: pthariensflame/JuliaCall

max_forks_repo_name: pthariensflame/JuliaCall


text:
"""
    linearfilter(N::BinaryNetwork, α::Vector{Float64}=[0.25, 0.25, 0.25, 0.25])

Given a network `N` compute the linear filter scores according to Stock et al.
(2017). High scores for negative interactions indicate potential false negative
or missing interactions. Though this it returned as a probabilistic network,
score do not necessary convey a probabilistic interpretation.

The values of α give the *relative* weight of, in order, the measured
interaction value, the out-degree of the species, the in-degree of the species,
and of network connectance. The default parameterization is to have all four at
equal importance.

#### References

Stock, M., Pahikkala, T., Airola, A., Waegeman, W., Baets, B.D., 2018. Algebraic
Shortcuts for Leave-One-Out Cross-Validation in Supervised Network Inference.
bioRxiv 242321. https://doi.org/10.1101/242321

Stock, M., Poisot, T., Waegeman, W., Baets, B.D., 2017. Linear filtering reveals
false negatives in species interaction data. Scientific Reports 7, 45908.
https://doi.org/10.1038/srep45908
"""
function linearfilter(N::T; α::Vector{Float64}=fill(1.0, 4)) where {T <: BinaryNetwork}
    @assert length(α) == 4
    @assert all(α .≥ 0.0)
    α = α ./ sum(α) # This ensures that α sums to 1.0

    # Get the new weights
    W = α[1]N.edges .+ α[2]mean(N.edges; dims=2) .+ α[3]mean(N.edges; dims=1) .+ α[4]mean(N.edges)

    # Prepare a return object
    return_type = T <: AbstractBipartiteNetwork ? BipartiteProbabilisticNetwork : UnipartiteProbabilisticNetwork
    F = return_type(W, EcologicalNetworks._species_objects(N)...)

    return F
end

"""
    linearfilterzoo(N::BinaryNetwork, α::Vector{Float64}=[0.25, 0.25, 0.25, 0.25])

Compute the zero-one-out version of the linear filter (`linearfilter`), i.e.
the score for each interaction if that interaction would not occur in the network.
For example, if N[4, 6] = 1 (interaction between species 4 and 6), the result at
postion (4, 6) is the score of the filter using a network in which that interaction
did not occur. This function is useful for validating the filter whether it can
detect false negative (missing) interactions.

#### References

Stock, M., Pahikkala, T., Airola, A., Waegeman, W., Baets, B.D., 2018. Algebraic
Shortcuts for Leave-One-Out Cross-Validation in Supervised Network Inference.
bioRxiv 242321. https://doi.org/10.1101/242321

Stock, M., Poisot, T., Waegeman, W., Baets, B.D., 2017. Linear filtering reveals
false negatives in species interaction data. Scientific Reports 7, 45908.
https://doi.org/10.1038/srep45908
"""
function linearfilterzoo(N::T; α::Vector{Float64}=fill(0.25, 4)) where {T <: BinaryNetwork}
    F = linearfilter(N, α=α)  # depart from scores of the filter
    Fzoo = F  # update the values to zero-one-out scores, reference is to keep
    # the code readable

    # Get the size of the network
    n = richness(N; dims=1)
    m = richness(N; dims=2)
    # Get probabilities
    for s1 in species(N; dims=1)
        for s2 in species(N; dims=2)
            t_val = F[s1,s2] - (α[1] + (α[2] / m) + (α[3] / n) + α[4] / (n * m)) * N[s1,s2]
            Fzoo[s1,s2] = max(min(t_val, one(t_val)), zero(t_val))
        end
    end

    return Fzoo
end


text:
module Bazykin85

using DiffEqBase: ODEProblem
using Parameters: @with_kw, @unpack
using StaticArrays: SVector
using Setfield: @lens
import Setfield

using ...Bifurcations: BifurcationProblem
using ...Codim2: DiffEqCodim2Problem

@with_kw struct Bazykin85Param{A, E, G, D}
    α::A = 0.1
    ϵ::E = 0.01  # ≪ 1 ?
    γ::G = 1.0
    δ::D = 0.1
end

function f(x::SVector, p, t)
    @unpack α, ϵ, γ, δ = p
    return SVector(
             x[1] - x[1] * x[2] / (1 + α * x[1]) - ϵ * x[1]^2,
        -γ * x[2] + x[1] * x[2] / (1 + α * x[1]) - δ * x[2]^2,
    )
end

function f(du, u, p, t)
    du .= f(SVector{2}(u), p, t)
    nothing
end

make_prob(
        p = Bazykin85Param();
        u0 = SVector(1 / p.ϵ, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis = (@lens _.α),
        t_domain = (0.01, 1.5),
        kwargs...) =
    BifurcationProblem(ode, param_axis, t_domain;
                       phase_space = (SVector(-0.1, -0.1),  # u_min
                                      SVector(Inf, Inf)),   # u_max
                       kwargs...)

prob = make_prob()
ode = prob.p.de_prob
u0 = ode.u0
param_axis = prob.p.param_axis
t_domain = prob.t_domain

function make_codim2_prob(
        p = Bazykin85Param();
        u0 = SVector(1 / p.ϵ, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis1 = (@lens _.α),
        param_axis2 = (@lens _.δ),
        t_domain = ([0.01, 0.0], [1.5, 10.0]),
        kwargs...)
    return DiffEqCodim2Problem(
        ode,
        param_axis1,
        param_axis2,
        t_domain;
        kwargs...)
end

end  # module


text:
module BlueNoise

function _parse_generated_data(
    (raw_data, dimensions)::Tuple{Base.CodeUnits{UInt8},NTuple{D,Int}}
)::Array{Float64,D} where {D}
    result = Array{Float64,D}(undef, dimensions)
    indices = prod(dimensions)
    digit_values = Dict(
        b"0"[1] => 0,
        b"1"[1] => 1,
        b"2"[1] => 2,
        b"3"[1] => 3,
        b"4"[1] => 4,
        b"5"[1] => 5,
        b"6"[1] => 6,
        b"7"[1] => 7,
        b"8"[1] => 8,
        b"9"[1] => 9,
        b"A"[1] => 0xA,
        b"B"[1] => 0xB,
        b"C"[1] => 0xC,
        b"D"[1] => 0xD,
        b"E"[1] => 0xE,
        b"F"[1] => 0xF,
    )
    for index in 1:indices
        int_val = 0
        for offset in 1:16
            hex_digit = raw_data[(index - 1) * 16 + offset]
            int_val = (int_val * 0x10) + digit_values[hex_digit]
        end
        result[index] = int_val / Float64(indices ÷ last(dimensions))
    end
    result
end

const PREGENERATED_1D4V_RES128 = _parse_generated_data(include("generated_1d4v_res128.jl"))
const PREGENERATED_1D4V_RES1024 = _parse_generated_data(include("generated_1d4v_res1024.jl"))
const PREGENERATED_1D4V_RES8192 = _parse_generated_data(include("generated_1d4v_res8192.jl"))
const PREGENERATED_1D4V_RES32768 = _parse_generated_data(include("generated_1d4v_res32768.jl"))
const PREGENERATED_1D4V = PREGENERATED_1D4V_RES32768

const PREGENERATED_2D4V_RES32 = _parse_generated_data(include("generated_2d4v_res32.jl"))
const PREGENERATED_2D4V_RES128 = _parse_generated_data(include("generated_2d4v_res128.jl"))
const PREGENERATED_2D4V_RES512 = _parse_generated_data(include("generated_2d4v_res512.jl"))
const PREGENERATED_2D4V_RES1024 = _parse_generated_data(include("generated_2d4v_res1024.jl"))
const PREGENERATED_2D4V = PREGENERATED_2D4V_RES1024

const PREGENERATED_3D4V_RES16 = _parse_generated_data(include("generated_3d4v_res16.jl"))
const PREGENERATED_3D4V_RES32 = _parse_generated_data(include("generated_3d4v_res32.jl"))
const PREGENERATED_3D4V_RES64 = _parse_generated_data(include("generated_3d4v_res64.jl"))
const PREGENERATED_3D4V = PREGENERATED_3D4V_RES64

blue_noise_1d4v(x::Int)::Vector{Float64} = PREGENERATED_1D4V[
    mod1(x, size(PREGENERATED_1D4V)[1]),
    1:4
]

blue_noise_1d3v(x::Int)::Vector{Float64} = PREGENERATED_1D4V[
    mod1(x, size(PREGENERATED_1D4V)[1]),
    1:3
]

blue_noise_1d2v(x::Int)::Vector{Float64} = PREGENERATED_1D4V[
    mod1(x, size(PREGENERATED_1D4V)[1]),
    1:2
]

blue_noise_1d1v(x::Int)::Float64 = PREGENERATED_1D4V[
    mod1(x, size(PREGENERATED_1D4V)[1]),
    1
]

blue_noise_2d4v(x::Int, y::Int)::Vector{Float64} = PREGENERATED_2D4V[
    mod1(x, size(PREGENERATED_2D4V)[1]),
    mod1(y, size(PREGENERATED_2D4V)[2]),
    1:4
]

blue_noise_2d3v(x::Int, y::Int)::Vector{Float64} = PREGENERATED_2D4V[
    mod1(x, size(PREGENERATED_2D4V)[1]),
    mod1(y, size(PREGENERATED_2D4V)[2]),
    1:3
]

blue_noise_2d2v(x::Int, y::Int)::Vector{Float64} = PREGENERATED_2D4V[
    mod1(x, size(PREGENERATED_2D4V)[1]),
    mod1(y, size(PREGENERATED_2D4V)[2]),
    1:2
]

blue_noise_2d1v(x::Int, y::Int)::Float64 = PREGENERATED_2D4V[
    mod1(x, size(PREGENERATED_2D4V)[1]),
    mod1(y, size(PREGENERATED_2D4V)[2]),
    1
]

blue_noise_3d4v(x::Int, y::Int, z::Int)::Vector{Float64} = PREGENERATED_3D4V[
    mod1(x, size(PREGENERATED_3D4V)[1]),
    mod1(y, size(PREGENERATED_3D4V)[2]),
    mod1(z, size(PREGENERATED_3D4V)[3]),
    1:4
]

blue_noise_3d3v(x::Int, y::Int, z::Int)::Vector{Float64} = PREGENERATED_3D4V[
    mod1(x, size(PREGENERATED_3D4V)[1]),
    mod1(y, size(PREGENERATED_3D4V)[2]),
    mod1(z, size(PREGENERATED_3D4V)[3]),
    1:3
]

blue_noise_3d2v(x::Int, y::Int, z::Int)::Vector{Float64} = PREGENERATED_3D4V[
    mod1(x, size(PREGENERATED_3D4V)[1]),
    mod1(y, size(PREGENERATED_3D4V)[2]),
    mod1(z, size(PREGENERATED_3D4V)[3]),
    1:2
]

blue_noise_3d1v(x::Int, y::Int, z::Int)::Float64 = PREGENERATED_3D4V[
    mod1(x, size(PREGENERATED_3D4V)[1]),
    mod1(y, size(PREGENERATED_3D4V)[2]),
    mod1(z, size(PREGENERATED_3D4V)[3]),
    1
]

export blue_noise_1d1v, blue_noise_1d2v, blue_noise_1d3v, blue_noise_1d4v
export blue_noise_2d1v, blue_noise_2d2v, blue_noise_2d3v, blue_noise_2d4v
export blue_noise_3d1v, blue_noise_3d2v, blue_noise_3d3v, blue_noise_3d4v

end # module

max_stars_repo_name: joshua-maros/julia-blue-noise

max_issues_repo_name: joshua-maros/julia-blue-noise

max_forks_repo_name: joshua-maros/julia-blue-noise


text:
using Plots
# pyplot()
do_viz = true

@views function heat_1D_nprocs()
    # Physics
    lx  = 10.0
    λ   = 1.0
    ρCp = 1.0
    nt  = 200
    # Numerics
    np  = 4             # number of procs
    nx  = 32            # local number of gridpoints
    # Derived numerics
    nxg = (nx-2)*np+2   # global number of grid points
    dxg = lx/nxg
    dt  = dxg^2/ρCp/λ/2.1
    # Array allocation
    x   = zeros(nx,np)  # local coord array
    T   = zeros(nx,np)  # local Temp array
    xt  = zeros(nxg)    # global coord array
    Tt  = zeros(nxg)    # global initial Temp array
    Tg  = zeros(nxg)    # global Temp array
    # Initial condition
    for ip = 1:np
        i1 = 1 + (ip-1)*(nx-2)
        for ix = 1:nx
            x[ix,ip] = ( (ip-1)*(nx-2) + (ix-0.5) )*dxg - 0.5*lx
            T[ix,ip] = exp(-x[ix,ip]^2)
        end
        xt[i1:i1+nx-2] .= x[1:end-1,ip]; if (ip==np) xt[i1+nx-1] = x[end,ip] end
        Tt[i1:i1+nx-2] .= T[1:end-1,ip]; if (ip==np) Tt[i1+nx-1] = T[end,ip] end
    end
    # Time loop
    for it = 1:nt
        for ip = 1:np # compute physics locally
            T[2:end-1,ip] .= T[2:end-1,ip] .+ dt.*λ./ρCp.*diff(diff(T[:,ip])./dxg)./dxg
        end
        for ip = 1:np-1 # update boundaries
            T[end,ip  ] = T[    2,ip+1]
            T[  1,ip+1] = T[end-1,ip  ]
        end
        for ip = 1:np # global picture
            i1 = 1 + (ip-1)*(nx-2)
            Tg[i1:i1+nx-2] .= T[1:end-1,ip]
        end
        # Visualise
        if do_viz
            plot(xt, Tt, legend=false, linewidth=5, xlabel="lx", ylabel="heat")
            display(plot!(xt, Tg, legend=false, linewidth=5, framestyle=:box, xlabel="lx", ylabel="heat", title="diffusion it=$(it)"))
        end
    end
end

@time heat_1D_nprocs()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
using QuantumComputing
using Base.Test

# write your own tests here
a = Ket(0,1,0,1,0)
b = Ket(1,0,1,0,1)

n = length(a)
@test n==5
@test a == a
@test a != b

c = Swap(n,1,1)*a
@test c == a
c = Swap(n,1,2)*a
@test c == Ket(1,0,0,1,0)
c = Swap(n,1,n)*a
@test c == a

c = Swap(2,2,2)
for i=0:1, j=0:1
@test c*Ket(i,j)==Ket(i,j)
end

c = Not(5,1)
@test c*Ket(0,0,0,0,0)==Ket(1,0,0,0,0)
c = Not(n,n)
@test c*Ket(1,1,1,1,1)==Ket(1,1,1,1,0)

a = Ket(1,1)
c = CNot(2,1,2)*a
println(c)
@test c == Ket(0,1)

c = GCNot(2,2,1)
for i=0:1, j=0:1
@test c*Ket(i,j) == GGate(CNOT)*Ket(i,j)
end


text:
### A Pluto.jl notebook ###
# v0.12.15

using Markdown
using InteractiveUtils

# ╔═╡ 40d8ceb2-4b66-11eb-19b7-b52e55c7659b
begin
	using ImageInpainting,FileIO,Images,ImageDraw,Plots
	using TestImages
end

# ╔═╡ bd22a0d0-4b67-11eb-0e58-974d11ff83d0
image=load("C:\\Tempy\\2077.jpg")

# ╔═╡ c7c76710-4b6b-11eb-3c08-a124a24b84aa
split=channelview(image);

# ╔═╡ 7bfb4cb0-4b6c-11eb-0c4a-cdfbd7020989
R=colorview(RGB,StackedView(split[1,:,:], zeroarray, zeroarray))

# ╔═╡ b97c14c0-4b6c-11eb-1fb8-015c71819eec
G=colorview(RGB,StackedView(zeroarray,split[2,:,:], zeroarray))

# ╔═╡ ba26fac0-4b6c-11eb-33f9-254abd8cedf7
B=colorview(RGB,StackedView(zeroarray, zeroarray,split[3,:,:] ))

# ╔═╡ 3954f680-4b6d-11eb-3878-ad32170054eb
masking=draw(image,Ellipse(CirclePointRadius(140,666,90)))

# ╔═╡ ba6df960-4b6d-11eb-2a58-41e7a9fd7140
mask= masking.!=image

# ╔═╡ f2727d40-4b6d-11eb-2c0d-f7711c97884c
begin
	function doIt(image)
	split=channelview(image);
	r=inpaint(convert(Array{Float64},(Gray.(split[1,:,:]))),mask,Criminisi(150,130))
	g=inpaint(convert(Array{Float64},(Gray.(split[2,:,:]))),mask,Criminisi(150,110))
	b=inpaint(convert(Array{Float64},(Gray.(split[3,:,:]))),mask,Criminisi(150,130))
	return colorview(RGB,StackedView(r,g,b))
	end
end

# ╔═╡ 22e3d8f0-4b70-11eb-018f-37c28b93a851
doIt(image)

# ╔═╡ Cell order:
# ╠═40d8ceb2-4b66-11eb-19b7-b52e55c7659b
# ╠═bd22a0d0-4b67-11eb-0e58-974d11ff83d0
# ╠═c7c76710-4b6b-11eb-3c08-a124a24b84aa
# ╠═7bfb4cb0-4b6c-11eb-0c4a-cdfbd7020989
# ╠═b97c14c0-4b6c-11eb-1fb8-015c71819eec
# ╟─ba26fac0-4b6c-11eb-33f9-254abd8cedf7
# ╟─3954f680-4b6d-11eb-3878-ad32170054eb
# ╠═ba6df960-4b6d-11eb-2a58-41e7a9fd7140
# ╠═f2727d40-4b6d-11eb-2c0d-f7711c97884c
# ╠═22e3d8f0-4b70-11eb-018f-37c28b93a851

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc