text:
using LinearAlgebra
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
using CUDA

include("part2_utils.jl")
include("krylov.jl")

""" Iterative methods that can be used as coarse solvers in Multigrid. """
@enum CoarseSolver_t begin
    jacobi
    conjugate_gradient
end

""" Input parameters for Multigrid. """
mutable struct MGOpt
    coarse_solve_size :: Int
    coarse_solver :: CoarseSolver_t
    execution_policy :: ExecutionPolicy_t

    MGOpt() = new(5, jacobi, parallel_shmem)
end

""" Preallocates all buffers/arrays required for all levels of Multigrid. """
function preallocate_buffers(nx, ny) :: Union{Dict{Symbol, Dict{Int, AbstractArray}}, Nothing}
     # prealloc all the buffers
     λx = (nx > ny) ? (nx-1) ÷ (ny-1) : 1
     λy = (ny > nx) ? (ny-1) ÷ (nx-1) : 1
     k_max = Int(log2(min(nx, ny)-1))
     prealloc_dict = Dict(
         :res_f   => Dict( k => @zeros(λx*2^k    +1, λy*2^k    +1) for k in 1:k_max ),
         :corr_f  => Dict( k => @zeros(λx*2^k    +1, λy*2^k    +1) for k in 1:k_max ),
         :corr_c  => Dict( k => @zeros(λx*2^(k-1)+1, λy*2^(k-1)+1) for k in 1:k_max ),
         :res_c   => Dict( k => @zeros(λx*2^(k-1)+1, λy*2^(k-1)+1) for k in 1:k_max ),
    )

    return prealloc_dict
end

""" Solves (∇^2 - c) u = f using matrix-free V-cycle Multigrid. """
function MGsolve_2DPoisson!(u::AbstractArray{Float64}, f::AbstractArray{Float64}, h::Float64, c::Float64, tol::Float64, niters::Int, apply_BCs::Bool; opt=MGOpt(), verbose=false, prealloc_dict::Union{Dict{Symbol, Dict{Int, AbstractArray}}, Nothing} = nothing)
    # set nx, ny
    nx, ny = size(u)

    @assert opt.coarse_solve_size <= min(nx, ny)
    @assert isinteger(log2(opt.coarse_solve_size - 1))

    # prealloc all the buffers
    if prealloc_dict == nothing
        prealloc_dict = preallocate_buffers(nx, ny)
    end

    f_rms = sqrt(sum(f.^2)/(nx*ny))

    tolf = tol * f_rms

    r_rms = 0.0
    for iter = 1:niters
        # apply Dirichlet and Neumann BCs for temperature T
        if (apply_BCs)
            apply_boundary_conditions!(u)
        end
        # execute V-cycle iteration
        r_rms = Vcycle_2DPoisson!(u, f, h, c, tol, opt.coarse_solve_size, opt.coarse_solver, opt.execution_policy, apply_BCs, prealloc_dict=prealloc_dict)
        @synchronize()

        if verbose
            println("$(iter) $(r_rms / f_rms)")
        end
        if (r_rms < tolf)
            if verbose
                println("V-cycle multigrid converged in $(iter) iterations.")
            end
            break
        end
    end

    if (r_rms > tolf)
        @warn "V-cycle multigrid failed to converge within" niters "iterations."
    end

    return r_rms

end

""" Runs a single V-cycle iteration of Multigrid.

Note: nx, ny must be = (2^k)+1, where k ∈ N.

"""
function Vcycle_2DPoisson!(u_f::AbstractArray{Float64}, rhs::AbstractArray{Float64}, h::Float64, c::Float64, tol::Float64, coarse_solve_size::Int, coarse_solver::CoarseSolver_t, execution_policy::ExecutionPolicy_t, apply_BCs::Bool; prealloc_dict::Union{Dict{Symbol, Dict{Int, AbstractArray}}, Nothing} = nothing)

    nx, ny = size(u_f)

    if ( nx-1 != 2*((nx-1)÷2) || ny-1 != 2*((ny-1)÷2) )
        error("ERROR:not a power of 2")
    end

    # coarse grid size
    nxc = 1+(nx-1)÷2
    nyc = 1+(ny-1)÷2

    k_curr = Int(log2(min(nx, ny)-1))

    if prealloc_dict == nothing
        prealloc_dict = preallocate_buffers(nx, ny)
    end

    res_f   = prealloc_dict[:res_f][k_curr]
    corr_f  = prealloc_dict[:corr_f][k_curr]
    corr_c  = prealloc_dict[:corr_c][k_curr]
    res_c   = prealloc_dict[:res_c][k_curr]

    res_f   .= 0.0
    corr_f  .= 0.0
    corr_c  .= 0.0
    res_c   .= 0.0

    res_rms = 0.0

    if (min(nx, ny) > coarse_solve_size)   # not the coarsest level

        #---------- take 2 iterations on the fine grid--------------
        res_rms = iteration_2DPoisson!(u_f, rhs, h, c, res_f, execution_policy)
        res_rms = iteration_2DPoisson!(u_f, rhs, h, c, res_f, execution_policy)

        #--------- restrict the residual to the coarse grid --------
        residual_2DPoisson_wrapper!(u_f, rhs, h, c, res_f, execution_policy)
        restrict_wrapper!(res_f, res_c, apply_BCs, execution_policy)

        #---------- solve for the coarse grid correction -----------
        corr_c .= 0.
        res_rms = Vcycle_2DPoisson!(corr_c, res_c, h*2, c, tol, coarse_solve_size, coarse_solver, execution_policy, apply_BCs) # *RECURSIVE CALL*

        #---- prolongate (interpolate) the correction to the fine grid
        prolongate_wrapper!(corr_c, corr_f, apply_BCs, execution_policy)

        #---------- correct the fine-grid solution -----------------
        u_f .= u_f - corr_f

        #---------- two more smoothing iterations on the fine grid---
        res_rms = iteration_2DPoisson!(u_f, rhs, h, c, res_f, execution_policy)
        res_rms = iteration_2DPoisson!(u_f, rhs, h, c, res_f, execution_policy)

    else

        if coarse_solver == jacobi
            #----- coarsest level (ny=5): iterate to get 'exact' solution
            coarse_solve_iters = 20*coarse_solve_size
            tol_rhs = tol * sqrt(sum(rhs.^2)/(nx*ny))
            for i = 1:coarse_solve_iters  # heuristic
                res_rms = iteration_2DPoisson!(u_f, rhs, h, c, res_f, execution_policy)
                if res_rms < tol_rhs
                    break
                end
            end
            if res_rms > tol_rhs
                @debug "Coarse solve of V-cycle multigrid did not converge to " tol " within " coarse_solve_iters "iterations. This is usually fine as long as it converges in a fixed amount of Vcycles."
            end
        elseif coarse_solver == conjugate_gradient
            coarse_solve_iters = 20*coarse_solve_size
            res_rms = cg!(u_f, rhs, h, h, c, tol, coarse_solve_iters, execution_policy=execution_policy)
        else
            error()
        end

    end

    return res_rms   # returns the rms. residual
end

""" Computes the residual res = (∇^2 - c) u - f, without using shared memory. """
@parallel_indices (ix, iy) function residual_2DPoisson!(u::AbstractArray{Float64}, f::AbstractArray{Float64}, h::Float64, c::Float64, res::AbstractArray{Float64})
    nx, ny = size(u)
    C = (4.0 + c * h^2)
    _h2 = 1 / h^2
    if (1 < ix < nx && 1 < iy < ny)
        res[ix, iy] = (
                       (  u[ix+1, iy]
                        + u[ix-1, iy]
                        + u[ix, iy+1]
                        + u[ix, iy-1]
                        - C * u[ix, iy] ) * _h2
                       - f[ix, iy]
                      )
    end
    return nothing
end

""" Computes the residual res = (∇^2 - c) u - f, using shared memory. """
@parallel_indices (ix, iy) function residual_2DPoisson_shmem!(u_::AbstractArray{Float64}, f::AbstractArray{Float64}, h::Float64, c::Float64, res::AbstractArray{Float64})
    nx, ny = size(u_)
    C = (4.0 + c * h^2)
    _h2 = 1 / h^2

    tx = @threadIdx().x + 1
    ty = @threadIdx().y + 1
    u = @sharedMem(eltype(u_), (@blockDim().x+2, @blockDim().y+2))
    if (1 <= ix <= nx && 1 <= iy <= ny)
        u[tx, ty] = u_[ix, iy]
    end
    @sync_threads()
    if (1 < ix < nx && 1 < iy < ny)
        u[tx, ty] = u_[ix, iy]
        if (@threadIdx().x == 1)             u[tx-1, ty] = u_[ix-1, iy] end
        if (@threadIdx().x == @blockDim().x) u[tx+1, ty] = u_[ix+1, iy] end
        if (@threadIdx().y == 1)             u[tx, ty-1] = u_[ix, iy-1] end
        if (@threadIdx().y == @blockDim().y) u[tx, ty+1] = u_[ix, iy+1] end
        @sync_threads()
        res[ix, iy] = (
                       (  u[tx+1, ty]
                        + u[tx-1, ty]
                        + u[tx, ty+1]
                        + u[tx, ty-1]
                        - C * u[tx, ty] ) * _h2
                       - f[ix, iy]
                      )
    end
    return nothing
end

""" Wrapper function computing residual res = (∇^2 - c) u - f, for various execution policies. """
function residual_2DPoisson_wrapper!(u_f, rhs, h, c, res_f, execution_policy)
    nx, ny = size(u_f)

    if execution_policy == parallel_shmem
        threads = (32, 8)
        blocks = ((nx, ny) .- 1) .÷ (threads) .+ 1
        shmem = prod(threads.+2) * sizeof(Float64)
        @parallel blocks threads shmem=shmem residual_2DPoisson_shmem!(u_f, rhs, h, c, res_f)
    elseif execution_policy == parallel
        @parallel residual_2DPoisson!(u_f, rhs, h, c, res_f)
    elseif execution_policy == serial
        error() # Not implemented yet
    else
        error()
    end
end

""" Performs a single Jacobi iteration on field u.

Note: alpha = 4.0 / 5.0 is the optimal damping factor for smoothing using Jacobi iterations.

"""
function iteration_2DPoisson!(u, f, h, c, res, execution_policy; alpha = 4.0 / 5.0)
    nx, ny = size(u)

    # Do one Jacobi iteration
    residual_2DPoisson_wrapper!(u, f, h, c, res, execution_policy)

    # update r_rms
    r_rms = sqrt(sum(res.^2) / (nx * ny))

    # update u
    u .+= alpha * (h^2 / (4.0 + c * h^2)) .* res

    return r_rms
end

""" Performs a single Gauss-Seidel iteration on field u.

Note: alpha = 1.0 is the optimal damping factor for smoothing using Gauss-Seidel iterations.

Gauss-Seidel is highly serial so in our current Multigrid implementation we replaced it by Jacobi.

We only use this function for debugging.

"""
@views function iteration_2DPoisson_gs!(u, f, h, c; alpha = 1.0)

    nx, ny = size(u)

    # Gauss-Seidel iteration
    r_rms = 0.0
    for j = 2:ny-1
        for i = 2:nx-1

            # compute residual
            r = ( u[i+1, j]
                 + u[i-1, j]
                 + u[i, j+1]
                 + u[i, j-1]
                 - (4.0 + c * h^2) * u[i, j] ) / h^2 - f[i, j]

            # update u
            u[i, j] = u[i, j] + alpha * (h^2 / (4.0 + c * h^2)) * r

            # update r_rms
            r_rms = r_rms + r^2

        end
    end

    r_rms = sqrt(r_rms / (nx * ny))

    return r_rms
end

""" Restricts a field from a fine grid to a coarse grid of half the resolution.

This is a serial version, mainly used for debugging.

"""
@views function  restrict_serial!(fine, coarse, apply_BCs)
    nx = size(fine, 1)
    ny = size(fine, 2)

    # initalize coarse to zero (+ Dirichlet(0) BCs)
    coarse .= 0.0

    # apply restriction stencil
    jc = 2
    for j = 3:2: ny-2
        ic = 2
        for i = 3:2: nx-2
            coarse[ic, jc] = fine[i, j]
            ic = ic + 1
        end
        jc = jc + 1
    end

    # apply Neumann BCs for temperature T
    if apply_BCs
        apply_boundary_conditions_neumann!(coarse)
    end

end

""" Restricts a field from a fine grid to a coarse grid of half the resolution. """
@parallel_indices (ix, iy) function restrict!(fine, coarse)
    nx, ny = size(fine)
    if (ix % 2 == 1 && iy % 2 == 1) && (3 <= ix <= nx-2 && 3 <= iy <= ny-2)
        ix_c, iy_c = ((ix-1)÷2)+1, ((iy-1)÷2)+1
        coarse[ix_c, iy_c] = fine[ix, iy]
    end
    return nothing
end

""" Restrict wrapper.

Calls the proper restrict function given an execution policy and applies boundary conditions on the coarser level.

"""
@views function restrict_wrapper!(fine, coarse, apply_BCs, execution_policy)
    # initalize coarse to zero (+ Dirichlet(0) BCs)
    coarse[:] .= 0.0

    if execution_policy in [parallel, parallel_shmem]
        @parallel restrict!(fine, coarse)
    else
        restrict_serial!(fine, coarse)
    end

    # apply Neumann BCs for temperature T
    if apply_BCs
        apply_boundary_conditions_neumann!(coarse)
    end
end

""" Prolongates a field from a coarse grid to a fine grid of double the resolution.

This is a serial version, mainly used for debugging.

"""
@views function prolongate_serial!(coarse, fine, apply_BCs)
    a2, a4 = 1.0 / 2.0, 1.0 / 4.0

    nx = size(fine, 1)
    ny = size(fine, 2)

    # initialize fine to zero (+ Dirichlet(0) BCs)
    fine .= 0.0

    for j = 1:ny
        for i = 1: nx
            ic, jc = ((i-1)÷2)+1, ((j-1)÷2)+1
            if (i % 2 == 1 && j % 2 == 1) && (3 <= i <= nx-2 && 3 <= j <= ny-2)
                fine[i, j] = coarse[ic, jc]
                fine[i+1, j] = fine[i+1, j] + a2 * coarse[ic, jc]
                fine[i-1, j] = fine[i-1, j] + a2 * coarse[ic, jc]
                fine[i, j+1] = fine[i, j+1] + a2 * coarse[ic, jc]
                fine[i, j-1] = fine[i, j-1] + a2 * coarse[ic, jc]
                fine[i+1, j+1] = fine[i+1, j+1] + a4 * coarse[ic, jc]
                fine[i+1, j-1] = fine[i+1, j-1] + a4 * coarse[ic, jc]
                fine[i-1, j+1] = fine[i-1, j+1] + a4 * coarse[ic, jc]
                fine[i-1, j-1] = fine[i-1, j-1] + a4 * coarse[ic, jc]
            end
        end
    end

    # apply Neumann BCs for temperature T
    if apply_BCs
        apply_boundary_conditions_neumann!(fine)
    end

end

""" Prolongates a field from a coarse grid to a fine grid of double the resolution.

Uses atomic operation for avoiding issues when run on the GPU.

"""
@parallel_indices (ix, iy) function prolongate_with_atomic!(coarse, fine)
    nx, ny = size(fine)
    a2 = 1.0 / 2.0
    a4 = 1.0 / 4.0
    if (ix % 2 == 1 && iy % 2 == 1) && (3 <= ix <= nx-2 && 3 <= iy <= ny-2)
        ix_c, iy_c = ((ix-1)÷2)+1, ((iy-1)÷2)+1
        @atomic fine[ix, iy] = fine[ix, iy] + coarse[ix_c, iy_c]
        @atomic fine[ix+1, iy] = fine[ix+1, iy] + a2 * coarse[ix_c, iy_c]
        @atomic fine[ix-1, iy] = fine[ix-1, iy] + a2 * coarse[ix_c, iy_c]
        @atomic fine[ix, iy+1] = fine[ix, iy+1] + a2 * coarse[ix_c, iy_c]
        @atomic fine[ix, iy-1] = fine[ix, iy-1] + a2 * coarse[ix_c, iy_c]
        @atomic fine[ix+1, iy+1] = fine[ix+1, iy+1] + a4 * coarse[ix_c, iy_c]
        @atomic fine[ix+1, iy-1] = fine[ix+1, iy-1] + a4 * coarse[ix_c, iy_c]
        @atomic fine[ix-1, iy+1] = fine[ix-1, iy+1] + a4 * coarse[ix_c, iy_c]
        @atomic fine[ix-1, iy-1] = fine[ix-1, iy-1] + a4 * coarse[ix_c, iy_c]
    end
    return nothing
end

""" Prolongates a field from a coarse grid to a fine grid of double the resolution.

This version does not use atomics. We only used for execution on the CPU.

"""
@parallel_indices (ix, iy) function prolongate!(coarse, fine)
    nx, ny = size(fine)
    a2 = 1.0 / 2.0
    a4 = 1.0 / 4.0
    if (ix % 2 == 1 && iy % 2 == 1) && (3 <= ix <= nx-2 && 3 <= iy <= ny-2)
        ix_c, iy_c = ((ix-1)÷2)+1, ((iy-1)÷2)+1
        fine[ix, iy] = fine[ix, iy] + coarse[ix_c, iy_c]
        fine[ix+1, iy] = fine[ix+1, iy] + a2 * coarse[ix_c, iy_c]
        fine[ix-1, iy] = fine[ix-1, iy] + a2 * coarse[ix_c, iy_c]
        fine[ix, iy+1] = fine[ix, iy+1] + a2 * coarse[ix_c, iy_c]
        fine[ix, iy-1] = fine[ix, iy-1] + a2 * coarse[ix_c, iy_c]
        fine[ix+1, iy+1] = fine[ix+1, iy+1] + a4 * coarse[ix_c, iy_c]
        fine[ix+1, iy-1] = fine[ix+1, iy-1] + a4 * coarse[ix_c, iy_c]
        fine[ix-1, iy+1] = fine[ix-1, iy+1] + a4 * coarse[ix_c, iy_c]
        fine[ix-1, iy-1] = fine[ix-1, iy-1] + a4 * coarse[ix_c, iy_c]
    end
    return nothing
end

""" Prolongate wrapper.

Calls the proper prolongate function given an execution policy and applies boundary conditions on the finer level.

"""
@views function prolongate_wrapper!(coarse, fine, apply_BCs, execution_policy)
    # initialize fine to zero (+ Dirichlet(0) BCs)
    fine[:] .= 0.0

    # We need the @atomic macro to avoid race-conditions on the GPU,
    # but this is not supported for CPU :(
    if execution_policy in [parallel, parallel_shmem]
        if USE_GPU
            @parallel prolongate_with_atomic!(coarse, fine)
        else
            @parallel prolongate!(coarse, fine)
        end
    else
        prolongate_serial!(coarse, fine)
    end

    # apply Neumann BCs for temperature T
    if apply_BCs
        apply_boundary_conditions_neumann!(fine)
    end

end


text:
"""
Utils module containing helper functions and macros for, e.g., block
decompositions and visualization.
"""
module Utils

using LazySets

# Visualization
export print_sparsity,
       plot_sparsity,
       print_last_interval,
       range_last_x_percent,
       add_plot_labels

# Block Structure
export @block_id,
       add_dimension,
       block_to_set_map,
       convert_partition,
       compute_block_sizes

# Usability
export @filename_to_png,
       @relpath

# internal conversion
export interpret_template_direction_symbol

# temporary helper function
export decompose_helper

# Extension of MathematicalSystems for use inside Reachability.jl
include("systems.jl")

"""
    @filename_to_png

Expands into the current filename and transforms the file extension to `png`.

### Examples

If this macro is executed from a script named my_model.jl, then:

```julia
julia> plot_name = @filename_to_png
my_model.png
```
"""
macro filename_to_png()
    return :(split(split(@__FILE__, "/")[end], ".")[1] * ".png")
end

"""
    @block_id v

Returns the block number associated to a given variable/dimension.

It is assumed that the block size is two and that blocks are ordered from top
(first two rows) to bottom (last two rows) in a matrix.

### Examples

```julia
julia> @block_id 4
2
```
"""
macro block_id(v::Int64)
    return :(div($v, 2) + mod($v, 2))
end

"""
    range_last_x_percent(N, first=0, step=1)

Returns a StepRange for the last x percent of the range 1:N.

### Input

- `N`     -- range length
- `first` -- relative starting index (in percent); a number from [0, 99]
                (default: 0)
- `step`  -- step size (default: 1)

EXAMPLE:

```julia
julia> range_last_x_percent(200, 20, 5)
160:5:200
```
"""
function range_last_x_percent(N::Int64, first::Int64=0, step::Int64=1)
    if first == 0
        start = 1
    elseif first > 0 && first < 100
        start = (100-first)*(floor(Int, N/100))
        start = max(1, start + ((N-start) % step)) # last index should be N
    else
        throw(DomainError("Start index must be in [0, 99]."))
    end
    return start : step : N
end

"""
    print_last_interval(RsetsProj)

Prints the max/min values for the last element in the array of two-dimensional
sets in the second dimension.

### Input

- `RsetsProj` -- (projected) reach set given as an array of two-dimensional sets
"""
function print_last_interval(RsetsProj::Vector{<:LazySet})
    last_polygon = RsetsProj[end]
    min_val = -ρ([0., -1.], last_polygon)
    max_val = ρ([0., 1.], last_polygon)

    println("max/min values for the last time step:")
    @printf "[ %.6e ; %.6e ]\n" max_val min_val
end

"""
    print_sparsity(ϕ, name)

Prints the sparsity of a matrix ϕ.

### Input

- `ϕ`    -- matrix
- `name` -- (optional, default: "") matrix name (for the output message)
"""
function print_sparsity(ϕ::AbstractMatrix{Float64}, name::String="")
    zero_blocks = 0
    b = (Int64)(size(ϕ, 1) / 2)
    @inline F(bi::Int64, bj::Int64) = ϕ[(2*bi-1):(2*bi), (2*bj-1):(2*bj)]
    for bj in 1:b
        for bi in 1:b
            if findfirst(F(bi, bj)) == 0
                zero_blocks += 1
            end
        end
    end
    all_blocks = b^2

    print_sparsity_message(name, zero_blocks, all_blocks)
end

function print_sparsity(ϕ::SparseMatrixExp{Float64}, name::String="")
    zero_blocks = 0
    b = (Int64)(size(ϕ, 1) / 2)
    @inline F(bi::Int64) = get_rows(ϕ, (2*bi-1):2*bi)
    for bi in 1:b
        block_row_bi = F(bi)
        for bj in 1:b
            if findfirst(block_row_bi[1:2, (2*bj-1):(2*bj)]) == 0
                zero_blocks += 1
            end
        end
    end
    all_blocks = b^2

    print_sparsity_message(name, zero_blocks, all_blocks)
end

@inline function print_sparsity_message(name::String, zero_blocks::Int, all_blocks::Int)
    sparsity = zero_blocks / all_blocks * 100.
    @printf "sparsity %s: %.3f %% (%d/%d zero blocks)\n" name sparsity zero_blocks all_blocks
end

"""
    plot_sparsity(ϕ, name, [plot_backend])

Plots the sparsity of a matrix ϕ.

### Input

- `ϕ`            -- dense or sparse matrix
- `name`         -- the name of the matrix (for the file name)
- `plot_backend` -- (optional, default: `''`): name of the plotting backend;
                      valid values are:

                 -  `'pyplot_savefig'` -- use PyPlot package, save to a file

                 -  `'pyplot_inline'` -- use PyPlot package, showing in external program

                 - `''` -- (empty string), no plotting
"""
function plot_sparsity(ϕ::Union{AbstractMatrix{Float64}, SparseMatrixExp{Float64}}, name::String, plot_backend::String="")
    if isempty(plot_backend)
        # plot nothing
    elseif plot_backend in ["pyplot_inline", "pyplot_savefig"]
        if !isdefined(:PyPlot)
              error("this backend requires that your script loads the PyPlot module")
        end
        eval(Expr(:using, :PyPlot))
        if plot_backend == "pyplot_savefig"
            PyPlot.ioff()  # turn off interactive plotting
            fig = PyPlot.figure()
        end
#       PyPlot.gray()
        if ϕ isa SparseMatrixExp
            explicit_matrix = get_rows(ϕ, 1:size(ϕ, 1))
            PyPlot.spy(explicit_matrix)
        else
            PyPlot.spy(ϕ)
        end
        if plot_backend == "pyplot_savefig"
            PyPlot.savefig(name * "_sparsity.png")
            PyPlot.close(fig)
        end
    else
        warn("plotting backend $plot_backend is not supported for sparsity plots")
    end
end

"""
    add_plot_labels(plot_vars, [project_output], [plot_labels])

Creates the axis labels for plotting.

### Input

- `plot_vars`      -- variable indices for plotting (0 stands for time)
- `project_output` -- (optional, default: false) flag indicating if the y axis
                      is an output function
- `plot_labels`    -- (optional) vector of plot labels (empty strings are ignored)
"""
function add_plot_labels(plot_vars::Vector{Int64}, project_output::Bool=false, plot_labels::Vector{String}=["", ""])
    labels = copy(plot_labels)
    if isempty(labels[1])
        xaxis = plot_vars[1]
        labels[1] = (xaxis == 0) ? "t" : "x" * dec(xaxis)
    end
    if isempty(labels[2])
        yaxis = plot_vars[2]
        labels[2] = (yaxis == 0) ? "t" : (project_output ? "y" : "x" * dec(yaxis))
    end
    return labels
end

"""
   @relpath(name)

Returns the path of the current code file.
This is handy, e.g., when calling a function from anywhere that wants to open a
file relative to its own location.

### Input

- `name` -- file name
"""
macro relpath(name::String)
    return quote
        f = @__FILE__
        if f == nothing
            pathdir = ""
        else
            pathdir = join(split(f, "/")[1:end-1], "/")
        end
        if !isempty(pathdir)
            pathdir = pathdir * "/"
        end
        pathdir * $name
    end
end

"""
    block_to_set_map(dict::Dict{Type{<:LazySet},
                                AbstractVector{<:AbstractVector{Int}}})

Invert a map (set type -> block structure) to a map (block index -> set type).

### Input

- `dict` -- map (set type -> block structure)

### Ouput

Vector mapping a block index to the respective set type.
"""
function block_to_set_map(dict::Dict{Type{<:LazySet},
                                     AbstractVector{<:AbstractVector{Int}}})
    # we assume that blocks do not overlap
    set_type = Vector{Type{<:LazySet}}()
    initial_block_indices = Vector{Int}()
    @inbounds for (key, val) in dict
        for bi in val
            push!(set_type, key)
            push!(initial_block_indices, bi[1])
        end
    end
    s = sortperm(initial_block_indices)
    return set_type[s]
end

"""
    convert_partition(partition::AbstractVector{<:AbstractVector{Int}})::Union{
        Vector{Int}, Vector{UnitRange{Int}}, Vector{Union{UnitRange{Int}, Int}}}

Convert a partition data structure to a more efficient representation where a 1D
block is represented by the index and a higher-dimensional block is represented
by a `UnitRange{Int}`.

### Input

- `partition` -- partition data structure (a vector of integer vectors)

### Output

New partition representation.

### Notes

We cannot mix `Int` and `UnitRange{Int}`, so we only create `Int` partitions if
all blocks are 1D.
"""
function convert_partition(partition::AbstractVector{<:AbstractVector{Int}})::Union{
        Vector{Int},
        Vector{UnitRange{Int}}}
    # are there kD blocks for k > 1?
    has_kD_blocks = false
    for (i, block) in enumerate(partition)
        if length(block) > 1
            has_kD_blocks = true
            break
        end
    end

    # use optimal partition type
    if !has_kD_blocks
        # only 1D blocks
        partition_out = Vector{Int}(length(partition))
        for (i, block) in enumerate(partition)
            partition_out[i] = block[1]
        end
    else
        # at least one kD block for k > 1
        partition_out = Vector{UnitRange{Int}}(length(partition))
        for (i, block) in enumerate(partition)
            partition_out[i] = block[1]:block[end]
        end
    end

    return partition_out
end

"""
    interpret_template_direction_symbol(symbol::Symbol)

Return a template direction type for a given symbol.

### Input

- `symbol` -- symbol

### Output

The template direction type if it is known, or `nothing` otherwise.
"""
function interpret_template_direction_symbol(symbol::Symbol)
    if symbol == :box
        dir = Approximations.BoxDirections
    elseif symbol == :oct
        dir = Approximations.OctDirections
    elseif symbol == :boxdiag
        dir = Approximations.BoxDiagDirections
    else
        if symbol != :nothing
            warn("ignoring unknown template direction $symbol")
        end
        dir = nothing
    end
    return dir
end

"""
    compute_block_sizes(partition::Union{Vector{Int}, Vector{UnitRange{Int}}}
                       )::Vector{Int}

Conversion of the blocks representation from partition to block sizes.

### Input

- `partition` -- partition, a vector of block index ranges

### Output

A vector where the ``i``th entry contains the size of the ``i``th block.
"""
function compute_block_sizes(partition::Union{Vector{Int}, Vector{UnitRange{Int}}}
                            )::Vector{Int}
    res = Vector{Int}(length(partition))
    for (i, block) in enumerate(partition)
        res[i] = length(block)
    end
    return res
end

"""
This is a temporary decomposition function that only applies a projection and
keeps the sets lazy.
Eventually this should be handled in LazySets.
"""
function decompose_helper(S::LazySet{N}, blocks::AbstractVector{Int},
                          n::Int=dim(S)) where {N}
    result = Vector{LazySet{N}}(length(blocks))
    block_start = 1
    @inbounds for (i, bi) in enumerate(blocks)
        M = sparse(1:bi, block_start:(block_start + bi - 1), ones(N, bi), bi, n)
        result[i] = M * S
        block_start += bi
    end
    return CartesianProductArray(result)
end

end # module

ext: jl


text:
using Revise, Distributions, StatsBase, Random, BenchmarkTools, Finometrics

Random.seed!(1111)

mbsampleblockstart(N, runlength) = sample(1:(N-runlength+1))

#samplerunlength(p::Float64) = rand(Geometric(p)) + 1

#draws a non-parametric bootstrap array
function mbbootstraparray!(dat::Vector{Tdat}, sim::Vector{Tsim}; w::Int) where {Tdat, Tsim}
  simrunstart::Int = 1
  N = length(dat[1])

  #need to shift by 1 due to the way the distribution is parametrized
  while simrunstart≤N
    #draw the length of the run from the target distribution
    simrunend = min(simrunstart+w-1,N)
    runlength = simrunend - simrunstart + 1

    #draw a set of bootstrap indices
    datstart = mbsampleblockstart(N,runlength)

    #for each vector, copy the relevant indices
    for (v,b) ∈ zip(dat,sim)
      copyto!(b, simrunstart, v, datstart, runlength)
    end

    simrunstart += runlength

  end

  return nothing
end

#run the aggregation functions on the bootstrap array
function mbbootstrappass!(dat::Vector{Tv}, statfunc::Tstatfunc, sim::Vector{Tv} = deepcopy(dat);
  w::Int,
  ) where {Tv<:AbstractVector, Tstatfunc}

  @assert (length.(dat) .=== length.(sim)) |> all

  mbbootstraparray!(dat, sim; w)

  return statfunc(sim)
end

#this runs the full bootstrap
#individual passes are NOT stored
function mbbootstrap(dat::Vector{Tv}, statfunc::Tstatfunc, aggfunc::Taggfunc;
    B::Int, w::Int) where {Tv, Tstatfunc, Taggfunc}

  #pre-allocate vectors used for the bootstrap
  sims = [deepcopy(dat) for t ∈ 1:Threads.nthreads()]
  results = Dict(s => Vector{MFloat64}(undef, B) for s ∈ keys(statfunc(dat)))

  #draw a sample and compute the aggregate function, sotring the reuslts in an array
  Threads.@threads for b ∈ 1:B
    t = Threads.threadid()
    res = mbbootstrappass!(dat, statfunc, sims[t]; w)

    for s ∈ keys(res)
      @assert results[s][b] === missing
      results[s][b] = res[s]
    end
  end

  return (; Dict(s=>aggfunc(results[s]) for s ∈ keys(results))...)
end

nonparametricp(v::AbstractVector{T}) where T<:Real = sum(v .≤ zero(T))/length(v)
nonparametricp(v::AbstractVector) = (
  sum(v .|> (x)->coalesce(x≤0.0, 0.0))/sum(v .!== missing))

function testmbbootstrap(N, B; K=2)
  dat = [rand(Normal(1,1), N) |> Vector{Union{Missing, Float64}} for i ∈ 1:K]

  #a placeholder aggregation function
  statfunc(sim) = Dict{Symbol, Float64}(
    :doublemean => (sum(sim[1]) + sum(sim[2]))/(length(sim[1]) + length(sim[2])),
    :pear=>cor(sim[1],sim[2]),
    #:kend=>corkendall(sim[1] |> Vector{Float64}, sim[2] |> Vector{Float64}),
    :spear=>corspearman(sim[1] |> Vector{Float64},sim[2] |> Vector{Float64}))

  aggfunc(v) = std(skipmissing(v))
  #aggfunc = nonparametricp

  w = 5

  sim = deepcopy(dat)
  for v ∈ sim
    v .= missing
  end

  mbbootstrappass!(dat, statfunc, sim; w)
  #@code_warntype bootstrappass!(dat, statfunc, sim; p)
  for b ∈ sim
    @assert all(b .!== missing)
  end

  results = mbbootstrap(dat, statfunc, aggfunc; B, w)
  isapprox(results.doublemean, 1/(2N)^0.5, atol=1e-2) || @warn(
    "!isapprox(results.doublemean, 1/(2N)^0.5, atol=1e-2) - could be an issue if this is unexpected")

  #create an autocorrelated series
  vauto = Vector{MFloat64}(undef, N)
  v = rand(Normal(2,2), N)
  for i ∈ 1:N
    vauto[i] = mean(v[max((i-11),1):i])
  end

  for v ∈ dat
    v .+= vauto
  end

  w=7
  nwΣ = neweywestΣfunc(5)
  X::Matrix{Float64} = hcat(dat[1] .* std(dat[2]) ./ std(dat[1]), ones(N))
  lm = FMLM(X, dat[2] |> Vector{Float64})
  pearhacz = nwΣ(lm)[1,1]^0.5

  results = mbbootstrap(dat, statfunc, aggfunc;B, w)
  @info "newey HAC: $pearhacz,
    moving se: $(results.pear),
    delta: $(abs(pearhacz-results.pear))
    cor: $(lm.β[1])
    std: $(((N-2)/(1-lm.β[1]^2))^(-0.5))"

  print("Time per bootstrap:")
  @btime mbbootstrap($dat, $statfunc, $aggfunc; B=$B, w=$w)

  return nothing
end

#@time testmbbootstrap(300,10000)

###################################
#The below is a variant that uses random blocks- shouldn't make much of a difference
#choose the moving block bootstrap since it seems more accepted in literature and is less complex
#draws a single run of an array
function sbsampleblockindices(N::Int, runlength::Int)::Vector{Int}
  datrunstart = sample(1:N)
  local out::Vector{Int}
  if datrunstart + runlength - 1 ≤ N
    out = datrunstart:(datrunstart + runlength - 1) |> collect
  else
    toplength = N - datrunstart + 1
    out = [datrunstart:N; 1:(runlength-toplength);]
  end

  #@assert length(out) == runlength "runlength=$runlength but out=$out"
  return out
end

samplerunlength(p::Float64) = rand(Geometric(p)) + 1

#draws a non-parametric bootstrap array
function sbbootstraparray!(dat::Vector{Tdat}, sim::Vector{Tsim}; p::Float64) where {Tdat, Tsim}
  simrunstart::Int = 1
  N = length(dat[1])

  #need to shift by 1 due to the way the distribution is parametrized
  while simrunstart≤N
    #draw the length of the run from the target distribution
    targetrunlength = samplerunlength(p)
    simrunend = min(simrunstart+targetrunlength-1,N)
    runlength = simrunend - simrunstart + 1

    #draw a set of bootstrap indices
    inds = sbsampleblockindices(N,runlength)

    #for each vector, copy the relevant indices
    for (v,b) ∈ zip(dat,sim)
      copyto!(b, simrunstart, @view v[inds])
    end

    simrunstart += runlength

  end

  return nothing
end

#run the aggregation functions on the bootstrap array
function sbbootstrappass!(dat::Vector{Tv}, statfunc::Tstatfunc, sim::Vector{Tv} = deepcopy(dat);
  p::Float64,
  ) where {Tv<:AbstractVector, Tstatfunc}

  @assert (length.(dat) .=== length.(sim)) |> all

  sbbootstraparray!(dat, sim; p)

  return statfunc(sim)
end

#this runs the full bootstrap
#individual passes are NOT stored
function sbbootstrap(dat::Vector{Tv}, statfunc::Tstatfunc, aggfunc::Taggfunc;
    B::Int, p::Float64) where {Tv, Tstatfunc, Taggfunc}

  #pre-allocate vectors used for the bootstrap
  sims = [deepcopy(dat) for t ∈ 1:Threads.nthreads()]
  results = Dict(s => Vector{MFloat64}(undef, B) for s ∈ keys(statfunc(dat)))

  #draw a sample and compute the aggregate function, sotring the reuslts in an array
  Threads.@threads for b ∈ 1:B
    t = Threads.threadid()
    res = sbbootstrappass!(dat, statfunc, sims[t]; p)

    for s ∈ keys(res)
      @assert results[s][b] === missing
      results[s][b] = res[s]
    end
  end

  return (; Dict(s=>aggfunc(skipmissing(results[s])) for s ∈ keys(results))...)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function getCoefficientsSRK3(T=Float64)
    a = @dec128 [
         [5/36         2/9        5/36-√15/10]
         [5/36         2/9        5/36       ]
         [5/36+√15/10  2/9        5/36       ]
        ]
    b = @dec128 [5/18,        4/9,       5/18       ]
    c = @dec128 [1/2-√15/10,  1/2,       1/2+√15/10 ]
    o = 4

    CoefficientsRK(T, :SRK3, o, a, b, c)
end


text:
using Weave
using Base.Test

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="tex")
result = readall(open("documents/pyplot_formats.tex"))
ref = readall(open("documents/pyplot_formats_ref.tex"))
result = replace(result, r"\s*PyObject.*\n", "") #Remove PyObjects, because they change
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="github")
result = readall(open("documents/pyplot_formats.md"))
ref = readall(open("documents/pyplot_formats_ref.md"))
result = replace(result, r"\s*PyObject.*\n", "")
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="rst", fig_ext=".svg")
result = readall(open("documents/pyplot_formats.tex"))
ref = readall(open("documents/pyplot_formats_ref.tex"))
result = replace(result, r"\s*PyObject.*\n", "")
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref


text:
"""
    UndefArray(size [, factory])
    UndefArray{T}(size [, factory])
    UndefArray{T,N}(size [, factory])

# Examples
```jldoctest
julia> using MicroCollections

julia> UndefArray((2,))
2-element UndefVector{Union{}}(2):
 #undef
 #undef

julia> UndefArray{Int}((2, 3))
2×3 UndefArray{2,Int64}((2, 3)):
 #undef  #undef  #undef
 #undef  #undef  #undef
```

The size of an `UndefArray` can be "changed" by using
[`Setfield.@set`](https://github.com/jw3126/Setfield.jl)

```jldoctest; setup = :(using MicroCollections)
julia> using Setfield

julia> x = UndefArray((2,))
2-element UndefVector{Union{}}(2):
 #undef
 #undef

julia> @set size(x) = (1, 3)
1×3 UndefArray{2,Union{}}((1, 3)):
 #undef  #undef  #undef
```
"""
UndefArray

"""
    UndefVector(length)

# Examples
```jldoctest
julia> using MicroCollections

julia> UndefVector(3)
3-element UndefVector{Union{}}(3):
 #undef
 #undef
 #undef
```
"""
UndefVector

struct UndefArray{T,N,F} <: AbstractArray{T,N}
    size::NTuple{N,Int}
    factory::F
end

const UndefVector{T,F} = UndefArray{T,1,F}

default_factory(T, size) = Array{T}(undef, size)

UndefArray{T,N}(size::NTuple{N,Int}, factory::F = default_factory) where {T,N,F} =
    UndefArray{T,N,F}(size, factory)
UndefArray{T}(size::Tuple{Vararg{Int}}, factory = default_factory) where {T} =
    UndefArray{T,length(size)}(size, factory)
UndefArray(size::Tuple{Vararg{Int}}, factory = default_factory) =
    UndefArray{Union{}}(size, factory)
UndefArray(size::Int...) = UndefArray(size)

UndefVector{T}(length::Int, factory = default_factory) where {T} =
    UndefArray{T,1}((length,), factory)
UndefVector(length::Int, factory = default_factory) = UndefArray((length,), factory)

UndefArray{T,N}(::Tuple{}, factory::F = default_factory) where {T,N,F} =
    UndefArray{T,N,F}(ntuple(_ -> 0, N), factory)

Base.size(a::UndefArray) = a.size
Base.isassigned(::UndefArray, ::Integer...) = false

Base.similar(a::UndefArray, T::Type, size::Tuple{Vararg{Int}}) = a.factory(T, size)

function Base.showarg(io::IO, x::UndefArray, _toplevel::Bool)
    if ndims(x) == 1
        print(io, "UndefVector{", eltype(x), "}")
    else
        print(io, "UndefArray{", ndims(x), ',', eltype(x), '}')
    end
    print(io, '(')
    if length(x.size) == 1
        show(io, x.size[1])
    else
        show(io, x.size)
    end
    if x.factory !== default_factory
        print(io, ", ")
        show(io, x.factory)
    end
    print(io, ')')
end

function Setfield.set(x::UndefArray, ::typeof(@lens size(_)), dims::Tuple{Vararg{Integer}})
    dims = convert(Tuple{Vararg{Int}}, dims)
    return UndefArray{eltype(x)}(dims, x.factory)
end


text:
include("../roblib.jl")

function mark_angels(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))
    putmarker!(r)
    through_rectangles_into_angle(r, (Up, Left))
    putmarker!(r)
    through_rectangles_into_angle(r, (Up, Right))
    putmarker!(r)
    through_rectangles_into_angle(r, (Down, Right))
    putmarker!(r)
    through_rectangles_into_angle(r, (Down, Left))
    movements_if_possible!(r, (Up, Right), reverse(num_steps))
end


text:
using LazyArrays, FillArrays, ArrayLayouts, Test
import LazyArrays: materialize, broadcasted, DefaultApplyStyle, Applied, arguments,
            ApplyArray, ApplyMatrix, ApplyVector, LazyArrayApplyStyle, ApplyLayout, call
import ArrayLayouts: StridedLayout

@testset "Applying" begin
    @testset "Applied" begin
        @test applied(exp,1) isa Applied{DefaultApplyStyle}
        @test apply(randn) isa Float64
        @test materialize(applied(*, 1)) == apply(*,1) == 1
        @test apply(exp, 1) === exp(1)
        @test apply(exp, broadcasted(+, 1, 2)) === apply(exp, applied(+, 1, 2)) === exp(3)
    end

    @testset "ApplyArray" begin
        A = randn(2,2)
        M = ApplyMatrix(exp, A)
        @test eltype(M) == eltype(Applied(M)) == Float64
        @test M == exp(A)
        @test ndims(M) == 2
        @test axes(M) == (Base.OneTo(2), Base.OneTo(2))

        b = randn(2)
        c = ApplyVector(*, ApplyMatrix(exp, A), b)

        @test axes(c) == (Base.OneTo(2),)

        @test c[1] == c[1,1]
        @test exp(A)*b ≈ c

        @test ApplyArray(+,[1,2],[3,4]) == ApplyVector(+,[1,2],[3,4]) == ApplyArray(+,[1,2],[3,4])

        M = ApplyArray(*, Ones(100_000_000,100_000_000), Diagonal(1:100_000_000))
        @test M[1,1] === 1.0
        @test M[1:10,1:10] == ones(10,10)*Diagonal(1:10)

        M = copy(Applied{LazyArrayApplyStyle}(exp, (A,)))
        @test M isa ApplyArray
        @test M == exp(A)
    end

    @testset "copy (#85)" begin
        Base.size(A::Applied) = (length(A.args[1]),)
        Base.eltype(A::Applied)  = eltype(A.args[1])
        v = ApplyVector(vec, ones(Int, 2, 2))
        vc = copy(float.(v))
        ve = convert(Vector, vc)
        @test eltype(ve) == Float64

        A = ApplyArray(exp, randn(5,5))
        @test copy(A) ≡ map(copy,A) ≡ A
    end

    @testset "copyto!" begin
        a = ApplyArray(+,[1,2],[3,4])
        b = ApplyArray(+,[3,4],[5,6])
        c = ApplyArray(+,[3.0,4.0],[5.0,6.0])
        @test copyto!(b, a) == a == b
        @test copyto!(c, a) == a == c
        @test copyto!(similar(a), a) == a
        if VERSION < v"1.8.0-"
            @test_throws ErrorException copyto!(c, Array(a))
        else
            @test_throws Base.CanonicalIndexError copyto!(c, Array(a))
        end
    end

    @testset "view" begin
        a = ApplyArray(+,[1,2,3],[3,4,5])
        v = view(a,1:2)
        @test MemoryLayout(typeof(v)) isa ApplyLayout{typeof(+)}
        @test call(v) == call(a) == +
        @test Array(v) == a[1:2] == Array(a)[1:2]
    end

    @testset "rot180" begin
        A = randn(3,2)
        R = ApplyArray(rot180, A)
        @test eltype(R) == Float64
        @test size(R) == size(A)
        @test MemoryLayout(R) isa StridedLayout
        @test strides(R) == (-1,-3)
        @test pointer(R) == pointer(view(A,3,2))
        @test R == rot180(A) == copyto!(similar(A),R)

        R = ApplyArray(rotl90, A)
        @test eltype(R) == Float64
        @test size(R) == (2,3)
        @test R == rotl90(A)

        R = ApplyArray(rotr90, A)
        @test eltype(R) == Float64
        @test size(R) == (2,3)
        @test R == rotr90(A)

        B = randn(2,3)
        R = ApplyArray(rot180, ApplyArray(*, A, B))
        @test MemoryLayout(R) isa ApplyLayout{typeof(*)}
        @test arguments(R) == (rot180(A), rot180(B))
        @test R ≈ rot180(A*B)
    end
end


text:
function MGCG!(A,b,x,H1,exact, my_solver;smooth_steps = 4,maxiter=length(b)^2,abstol=sqrt(eps(real(eltype(b)))))
    r = b - A * x;
    z = zeros(length(b))
    (M, R, H, I_p, A_2h, I_r, IN) = precond_matrix(A,b;m=smooth_steps,solver="jacobi")
    # Two_level_multigrid!(A,r,z;nu=smooth_steps, solver = my_solver)
    z = M*r #test
    p = z;
    rzold = r'*z

    num_iter_steps = 0
    norms = [sqrt(r'*H1*r)]

    diff = x-exact
    err = sqrt(diff'*A*diff)

    E = [err]
    for step = 1:maxiter
    # for step = 1:5
        Ap = A*p;
        @show norm(A*p)
        num_iter_steps += 1
        alpha = rzold/(p'*Ap)
        @show alpha
        x .= x .+ alpha * p;

        diff = x-exact
        err = sqrt(diff'*A*diff)
        @show err
        append!(E, err)

        r .= r .- alpha * Ap;

        norm_r = sqrt(r' * r)
        append!(norms,norm_r)

        if norm_r < abstol
            break
        end

        z .= 0
        # Two_level_multigrid!(A,r,z;nu=smooth_steps, solver = my_solver)
        z = M*r

        rznew = r' * z
        beta = rznew/rzold;
        @show beta
        p = z + beta * p;
        rzold = rznew
    end

    return E, num_iter_steps, norms
end

using Random
Random.seed!(777)
using LinearAlgebra
using SparseArrays
using Plots

function e(i,n)
    out = spzeros(n)
    out[i] = 1.0
    return out
end

function eyes(n)
    out = spzeros(n,n)
    for i in 1:n
        out[i,i] = 1.0
    end
    return out
end

function u(x, y)
    return sin.(π*x) * (cosh.(π*y))'
end

function u_x(x, y)
    return π * cos.(π*x) * (cosh.(π*y))'
end

function u_y(x, y)
    return π * sin.(π*x) * (sinh.(π*y))'
end

function u_xx(x, y)
    return -π^2 * sin.(π*x) * (cosh.(π*y))'
end

function u_yy(x, y)
    return π^2 * sin.(π*x) * (cosh.(π*y))'
end
function Diag(A)
    # Self defined function that is similar to Matlab Diag
    return Diagonal(A[:])
end

function get_operators(p, h)
    N = Integer(1/h)
    Identity = eyes(N+1)

    # Create some random positive definite coefficients

    x = range(0,step=h,1);
    y = range(0,step=h,1);
    B = ones(N+1)

    (D1, HI, H1, r1) = diagonal_sbp_D1(p,N,xc=(0,1));
    (D2, BS, HI2, H2, r2) = variable_diagonal_sbp_D2(p,N,B;xc=(0,1));

    # Forming 2d Operators
    Hx = kron(Identity, H1)
    HxI = kron(Identity, HI)
    BSx = kron(Identity, BS)

    Hy = kron(H1, Identity)
    HyI = kron(HI, Identity)
    BSy = kron(BS, Identity)

    HH = kron(H1, H1)

    D2xplusD2y = kron(Identity, D2) + kron(D2, Identity)
    e1 = e(1,N+1);
    eN = e(N+1,N+1);

    E_N = sparse(kron(Diag(eN), Identity));
    E_S = sparse(kron(Diag(e1), Identity));
    E_E = sparse(kron(Identity, Diag(eN)));
    E_W = sparse(kron(Identity, Diag(e1)));

    # Penalty Parameters
    tau_E = -13/h;
    tau_W = -13/h;
    tau_N = -1;
    tau_S = -1;
    beta = 1;
    SAT_W = tau_W*HxI*E_W + beta*HxI*BSx'*E_W;
    SAT_E = tau_E*HxI*E_E + beta*HxI*BSx'*E_E;
    #SAT_W_r = tau_W*crr[1]*HyI*E_W*e1 + beta*HyI*BSy'*E_W*e1;
    #SAT_E_r = tau_E*crr[end]*HyI*E_E*eN + beta*HyI*BSy'*E_E*eN;

    g_W = 0;
    g_E = 0;
    g_N = u_y(x, 1);
    g_S = u_y(x, 0);

    SAT_N = tau_N*HyI*E_N*BSy;
    SAT_S = tau_S*HyI*E_S*BSy;
    SAT_N_r = tau_N*HyI*kron(eN, Identity);
    SAT_S_r = tau_S*HyI*kron(e1,Identity);

    A = D2xplusD2y + SAT_W + SAT_E + SAT_N + SAT_S;

    b =  0 .* u(x, y')[:] + SAT_N_r*g_N  + SAT_S_r*g_S;

    A = -HH*A;
    b = Vector(-HH*b);

    return (A, b, x, y, HH)
end

function prolongation_matrix(N)
    # SBP preserving
    odata = spzeros(2*N-1,N)
    for i in 1:2*N-1
        if i % 2 == 1
            odata[i,div(i+1,2)] = 1
        else
            odata[i,div(i,2)] = 1/2
            odata[i,div(i,2)+1] = 1/2
        end
    end
    return odata
end

function restriction_matrix(N)
    # SBP preserving
    odata = spzeros(div(N+1,2),N)
    odata[1,1] = 1/2
    odata[1,2] = 1/2
    odata[end,end-1] = 1/2
    odata[end,end] = 1/2
    for i in 2:div(N+1,2)-1
        odata[i,2*i-2] = 1/4
        odata[i,2*i-1] = 1/2
        odata[i,2*i] = 1/4
    end
    return odata
end

function standard_prolongation_matrix(N)
    # SBP preserving
    m = (N-1)*2 + 1
    I_r = standard_restriction_matrix(m)
    odata = 2 .* I_r'
    return odata
end

function standard_restriction_matrix(N)
    # SBP preserving
    odata = spzeros(div(N+1,2),N)
    odata[1,1] = 1/2
    odata[1,2] = 1/4
    odata[end,end-1] = 1/4
    odata[end,end] = 1/2
    for i in 2:div(N+1,2)-1
        odata[i,2*i-2] = 1/4
        odata[i,2*i-1] = 1/2
        odata[i,2*i] = 1/4
    end
    return odata
end

function standard_prolongation_matrix_2D(N)

    odata = standard_prolongation_matrix(Integer(sqrt(N)))
    return kron(odata, odata)
end

function standard_restriction_matrix_2D(N)
    odata = standard_restriction_matrix(Integer(sqrt(N)))
    return kron(odata, odata)
end

function ssor!(x,A,b;no_iter=4, ω = 1.4)

    N = length(b)
    IN = sparse(Matrix(I, N, N))
    P = Diagonal(diag(A))
    L = A - triu(A)
    U = A - tril(A)

    B1 = (P + ω*U)\Matrix(-ω*L + (1-ω)*P)
    B2 = (P + ω*L)\Matrix(-ω*U + (1-ω)*P)
    X = (P+ω*L)\Matrix(IN)
    C = ω*(2-ω)*(P+ω*U)\Matrix(P*X)

    for j in 1:no_iter
        x[:] = B1*B2*x[:] + C*b
    end

end

function richardson!(x,A,b;no_iter=4, ω =ω_richardson)

    for j in 1:no_iter
        x[:] = x[:] + ω * (b - A*x[:]) # standard richardson
    end

end

function richardson_chebyshev!(x,A,b;no_iter=4, ω =ω_richardson)

    M = I - ω*A
    evals = eigen(Matrix(M))
    es1 =  extrema(evals.values)
    aa = es1[1]
    bb = es1[2]
    γ = 1-bb
    Γ = 1-aa
    κ = Γ/γ
    c = (sqrt(κ)-1)/(sqrt(κ)+1)
    θ = 4*c/(bb-aa)
    ν = -2*c*(aa+bb)/(bb-aa)

    y0 = copy(x)
    y1 = y0 + ω * (b - A*y0)

    x[:] = y1[:]

   for j in 2:no_iter
        x[:] = θ * (y1[:] + ω * (b - A * y1[:]) - y0) + ν * (y1 - y0) + y0
        y0[:] = y1[:]
        y1[:] = x[:]
    end

end

function chebyshev!(x,A,b;no_iter=4)

    p = zeros(size(b))
    α = 0*rand(1)

    for j in 1:no_iter
        r = b - A*x
        if j == 0
            α[1] = 1/d
        elseif j == 1
            α[1] = 2*d/(2*d^2 - c^2)
        else
            α[1] = 1/(d - α[1]*c^2/4)
        end
        β = α[1]*d - 1
        p = α[1]*r + β*p
        x[:] = x[:] + p[:]

    end

end

function jacobi!(x,A,b;no_iter=4, ω = 2/3)

    Pinv = Diagonal(1 ./ diag(A))
    P = Diagonal(diag(A))
    Q = A-P

    for j in 1:no_iter
        x[:] = ω * Pinv*(b .- Q*x[:]) + (1 - ω)*x[:]
    end

end

function Two_level_multigrid!(A,b,x0;nu=3, solver = "jacobi")
    N = length(b)

    #pre-smoothing with initial guess x0
    if solver == "jacobi"
        jacobi!(x0,A,b;no_iter=nu)
    elseif solver == "ssor"
        ssor!(x0,A,b;no_iter=nu)
    elseif solver == "richardson"
        richardson!(x0,A,b;no_iter=nu)
    elseif solver == "richardson_chebyshev"
        richardson_chebyshev!(x0,A,b;no_iter=nu)
    elseif solver == "chebyshev"
        chebyshev!(x0,A,b;no_iter=nu)
    else

        #add other solvers
    end

    #compute residual on fine grid and restrict
    rh = b - A*x0
    r2h = standard_restriction_matrix(N)*rh

    #Get operators on coarse grid
    (A_2h, b_2h, x_2h, H1_2h) = get_operators(p, 2*h);

    #Direct solve on coarse grid and prolong error
    e_2h = A_2h\r2h
    e_h = standard_prolongation_matrix(length(b_2h))*e_2h

    #Fine grid correction
    x0[:] .+= e_h

    #post-smoothing
    if solver == "jacobi"
        jacobi!(x0,A,b;no_iter=nu)
    elseif solver == "ssor"
        ssor!(x0,A,b;no_iter=nu)
    elseif solver == "richardson"
        richardson!(x0,A,b;no_iter=nu)
    elseif solver == "richardson_chebyshev"
        richardson_chebyshev!(x0,A,b;no_iter=nu)
    elseif solver == "chebyshev"
        chebyshev!(x0,A,b;no_iter=nu)

    else
        #add other solvers
    end

end


text:
module Discord

using Dates
using Distributed
using HTTP
using InteractiveUtils
using JSON
using OpenTrick
using Setfield
using TimeToLive

const DISCORD_JL_VERSION = v"0.1.0"
const API_VERSION = 6
const DISCORD_API = "https://discordapp.com/api"

# Shortcuts for common unions.
const Optional{T} = Union{T, Missing}
const Nullable{T} = Union{T, Nothing}
const OptionalNullable{T} = Union{T, Missing, Nothing}
const StringOrChar = Union{AbstractString, AbstractChar}

# Constant functions.
donothing(args...; kwargs...) = nothing
alwaystrue(args...; kwargs...) = true
alwaysfalse(args...; kwargs...) = false

# Run a function with an acquired semaphore.
function withsem(f::Function, s::Base.Semaphore)
    Base.acquire(s)
    try f() finally Base.release(s) end
end

# Precompile all methods of a function, running it if force is set.
function compile(f::Function, force::Bool; kwargs...)
    for m in methods(f)
        types = Tuple(m.sig.types[2:end])
        precompile(f, types)
        force && try f(mock.(types; kwargs...)...) catch end
    end
end

include(joinpath("types", "types.jl"))
include(joinpath("gateway", "events", "events.jl"))
include(joinpath("client", "client.jl"))
include(joinpath("gateway", "gateway.jl"))
include(joinpath("rest", "rest.jl"))
include(joinpath("utils", "helpers.jl"))
include(joinpath("utils", "commands.jl"))
include("Defaults.jl")

end


text:
################################################################################
#  Copyright 2021, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

@testset "Markov Process" begin

    @testset "Coal Fired Power Generating Unit" begin
        # Example from: A Multi-State Markov Model for a Short-Term Reliability
        # Analysis of a Power Generating Unit (p. 4)
        stdᶜᶠ = include(joinpath(_MSS.BASE_DIR,"test/elements/coal_fired_power_generating_unit.jl"))

        solve!(stdᶜᶠ, MarkovProcess(), tsim = 1.0u"yr")
        prb_m = [np[end] for np in _MSS.get_sprop(stdᶜᶠ, :prob)]

        solve!(stdᶜᶠ, SteadyStateProcess())
        prb_s = [np[end] for np in _MSS.get_sprop(stdᶜᶠ, :prob)]

        @test isapprox(sum(prb_m), 1.0, rtol=1e-6)
        @test all(isapprox.(prb_m, prb_s, rtol=1e-6))
    end

end


text:
"""
    translate!(camera::CameraModel, 𝐭::AbstractArray)

Translate the camera optical center by a vector 𝐭.
"""
function translate!(camera::CameraModel, 𝐭::AbstractArray)
    camera.𝐜 = camera.𝐜 + 𝐭
    camera
end

"""
    relocate!(camera::CameraModel, 𝐭::AbstractArray)

Rotates a camera around its optical center by a rotation matrix 𝐑 and
then translates the optical center by a vector 𝐭.
"""
function relocate!(camera::CameraModel, 𝐑::AbstractArray, 𝐭::AbstractArray)
    camera.𝐞₁ = 𝐑*camera.𝐞₁
    camera.𝐞₂ = 𝐑*camera.𝐞₂
    camera.𝐞₃ = 𝐑*camera.𝐞₃
    camera.𝐜 = camera.𝐜 + 𝐭
    camera
end


text:
void merge(int[] a, int low, int mid, int high)
{
    int[] b = new int[1000];;
    int i = low, j = mid + 1, k = 0;

    while (i <= mid && j <= high) {
        if (a[i] <= a[j]) {
            b[k] = a[i];
            k++; i++;
        }
        else {
            b[k] = a[j];
            k++; j++;
        }
    }
    while (i <= mid) {
        b[k] = a[i];
        k++; i++;
    }

    while (j <= high) {
        b[k] = a[j];
        k++; j++;
    }

    k--;
    while (k >= 0) {
        a[low + k] = b[k];
        k--;
    }
}

void mergesort(int[] a, int low, int high)
{
    if (low < high) {
        int m = (high + low)/2;
        mergesort(a, low, m);
        mergesort(a, m + 1, high);
        merge(a, low, m, high);
    }
}

int main() {
  int[] array = new int[1000];
  int i = 0;
  while (i<array.length) {
    array[i] = 1000-i;
    i++;
  }
  mergesort(array, 0, 999);
  for (int x : array) {
    printInt(x);
  }
return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
using NeuralOperators
using CUDA
using Flux
using GeometricFlux
using Graphs
using Zygote
using Test

CUDA.allowscalar(false)

cuda_tests = [
    "cuda.jl",
]

tests = [
    "Transform/Transform.jl",
    "operator_kernel.jl",
    "loss.jl",
    "model.jl",
    "deeponet.jl",
]

if CUDA.functional()
    append!(tests, cuda_tests)
else
    @warn "CUDA unavailable, not testing GPU support"
end

@testset "NeuralOperators.jl" begin
    for t in tests
        include(t)
    end
end

#=
    　 ／l、    meow?
    ﾞ（ﾟ､ ｡ ７
    　l、ﾞ ~ヽ
    　じしf_, )ノ
=#


text:
using DCEMRI
using Test

ccc4, cccnoisy4 = validate(4, makeplots=false)
@testset "QIBA v4 Extended Tofts Phantom" begin
  @test ccc4[1] >= 0.9998
  @test ccc4[2] >= 0.8904
  @test ccc4[3] >= 0.9999
  @test round(cccnoisy4[1], sigdigits=2) >= 0.97
  @test round(cccnoisy4[2], sigdigits=2) >= 0.70
  @test round(cccnoisy4[3], sigdigits=2) >= 0.97
end

ccc6, cccnoisy6 = validate(6, makeplots=false)
@testset "QIBA v6 Standard Tofts Phantom" begin
  @test ccc6[1] >= 0.9999
  @test ccc6[2] >= 0.9999
  @test round(cccnoisy6[1], sigdigits=2) >= 0.84
  @test round(cccnoisy6[2], sigdigits=2) >= 0.86
end


text:
"""
    ChirpSignal
Chirp Swept-frequency sine generator.
"""
module ChirpSignal

export chirp

"""
    chirp(T, fs, fl, fh; method="linear", phase=0.0) -> Array{Float64,1}

Chirp Swept-frequency sine generator.
# Arguments
- `T`     : signal's length in second
- `fs`    : sampling rate
- `fl`    : lower frequency at the ending time step, fl <= fs
- `fh `   : higher frequency at the ending time step, fh <= fs
- `phase` : phase of the chirp sigal, like sin(wt + phase)
- `method`: available methods are 'linear','quadratic','exponential' and 'logarithmic'; the default is 'linear'
"""
function chirp(T, fs, fl, fh; method="linear", phase=0.0)
    method == "linear" && return chirpLinear(T, fs, fl, fh; phase=phase)
    method == "quadratic" && return chirpQuadratic(T, fs, fl, fh; phase=phase)
    method == "logarithmic" && return chirpLogarithmic(T, fs, fl, fh; phase=phase)
    method == "exponential" && return chirpExponential(T, fs, fl, fh; phase=phase)
end

function chirpLinear(T, fs, fl, fh; phase=0.0)
    Δt = 1.0/fs;
    fl = min(fs/2,max(0,fl));
    fh = min(fs/2,max(0,fh));
    k  = (fh - fl)/T/2;
    t  = 0:Δt:T-Δt;
    y  = @. sin(6.2831853*(k*t + fl)*t + phase);
end

function chirpQuadratic(T, fs, fl, fh; phase=0.0)
    # f(t) = k * t^2 + fl, k = (fh-fl)/T^2
    # F(t) = k/3 * t^3 + fl*t
    Δt = 1.0/fs;
    fl = min(fs/2,max(0,fl));
    fh = min(fs/2,max(0,fh));
    t  = 0:Δt:T-Δt;
    y  = @. sin(6.2831853*( (fh-fl)/(3*T^2) * t^3 + fl * t) + phase)
end

function chirpLogarithmic(T, fs, fl, fh; phase=0.0)
    # f(t) = k*log(t+1) + fl, k = (fh-fl)/log(T+1)
    # F(t) = k*(t+1)*log(t+1) + (fl-k)*t
    (fl <= fh) || error("fl <= fh not met")
    Δt = 1.0/fs;
    fl = min(fs/2,max(0,fl));
    fh = min(fs/2,max(0,fh));
    t  = 0:Δt:T-Δt;
    k  = (fh-fl)/log(T+1);
    y  = @. sin(6.2831853*( k*(t+1)*log(t+1) + (fl-k)*t ) + phase);
end

function chirpExponential(T, fs, fl, fh; phase=0.0)
    # f(t) = exp(k * t) + (fl-1), k = 1/T*log(fh-fl+1)
    # F(t) = 1/k * exp(k * t) + (fl-1)*t
    (fl <= fh) || error("fl <= fh not met")
    Δt = 1.0/fs;
    fl = min(fs/2,max(0,fl));
    fh = min(fs/2,max(0,fh));
    t  = 0:Δt:T-Δt;
    k  = 1/T*log(fh-fl+1);
    y  = @. sin(6.2831853*( (1/k) * exp(k*t) + (fl-1)*t ) + phase);
end

"""
    chirp(T, fs, f::Function; phase=0.0) -> Array{Float64,1}

Function customizable chirp swept-frequency sine generator.
# Arguments
- `T`     : signal's length in second
- `fs`    : sampling rate
- `f`     : function that defines how the frequency changes vs time
- `phase` : phase of the chirp sigal, like sin(wt + phase)
"""
function chirp(T, fs, f::Function; phase=0.0)
    Δt = 1.0/fs;
    τ  = 0:Δt:T-Δt;
    ϕ  = zeros(length(τ),1);
    x  = zeros(length(τ),1);
    ϕ[1] = f(0)*Δt;
    x[1] = sin(6.2831853*ϕ[1] + phase);
    for t = 2:length(τ)
        ϕ[t] = ϕ[t-1] + f(τ[t-1])*Δt;
        x[t] = sin(6.2831853*ϕ[t] + phase);
    end
    return x
end

end  # module

max_stars_repo_name: sonosole/ChirpSignal

max_issues_repo_name: sonosole/ChirpSignal

max_forks_repo_name: sonosole/ChirpSignal


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

"""
    SimpleKriging(γ, μ)
    SimpleKriging(data, γ, μ)

Simple Kriging with variogram model `γ` and constant mean `μ`.

Optionally, pass the geospatial `data` to the [`fit`](@ref) function.

### Notes

* Simple Kriging requires stationary variograms
"""
struct SimpleKriging{G<:Variogram,V} <: KrigingEstimator
  # input fields
  γ::G
  μ::V

  function SimpleKriging{G,V}(γ, μ) where {G<:Variogram,V}
    @assert isstationary(γ) "Simple Kriging requires stationary variogram"
    new(γ, μ)
  end
end

SimpleKriging(γ, μ) = SimpleKriging{typeof(γ),typeof(μ)}(γ, μ)

SimpleKriging(data, γ, μ) = GeoStatsBase.fit(SimpleKriging(γ, μ), data)

nconstraints(::SimpleKriging) = 0

set_constraints_lhs!(::SimpleKriging, LHS::AbstractMatrix, domain) = nothing

factorize(::SimpleKriging, LHS::AbstractMatrix) = cholesky(Symmetric(LHS), check=false)

set_constraints_rhs!(::FittedKriging{<:SimpleKriging}, pₒ) = nothing

function combine(fitted::FittedKriging{<:SimpleKriging},
                 weights::KrigingWeights, z::AbstractVector)
  γ = fitted.estimator.γ
  μ = fitted.estimator.μ
  b = fitted.state.RHS
  λ = weights.λ
  y = z .- μ

  μ + y⋅λ, sill(γ) - b⋅λ
end


text:
################################################################################
# Displacements
################################################################################
function displacement(joint::Translational,
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion;
    rotate::Bool = true)

    vertices = joint.vertices
    d = xb + vector_rotate(vertices[2], qb) - (xa + vector_rotate(vertices[1], qa))
    rotate && (return vector_rotate(d, inv(qa))) : (return d)
end

function displacement_jacobian_configuration(relative::Symbol, joint::Translational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion;
    attjac=true) where T

    vertices = joint.vertices

    if relative == :parent
        d = xb + vector_rotate(vertices[2], qb) - (xa + vector_rotate(vertices[1], qa)) # in the world frame
        X = -rotation_matrix(inv(qa))
        Q = -rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qa, vertices[1])
        Q += ∂rotation_matrix_inv∂q(qa, d)
        attjac && (Q *= LVᵀmat(qa))
    elseif relative == :child
        X = rotation_matrix(inv(qa))
        Q = rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qb, vertices[2])
        attjac && (Q *= LVᵀmat(qb))
    end

    return X, Q
end

################################################################################
# Coordinates
################################################################################
function minimal_coordinates(joint::Translational, xa::AbstractVector, qa::Quaternion, xb::AbstractVector, qb::Quaternion)
    return nullspace_mask(joint) * displacement(joint, xa, qa, xb, qb)
end

function minimal_coordinates_jacobian_configuration(relative::Symbol, joint::Translational,
        xa::AbstractVector, qa::Quaternion, xb::AbstractVector, qb::Quaternion; attjac::Bool=true)
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=attjac)
	return nullspace_mask(joint) * [X Q]
end

################################################################################
# Coordinates
################################################################################
function set_minimal_coordinates!(joint::Translational,
    pnode::Node, cnode::Node,
    timestep;
    Δx::AbstractVector=szeros(input_dimension(joint)))

    pa = joint.vertices[1]
    pb = joint.vertices[2]

    qa = pnode.state.q2
    xa = pnode.state.x2

    qb = cnode.state.q2

    Aᵀ = zerodimstaticadjoint(nullspace_mask(joint))
    xb = xa + vector_rotate(pa + Aᵀ * Δx, qa) - vector_rotate(pb, qb)
    set_maximal_configurations!(cnode; x=xb, q=cnode.state.q2)
    set_previous_configuration!(cnode, timestep)
    return nothing
end

################################################################################
# Velocities
################################################################################
function minimal_velocities(joint::Translational,
		xa::AbstractVector, va::AbstractVector, qa::Quaternion, ϕa::AbstractVector,
		xb::AbstractVector, vb::AbstractVector, qb::Quaternion, ϕb::AbstractVector,
		timestep)
	A = nullspace_mask(joint)

	# 1 step backward in time
	xa1 = next_position(xa, -va, timestep)
	qa1 = next_orientation(qa, -ϕa, timestep)
	xb1 = next_position(xb, -vb, timestep)
	qb1 = next_orientation(qb, -ϕb, timestep)

	# Coordinates
	Δx = A * displacement(joint, xa, qa, xb, qb)
	# Previous step coordinates
	Δx1 = A * displacement(joint, xa1, qa1, xb1, qb1)

	# Finite difference
	Δv = (Δx - Δx1) / timestep
	return Δv
end

function minimal_velocities_jacobian_configuration(relative::Symbol, joint::Translational,
    xa::AbstractVector, va::AbstractVector, qa::Quaternion, ϕa::AbstractVector,
    xb::AbstractVector, vb::AbstractVector, qb::Quaternion, ϕb::AbstractVector,
    timestep)
    A = nullspace_mask(joint)

    # 1 step backward in time
    xa1 = next_position(xa, -va, timestep)
    qa1 = next_orientation(qa, -ϕa, timestep)
    xb1 = next_position(xb, -vb, timestep)
    qb1 = next_orientation(qb, -ϕb, timestep)

    # Coordinates
    Δx = A * displacement(joint, xa, qa, xb, qb)
    # Previous step coordinates
    Δx1 = A * displacement(joint, xa1, qa1, xb1, qb1)

    # Finite difference
    # Δv = (Δx - Δx1) / timestep

    if relative == :parent
        X, Q = displacement_jacobian_configuration(:parent, joint, xa, qa, xb, qb, attjac=false)
        X1, Q1 = displacement_jacobian_configuration(:parent, joint, xa1, qa1, xb1, qb1, attjac=false)
        X1 *= -1.0
        Q1 *= -1.0 * rotational_integrator_jacobian_orientation(qa, -ϕa, timestep, attjac=false)
        Q *= LVᵀmat(qa)
        Q1 *= LVᵀmat(qa)
        J = 1.0 / timestep * A * [X Q]
        J += 1.0 / timestep * A * [X1 Q1]
    elseif relative == :child
        1.0 / timestep * (Δx - Δx1)
        X, Q = displacement_jacobian_configuration(:child, joint, xa, qa, xb, qb, attjac=false)
        X1, Q1 = displacement_jacobian_configuration(:child, joint, xa1, qa1, xb1, qb1, attjac=false)
        X1 *= -1.0
        Q1 *= -1.0 * rotational_integrator_jacobian_orientation(qb, -ϕb, timestep, attjac=false)
        Q *= LVᵀmat(qb)
        Q1 *= LVᵀmat(qb)
        J = 1.0 / timestep * A * [X Q]
        J += 1.0 / timestep * A * [X1 Q1]
    end

    return J
end

function minimal_velocities_jacobian_velocity(relative::Symbol, joint::Translational,
    xa::AbstractVector, va::AbstractVector, qa::Quaternion, ϕa::AbstractVector,
    xb::AbstractVector, vb::AbstractVector, qb::Quaternion, ϕb::AbstractVector,
    timestep)
    A = nullspace_mask(joint)

    # 1 step backward in time
    xa1 = next_position(xa, -va, timestep)
    qa1 = next_orientation(qa, -ϕa, timestep)
    xb1 = next_position(xb, -vb, timestep)
    qb1 = next_orientation(qb, -ϕb, timestep)

    # Coordinates
    Δx = A * displacement(joint, xa, qa, xb, qb)

    # Previous step coordinates
    Δx1 = A * displacement(joint, xa1, qa1, xb1, qb1)

    # Finite difference
    Δv = (Δx - Δx1) / timestep

    if relative == :parent
        X1, Q1 = displacement_jacobian_configuration(:parent, joint, xa1, qa1, xb1, qb1, attjac=false)
        X1 *= -timestep
        Q1 *= -rotational_integrator_jacobian_velocity(qa, -ϕa, timestep)
        J = -1.0 / timestep * A * [X1 Q1]
    elseif relative == :child
        X1, Q1 = displacement_jacobian_configuration(:child, joint, xa1, qa1, xb1, qb1, attjac=false)
        X1 *= -timestep
        Q1 *= -rotational_integrator_jacobian_velocity(qb, -ϕb, timestep)
        J = -1.0 / timestep * A * [X1 Q1]
    end

    return J
end


text:
using MaskedArrays
using Test

@testset "MaskedArrays.jl" begin
    # Write your tests here.
end


text:
@show ENV["PATH"]
@show ENV["HOME"]
@show ENV["USER"]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
toString(x ::Real)                                                     = string(typeof(x), "(", isa(x, Bool) ? Int(x) : x, ")")
toString(x ::String)                                                   = string(typeof(x), "(\"", x, "\")")
toString(x ::Array)                                                    = isempty(x) ? string(x) : string(eltype(x), x)
toData(s   ::String)                                                   = eval(parse(s))
writeKV(io::IOStream, key::String, value::Union{Real, String, Array}) = write(io, key * "@" * toString(value) * "\n")
readKV(text::String) = begin
  key, value = split(text, '@')
  String(key), eval(parse(value[1:end-1]))
end
checkPath(path::String) = isabspath(path) || !('/' in path || '\\' in path) || error(path * "'s format is error")


text:
@inline match_target_exposure(target_exposure::Price, dir::TradeDir, ba::BidAsk) = target_exposure / open_price(dir, ba)

@inline function update_pnl!(pos::Position, ba::BidAsk)
    # update last market price
    pos.last_quote = ba
    pos.last_dt = ba.dt
    pos.last_price = close_price(pos.dir, ba)

    # update P&L
    pos.pnl = pnl_net(pos)
    return
end


text:
module Particles

include("particle_list.jl")

export ParticleList

include("poisson.jl")

export PoissonSolver
export solve!, eval_density, eval_potential, eval_field

include("poisson_fft.jl")

export PoissonSolverFFT

include("splines.jl")
include("poisson_splines.jl")

export PoissonSolverPBSplines

include("vlasov_poisson.jl")

export VPIntegratorParameters, VPIntegratorCache, integrate_vp!

include("sampling.jl")

export draw_g_accept_reject, draw_g_importance_sampling, weight_f

include("visualisation.jl")

export plot_particles, plot_distribution

include("bump_on_tail.jl")

end # module


text:
function gschurtest(A::Matrix{T}, B::Matrix{T}, tol;
                   normal=false, commuting=false) where {T<:Complex}
    n = size(A,1)
    ulp = eps(real(T))
    if real(T) <: BlasFloat
        S = GenericSchur.ggschur!(copy(A), copy(B))
    else
        S = schur(A, B)
    end
    # test 1: S.T is upper triangular
    @test norm(tril(S.S,-1)) / (n * norm(A) * ulp) < tol
    @test norm(tril(S.T,-1)) / (n * norm(A) * ulp) < tol
    # @test all(tril(S.S,-1) .== 0)
    # @test all(tril(S.T,-1) .== 0)
    # test 2: norm(A - S.Q * S.T * S.Z') / (n * norm(A) * ulp) < tol
    @test norm(A - S.Q * S.S * S.Z') / (n * norm(A) * ulp) < tol
    @test norm(B - S.Q * S.T * S.Z') / (n * norm(A) * ulp) < tol
    # test 3: S.Z is unitary: norm(I - S.Z * S.Z') / (n * ulp) < tol
    @test norm(I - S.Z * S.Z') / (n * ulp) < tol
    @test norm(I - S.Q * S.Q') / (n * ulp) < tol
    # test 4: S.values are e.v. of T
#    @test all(csort(S.values) .== csort(diag(S.T)))

    # It is tempting to check eigenvalues against LAPACK eigvals(A)
    # for suitable types, but the comparison is misleading without
    # detailed analysis because many of the test matrices are severly
    # non-normal.

    if normal && commuting
        # verify that Schur diagonalizes normal commuting matrices
        @test norm(triu(S.S,1)) / (n * norm(B) * ulp) < tol
        @test norm(triu(S.T,1)) / (n * norm(A) * ulp) < tol
    end

    VR = eigvecs(S)
    evcheck = norm((A * VR) * diagm(0 => S.β) - (B * VR) * diagm(0 => S.α))
    @test evcheck / (n*norm(A) * ulp) < tol
    # someday maybe left eigenvectors too
end

@testset "generalized sanity $T" for T in [ComplexF64, Complex{BigFloat}]

unfl = floatmin(real(T))
ovfl = one(real(T)) / unfl
ulp = eps(real(T))
ulpinv = one(real(T)) / ulp
rtulp = sqrt(ulp)
rtulpi = one(real(T)) / rtulp

# sizes of matrices to test
ens = [4,32]

    for n in ens
        A = rand(T,n,n)
        B = rand(T,n,n)
        gschurtest(A, B, 20)
    end
end # testset

#end # module


text:
"""
Abstract Geometry in R{Dim} with Number type T
"""
abstract type AbstractGeometry{Dim,T<:Number} end
abstract type GeometryPrimitive{Dim,T} <: AbstractGeometry{Dim,T} end
Base.ndims(x::AbstractGeometry{Dim}) where {Dim} = Dim

"""
Geometry made of N connected points. Connected as one flat geometry, it makes a Ngon / Polygon.
Connected as volume it will be a Simplex / Tri / Cube.
Note That `Polytope{N} where N == 3` denotes a Triangle both as a Simplex or Ngon.
"""
abstract type Polytope{Dim,T} <: AbstractGeometry{Dim,T} end
abstract type AbstractPolygon{Dim,T} <: Polytope{Dim,T} end

abstract type AbstractPoint{Dim,T} <: StaticVector{Dim,T} end
abstract type AbstractFace{N,T} <: StaticVector{N,T} end
abstract type AbstractSimplexFace{N,T} <: AbstractFace{N,T} end
abstract type AbstractNgonFace{N,T} <: AbstractFace{N,T} end

abstract type AbstractSimplex{Dim,N,T} <: StaticVector{Dim,T} end

"""
Face index, connecting points to form a simplex
"""

@fixed_vector SimplexFace AbstractSimplexFace
const TetrahedronFace{T} = SimplexFace{4,T}
Face(::Type{<:SimplexFace{N}}, ::Type{T}) where {N,T} = SimplexFace{N,T}

"""
Face index, connecting points to form an Ngon
"""

@fixed_vector NgonFace AbstractNgonFace
const LineFace{T} = NgonFace{2,T}
const TriangleFace{T} = NgonFace{3,T}
const QuadFace{T} = NgonFace{4,T}

function Base.show(io::IO, x::TriangleFace{T}) where {T}
    return print(io, "TriangleFace(", join(x, ", "), ")")
end

Face(::Type{<:NgonFace{N}}, ::Type{T}) where {N,T} = NgonFace{N,T}
Face(F::Type{NgonFace{N,FT}}, ::Type{T}) where {FT,N,T} = F

@propagate_inbounds Base.getindex(x::Polytope, i::Integer) = coordinates(x)[i]
@propagate_inbounds Base.iterate(x::Polytope) = iterate(coordinates(x))
@propagate_inbounds Base.iterate(x::Polytope, i) = iterate(coordinates(x), i)

"""
Fixed Size Polygon, e.g.
- N 1-2 : Illegal!
- N = 3 : Triangle
- N = 4 : Quadrilateral (or Quad, Or tetragon)
- N = 5 : Pentagon
- ...
"""
struct Ngon{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: AbstractPolygon{Dim,T}

    points::SVector{N,Point}
end

const NNgon{N} = Ngon{Dim,T,N,P} where {Dim,T,P}

function (::Type{<:NNgon{N}})(points::Vararg{P,N}) where {P<:AbstractPoint{Dim,T},
                                                          N} where {Dim,T}
    return Ngon{Dim,T,N,P}(SVector(points))
end
Base.show(io::IO, x::NNgon{N}) where {N} = print(io, "Ngon{$N}(", join(x, ", "), ")")

# Interfaces
coordinates(x::Ngon) = x.points
# Base Array interface
Base.length(::Type{<:NNgon{N}}) where {N} = N
Base.length(::NNgon{N}) where {N} = N

"""
The Ngon Polytope element type when indexing an array of points with a SimplexFace
"""
function Polytope(P::Type{<:AbstractPoint{Dim,T}},
                  ::Type{<:AbstractNgonFace{N,IT}}) where {N,Dim,T,IT}
    return Ngon{Dim,T,N,P}
end

"""
The fully concrete Ngon type, when constructed from a point type!
"""
function Polytope(::Type{<:NNgon{N}}, P::Type{<:AbstractPoint{NDim,T}}) where {N,NDim,T}
    return Ngon{NDim,T,N,P}
end

const LineP{Dim,T,P<:AbstractPoint{Dim,T}} = Ngon{Dim,T,2,P}
const Line{Dim,T} = LineP{Dim,T,Point{Dim,T}}

# Simplex{D, T, 3} & Ngon{D, T, 3} are both representing a triangle.
# Since Ngon is supposed to be flat and a triangle is flat, lets prefer Ngon
# for triangle:
const TriangleP{Dim,T,P<:AbstractPoint{Dim,T}} = Ngon{Dim,T,3,P}
const Triangle{Dim,T} = TriangleP{Dim,T,Point{Dim,T}}
const Triangle3d{T} = Triangle{3,T}

Base.show(io::IO, x::TriangleP) = print(io, "Triangle(", join(x, ", "), ")")
Base.summary(io::IO, x::Type{<:TriangleP}) = print(io, "Triangle")

const Quadrilateral{Dim,T} = Ngon{Dim,T,4,P} where {P<:AbstractPoint{Dim,T}}

Base.show(io::IO, x::Quadrilateral) = print(io, "Quad(", join(x, ", "), ")")
Base.summary(io::IO, x::Type{<:Quadrilateral}) = print(io, "Quad")

function coordinates(lines::AbstractArray{LineP{Dim,T,PointType}}) where {Dim,T,PointType}
    return if lines isa Base.ReinterpretArray
        return coordinates(lines.parent)
    else
        result = PointType[]
        for line in lines
            append!(result, coordinates(line))
        end
        return result
    end
end

"""
A `Simplex` is a generalization of an N-dimensional tetrahedra and can be thought
of as a minimal convex set containing the specified points.

* A 0-simplex is a point.
* A 1-simplex is a line segment.
* A 2-simplex is a triangle.
* A 3-simplex is a tetrahedron.

Note that this datatype is offset by one compared to the traditional
mathematical terminology. So a one-simplex is represented as `Simplex{2,T}`.
This is for a simpler implementation.

It applies to infinite dimensions. The structure of this type is designed
to allow embedding in higher-order spaces by parameterizing on `T`.
"""
struct Simplex{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: Polytope{Dim,T}

    points::SVector{N,Point}
end

const NSimplex{N} = Simplex{Dim,T,N,P} where {Dim,T,P}
const TetrahedronP{T,P<:AbstractPoint{3,T}} = Simplex{3,T,4,P}
const Tetrahedron{T} = TetrahedronP{T,Point{3,T}}

Base.show(io::IO, x::TetrahedronP) = print(io, "Tetrahedron(", join(x, ", "), ")")

coordinates(x::Simplex) = x.points

function (::Type{<:NSimplex{N}})(points::Vararg{P,N}) where {P<:AbstractPoint{Dim,T},
                                                             N} where {Dim,T}
    return Simplex{Dim,T,N,P}(SVector(points))
end

# Base Array interface
Base.length(::Type{<:NSimplex{N}}) where {N} = N
Base.length(::NSimplex{N}) where {N} = N

"""
The Simplex Polytope element type when indexing an array of points with a SimplexFace
"""
function Polytope(P::Type{<:AbstractPoint{Dim,T}},
                  ::Type{<:AbstractSimplexFace{N}}) where {N,Dim,T}
    return Simplex{Dim,T,N,P}
end

"""
The fully concrete Simplex type, when constructed from a point type!
"""
function Polytope(::Type{<:NSimplex{N}}, P::Type{<:AbstractPoint{NDim,T}}) where {N,NDim,T}
    return Simplex{NDim,T,N,P}
end
Base.show(io::IO, x::LineP) = print(io, "Line(", x[1], " => ", x[2], ")")

"""
    LineString(points::AbstractVector{<:AbstractPoint})

A LineString is a geometry of connected line segments
"""
struct LineString{Dim,T<:Real,P<:AbstractPoint,V<:AbstractVector{<:LineP{Dim,T,P}}} <:
       AbstractVector{LineP{Dim,T,P}}
    points::V
end

coordinates(x::LineString) = coordinates(x.points)

Base.copy(x::LineString) = LineString(copy(x.points))
Base.size(x::LineString) = size(getfield(x, :points))
Base.getindex(x::LineString, i) = getindex(getfield(x, :points), i)

function LineString(points::AbstractVector{LineP{Dim,T,P}}) where {Dim,T,P}
    return LineString{Dim,T,P,typeof(points)}(points)
end

"""
    LineString(points::AbstractVector{<: AbstractPoint}, skip = 1)

Creates a LineString from a vector of points
With `skip == 1`, the default, it will connect the line like this:
```julia
points = Point[a, b, c, d]
linestring = LineString(points)
@assert linestring == LineString([a => b, b => c, c => d])
```
"""
function LineString(points::AbstractVector{<:AbstractPoint}, skip=1)
    return LineString(connect(points, LineP, skip))
end

function LineString(points::AbstractVector{<:Pair{P,P}}) where {P<:AbstractPoint{N,T}} where {N,
                                                                                T}
    return LineString(reinterpret(LineP{N,T,P}, points))
end

function LineString(points::AbstractVector{<:AbstractPoint},
                    faces::AbstractVector{<:LineFace})
    return LineString(connect(points, faces))
end

"""
    LineString(points::AbstractVector{<: AbstractPoint}, indices::AbstractVector{<: Integer}, skip = 1)

Creates a LineString from a vector of points and an index list.
With `skip == 1`, the default, it will connect the line like this:

```julia
points = Point[a, b, c, d]; faces = [1, 2, 3, 4]
linestring = LineString(points, faces)
@assert linestring == LineString([a => b, b => c, c => d])
```
To make a segmented line, set skip to 2
```julia
points = Point[a, b, c, d]; faces = [1, 2, 3, 4]
linestring = LineString(points, faces, 2)
@assert linestring == LineString([a => b, c => d])
```
"""
function LineString(points::AbstractVector{<:AbstractPoint},
                    indices::AbstractVector{<:Integer}, skip=1)
    faces = connect(indices, LineFace, skip)
    return LineString(points, faces)
end

"""
    Polygon(exterior::AbstractVector{<:Point})
    Polygon(exterior::AbstractVector{<:Point}, interiors::Vector{<:AbstractVector{<:AbstractPoint}})

"""
struct Polygon{Dim,T<:Real,P<:AbstractPoint{Dim,T},L<:AbstractVector{<:LineP{Dim,T,P}},
               V<:AbstractVector{L}} <: AbstractPolygon{Dim,T}
    exterior::L
    interiors::V
end

Base.copy(x::Polygon) = Polygon(copy(x.exterior), copy(x.interiors))

function Base.:(==)(a::Polygon, b::Polygon)
    return (a.exterior == b.exterior) && (a.interiors == b.interiors)
end

function Polygon(exterior::E,
                 interiors::AbstractVector{E}) where {E<:AbstractVector{LineP{Dim,T,P}}} where {Dim,
                                                                                T,
                                                                                P}
    return Polygon{Dim,T,P,typeof(exterior),typeof(interiors)}(exterior, interiors)
end

Polygon(exterior::L) where {L<:AbstractVector{<:LineP}} = Polygon(exterior, L[])

function Polygon(exterior::AbstractVector{P},
                 skip::Int=1) where {P<:AbstractPoint{Dim,T}} where {Dim,T}
    return Polygon(LineString(exterior, skip))
end

function Polygon(exterior::AbstractVector{P}, faces::AbstractVector{<:Integer},
                 skip::Int=1) where {P<:AbstractPoint{Dim,T}} where {Dim,T}
    return Polygon(LineString(exterior, faces, skip))
end

function Polygon(exterior::AbstractVector{P},
                 faces::AbstractVector{<:LineFace}) where {P<:AbstractPoint{Dim,T}} where {Dim,
                                                                                T}
    return Polygon(LineString(exterior, faces))
end

function Polygon(exterior::AbstractVector{P},
                 interior::AbstractVector{<:AbstractVector{P}}) where {P<:AbstractPoint{Dim,
                                                                                T}} where {Dim,
                                                                                T}
    ext = LineString(exterior)
    # We need to take extra care for empty interiors, since
    # if we just map over it, it won't infer the element type correctly!
    int = typeof(ext)[]
    foreach(x -> push!(int, LineString(x)), interior)
    return Polygon(ext, int)
end

function coordinates(polygon::Polygon{N,T,PointType}) where {N,T,PointType}
    exterior = coordinates(polygon.exterior)
    if isempty(polygon.interiors)
        return exterior
    else
        result = PointType[]
        append!(result, exterior)
        foreach(x -> append!(result, coordinates(x)), polygon.interiors)
        return result
    end
end

"""
    MultiPolygon(polygons::AbstractPolygon)
"""
struct MultiPolygon{Dim,T<:Real,Element<:AbstractPolygon{Dim,T},
                    A<:AbstractVector{Element}} <: AbstractVector{Element}
    polygons::A
end

function MultiPolygon(polygons::AbstractVector{P};
                      kw...) where {P<:AbstractPolygon{Dim,T}} where {Dim,T}
    return MultiPolygon(meta(polygons; kw...))
end

Base.getindex(mp::MultiPolygon, i) = mp.polygons[i]
Base.size(mp::MultiPolygon) = size(mp.polygons)

struct MultiLineString{Dim,T<:Real,Element<:LineString{Dim,T},A<:AbstractVector{Element}} <:
       AbstractVector{Element}

    linestrings::A
end

function MultiLineString(linestrings::AbstractVector{L};
                         kw...) where {L<:AbstractVector{LineP{Dim,T,P}}} where {Dim,T,P}
    return MultiLineString(meta(linestrings; kw...))
end

Base.getindex(ms::MultiLineString, i) = ms.linestrings[i]
Base.size(ms::MultiLineString) = size(ms.linestrings)

"""
    MultiPoint(points::AbstractVector{AbstractPoint})

A collection of points
"""
struct MultiPoint{Dim,T<:Real,P<:AbstractPoint{Dim,T},A<:AbstractVector{P}} <:
       AbstractVector{P}

    points::A
end

function MultiPoint(points::AbstractVector{P};
                    kw...) where {P<:AbstractPoint{Dim,T}} where {Dim,T}
    return MultiPoint(meta(points; kw...))
end

Base.getindex(mpt::MultiPoint, i) = mpt.points[i]
Base.size(mpt::MultiPoint) = size(mpt.points)

"""
    AbstractMesh

An abstract mesh is a collection of Polytope elements (Simplices / Ngons).
The connections are defined via faces(mesh), the coordinates of the elements are returned by
coordinates(mesh). Arbitrary meta information can be attached per point or per face
"""
abstract type AbstractMesh{Element<:Polytope} <: AbstractVector{Element} end

"""
    Mesh <: AbstractVector{Element}
The conrecte AbstractMesh implementation
"""
struct Mesh{Dim,T<:Number,Element<:Polytope{Dim,T},V<:AbstractVector{Element}} <:
       AbstractMesh{Element}
    simplices::V # usually a FaceView, to connect a set of points via a set of faces.
end

Tables.schema(mesh::Mesh) = Tables.schema(getfield(mesh, :simplices))

function Base.getproperty(mesh::Mesh, name::Symbol)
    if name === :position
        # a mesh always has position defined by coordinates...
        return metafree(coordinates(mesh))
    else
        return getproperty(getfield(mesh, :simplices), name)
    end
end

function Base.propertynames(mesh::Mesh)
    names = propertynames(getfield(mesh, :simplices))
    if :position in names
        return names
    else
        # a mesh always has positions!
        return (names..., :position)
    end
end

function Base.summary(io::IO, ::Mesh{Dim,T,Element}) where {Dim,T,Element}
    print(io, "Mesh{$Dim, $T, ")
    summary(io, Element)
    return print(io, "}")
end

Base.size(mesh::Mesh) = size(getfield(mesh, :simplices))
Base.getindex(mesh::Mesh, i::Integer) = getfield(mesh, :simplices)[i]

function Mesh(elements::AbstractVector{<:Polytope{Dim,T}}) where {Dim,T}
    return Mesh{Dim,T,eltype(elements),typeof(elements)}(elements)
end

function Mesh(points::AbstractVector{<:AbstractPoint},
              faces::AbstractVector{<:AbstractFace})
    return Mesh(connect(points, faces))
end

function Mesh(points::AbstractVector{<:AbstractPoint}, faces::AbstractVector{<:Integer},
              facetype=TriangleFace, skip=1)
    return Mesh(connect(points, connect(faces, facetype, skip)))
end


text:
"""

$(TYPEDEF)

Structure to contain data needed to compute the mddf for a single frame

$(TYPEDFIELDS)

"""
mutable struct FrameData{T<:Trajectory,V}

    trajectory::T
    volume_frame::Volume
    rdf_count_random_frame::Vector{Float64}
    md_count_random_frame::Vector{Float64}

    dc::CutoffDistances
    dmin_mol::Vector{DminMol}
    dref_mol::Vector{Float64}
    x_solvent_random::Vector{V}

    lc_solvent::LinkedCells

end

function FrameData(trajectory::Trajectory, R::Result)
    return FrameData(
        trajectory,                                       # trajectory
        Volume(R.nbins),                                  # volume_frame
        zeros(R.nbins),                                   # rdf_count_random_frame
        zeros(R.nbins),                                   # md_count_random_frame
        CutoffDistances(trajectory.solvent.natoms),       # dc
        [DminMol(+Inf, i, 0, 0) for i = 1:trajectory.solvent.nmols], # dmin_mol
        zeros(trajectory.solvent.nmols),                             # dref_mol
        similar(trajectory.x_solvent),                               # x_solvent_random
        LinkedCells(trajectory.solvent.natoms),
    )                      # lc_solvent
end