text:
# OpenCV Original License

#
#
#  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
#
#  By downloading, copying, installing or using the software you agree to this license.
#  If you do not agree to this license, do not download, install,
#  copy or use the software.
#
#
#                           License Agreement
#                For Open Source Computer Vision Library
#
# Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
# Copyright (C) 2009, Willow Garage Inc., all rights reserved.
# Third party copyrights are property of their respective owners.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
#   * Redistribution's of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#
#   * Redistribution's in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#
#   * The name of the copyright holders may not be used to endorse or promote products
#     derived from this software without specific prior written permission.
#
# This software is provided by the copyright holders and contributors "as is" and
# any express or implied warranties, including, but not limited to, the implied
# warranties of merchantability and fitness for a particular purpose are disclaimed.
# In no event shall the Intel Corporation or contributors be liable for any direct,
# indirect, incidental, special, exemplary, or consequential damages
# (including, but not limited to, procurement of substitute goods or services;
# loss of use, data, or profits; or business interruption) however caused
# and on any theory of liability, whether in contract, strict liability,
# or tort (including negligence or otherwise) arising in any way out of
# the use of this software, even if advised of the possibility of such damage.
#
#M

using Images

"""

creategridcoordinates(xnum::Integer, ynum::Integer)::Matrix{<:Integer}

Creates a vector of grid 0-indexed coordinates up to xnum in the first dimension and up to ynum in the second dimension.

# Examples
```julia-repl
julia> creategridcoordinates(3,5)
```
"""
function creategridcoordinates(xnum::Integer, ynum::Integer)
    n = xnum * ynum
    res = zeros(Int32, n, 2)
    for i = 0:(n - 1)
        res[i + 1, :] = [i ÷ ynum, i % ynum]
    end
    return res
end

"""

creategridcoordinates(xlength::Real, ylength::Real, xnum::Integer, ynum::Integer)

Creates a vector of grid coordinates which span the rowsize in the first dimension and the colsize in the second dimension.

# Arguments
- `xlength::Real`: Length to span in first dimension
- `ylength::Real`: Length to span in second dimension
- `xnum::Integer`: Number of points to use when spanning first dimension
- `ynum::Integer`: Number of points to use when spanning second dimension

# Examples
```julia-repl
julia> creategridcoordinates(640, 480, 10, 10)
```
"""
function creategridcoordinates(xlength::Real, ylength::Real, xnum::Integer, ynum::Integer)
    grid = creategridcoordinates(xnum, ynum)
    grid = grid .* [(xlength / (xnum - 1)) (ylength / (ynum - 1))]
    return grid
end

function distortpoints(src::AbstractMatrix{<:Real}; cameramatrix::AbstractMatrix{<:Real}, distcoeffs::AbstractVector{<:Real}, newcameramatrix::AbstractMatrix{<:Real})
    fx1 = newcameramatrix[1,1]
    fy1 = newcameramatrix[2,2]
    cx1 = newcameramatrix[1,3]
    cy1 = newcameramatrix[2,3]
    ifx1 = 1 / fx1
    ify1 = 1 / fy1

    fx0 = cameramatrix[1,1]
    fy0 = cameramatrix[2,2]
    cx0 = cameramatrix[1,3]
    cy0 = cameramatrix[2,3]

    k1 = distcoeffs[1]
    k2 = distcoeffs[2]
    k3 = distcoeffs[5]
    p1 = distcoeffs[3]
    p2 = distcoeffs[4]

    # undo linear distortion
    dst = (src .- [cx1 cy1]) .* [ifx1 ify1]

    # perform nonlinear distortion
    for i = 1:size(src)[1]
        x = dst[i, 1]
        y = dst[i, 2]

        x2 = x * x
        y2 = y * y

        r2 = x2 + y2
        r4 = r2 * r2
        r6 = r4 * r2

        rad = 1 + k1 * r2 + k2 * r4 + k3 * r6

        xp = x * rad + 2p1 * x * y + p2 * (r2 + 2x2)
        yp = y * rad + 2p2 * x * y + p1 * (r2 + 2y2)

        dst[i, :] = [xp yp]
    end

    # distort linearly
    dst = dst .* [fx0 fy0] .+ [cx0 cy0]
    return dst
end

"""

    undistortpoints(src::AbstractMatrix{<:Real}; cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real}, Nothing}, newcameramatrix::Union{AbstractMatrix{<:Real}, Nothing})

Computes the ideal point coordinates from the observed point coordinates.

# Arguments
- `src::AbstractMatrix{<:Real}`: Observerd point coordinates, (N,2)
- `cameramatrix::AbstractMatrix{<:Real}`: Camera matrix, (3,3)
- `distcoeffs::Union{AbstractVector{<:Real}, Nothing}`: Input vector of distortion coefficients, (5,) (k1, k2, p1, p2, k3)
- `newcameramatrix::Union{AbstractMatrix{<:Real}, Nothing}`: New camera matrix, (3,3). If the matrix is empty, the identity new camera matrix is used.

# Examples
```julia-repl
julia> undistortpoints([10 23; 100 100.3; 23.1 20.5], cameramatrix=[517.3 0 318.6; 0 516.5 244.3; 0 0 1], distcoeffs=[0.2624, -0.9531, -0.0054, 0.0026, 1.1633])
```
"""
function undistortpoints(src::AbstractMatrix{<:Real}; cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real},Nothing}, newcameramatrix::Union{AbstractMatrix{<:Real},Nothing})::AbstractMatrix{<:Real}
    # sanitize inputs
    if size(src)[2] != 2
        @error "src should have 2 columns"
        return
    end

    if size(cameramatrix) != (3, 3)
        @error "cameramatrix should be of size (3,3)"
        return
    end

    if newcameramatrix != nothing && size(newcameramatrix) != (3, 3)
        @error "newcameramatrix should be of size (3,3)"
        return
    end

    if distcoeffs != nothing && size(distcoeffs) != (5,)
        @error "distcoeffs must be of size (5,)"
        return
    end

    fx = cameramatrix[1,1]
    fy = cameramatrix[2,2]
    cx = cameramatrix[1,3]
    cy = cameramatrix[2,3]

    ifx = 1 / fx
    ify = 1 / fy

    # undo linear distortion
    dst = (src .- [cx cy]) .* [ifx ify]

    # undo nonlinear distortion numerically
    num_points = size(src)[1]
    if distcoeffs != nothing

        k1 = distcoeffs[1]
        k2 = distcoeffs[2]
        k3 = distcoeffs[5]
        p1 = distcoeffs[3]
        p2 = distcoeffs[4]

        for i = 1:num_points
            x0 = dst[i, 1]
            y0 = dst[i, 2]
            x = x0
            y = y0
            for j = 1:5
                r2 = x * x + y * y
                r4 = r2 * r2
                r6 = r4 * r2
                icdist = 1 / (1 + k1 * r2 + k2 * r4 + k3 * r6)
                deltaX = 2 * p1 * x * y + p2 * (r2 + 2 * x * x)
                deltaY = p1 * (r2 + 2 * y * y) + 2 * p2 * x * y
                x = (x0 - deltaX) * icdist
                y = (y0 - deltaY) * icdist
            end
            dst[i, :] = [x y]
        end
    end

    # linear distortion according to newcameramatrix
    if newcameramatrix != nothing
        fx = newcameramatrix[1,1]
        fy = newcameramatrix[2,2]
        cx = newcameramatrix[1,3]
        cy = newcameramatrix[2,3]

        dst = dst .* [fx fy] .+ [cx cy]
    end
    return dst

end

struct Rectangle
    x::Real
    y::Real
    w::Real
    h::Real
end

"""

    getrectangles(cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real}, Nothing}, newcameramatrixUnion::{AbstractVector{<:Real}, Nothing}, imgsize::Tuple{Integer, Integer})

Compute the inscribed and circumscribed rectangle of a set of points in an image after undistorting them.

Returns outerrectangle, innerrectangle

# Arguments
- `cameramatrix::AbstractMatrix{<:Real}`: Camera matrix, (3,3)
- `distcoeffs::Union{AbstractVector{<:Real}, Nothing}`: Input vector of distortion coefficients, (5,) (k1, k2, p1, p2, k3)
- `newcameramatrix::Union{AbstractMatrix{<:Real}, Nothing}`: New camera matrix, (3,3). If the matrix is empty, the identity new camera matrix is used.
- `imgsize::Tuple{Integer, Integer}`: size of image represented as a tuple

# Examples
```julia-repl
julia> getrectangles(cameramatrix=[517.3 0 318.6; 0 516.5 244.3; 0 0 1], distcoeffs=[0.2624, -0.9531, -0.0054, 0.0026, 1.1633], imgsize=(640, 480))
```
"""
function getrectangles(;cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real},Nothing} = nothing, newcameramatrix::Union{AbstractMatrix{<:Real},Nothing} = nothing, imgsize::Tuple{Integer,Integer})::Tuple{Rectangle,Rectangle}
    N = 9
    pts = creategridcoordinates(imgsize[1], imgsize[2], N, N)
    pts = undistortpoints(pts, cameramatrix = cameramatrix, distcoeffs = distcoeffs, newcameramatrix = newcameramatrix)

    iX0 = -Inf; iX1 = Inf; iY0 = -Inf; iY1 = Inf
    oX0 = Inf; oX1 = -Inf; oY0 = Inf; oY1 = -Inf

    i = 1
    for x = 1:N
        for y = 1:N
            px = pts[i,1]
            py = pts[i,2]

            oX0 = min(oX0, px)
            oX1 = max(oX1, px)
            oY0 = min(oY0, py)
            oY1 = max(oY1, py)

            if x == 1
                iX0 = max(iX0, px)
            end
            if x == N
                iX1 = min(iX1, px)
            end
            if y == 1
                iY0 = max(iY0, py)
            end
            if y == N
                iY1 = min(iY1, py)
            end
            i += 1
        end
    end
    return Rectangle(oX0, oY0, oX1 - oX0, oY1 - oY0), Rectangle(iX0, iY0, iX1 - iX0, iY1 - iY0)

end

"""

    getoptimal_newcameramatrix(;cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real}, Nothing} = nothing, imgsize::Tuple{Integer, Integer}, alpha::Real, newimgsize::Tuple{Integer, Integer})

Compute the optimal camera matrix such that after removing nonlinear distortion the image is remapped nicely.

# Arguments
- `cameramatrix::AbstractMatrix{<:Real}`: Camera matrix, (3,3)
- `distcoeffs::Union{AbstractVector{<:Real}, Nothing}`: Input vector of distortion coefficients, (5,) (k1, k2, p1, p2, k3)
- `newcameramatrix::Union{AbstractMatrix{<:Real}, Nothing}`: New camera matrix, (3,3). If the matrix is empty, the identity new camera matrix is used.
- `imgsize::Tuple{Integer, Integer}`: size of image represented as a tuple
- `alpha::Real`: interpolation factor between "remap without missing regions" to "remap fitting all data"
- `newimgsize::Tuple{Integer, Integer}`: desired size of image after remapping

"""
function get_optimalnewcameramatrix(;cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real},Nothing} = nothing, imgsize::Tuple{Integer,Integer}, alpha::Real, newimgsize::Tuple{Integer,Integer})

    # Get inscribed and circumscribed rectangles in normalized
    outer, inner = getrectangles(cameramatrix = cameramatrix, distcoeffs = distcoeffs, imgsize = imgsize)

    # Projection mapping inner rectangle to viewport
    fx0 = newimgsize[1] / inner.w
    fy0 = newimgsize[2] / inner.h
    cx0 = -fx0 * inner.x
    cy0 = -fy0 * inner.y

    # Projection mapping outer rectangle to viewport
    fx1 = newimgsize[1] / outer.w
    fy1 = newimgsize[2] / outer.h
    cx1 = -fx1 * outer.x
    cy1 = -fy1 * outer.y

    # Interpolate between the two optimal projections
    M = zeros(3, 3) # optimal camera matrix
    M[1,1] = fx0 * (1 - alpha) + fx1 * alpha
    M[2,2] = fy0 * (1 - alpha) + fy1 * alpha
    M[1,3] = cx0 * (1 - alpha) + cx1 * alpha
    M[2,3] = cy0 * (1 - alpha) + cy1 * alpha
    M[3,3] = 1
    return M
end

"""

    getoptimal_newcameramatrix(;cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real}, Nothing} = nothing, imgsize::Tuple{Integer, Integer}, alpha::Real, newimgsize::Tuple{Integer, Integer})

Computes the undistortion and rectification transformation map. Returns a tuple of maps in the form of matrices.

# Returns
- mapx: (u,v) -> x
- mapy: (u,v) -> y

# Arguments
- `cameramatrix::AbstractMatrix{<:Real}`: Camera matrix, (3,3)
- `distcoeffs::Union{AbstractVector{<:Real}, Nothing}`: Input vector of distortion coefficients, (5,) (k1, k2, p1, p2, k3)
- `newcameramatrix::Union{AbstractMatrix{<:Real}, Nothing}`: New camera matrix, (3,3). If the matrix is empty, the identity new camera matrix is used.
- `imgsize::Tuple{Integer, Integer}`: size of image represented as a tuple

"""
function init_undistortrectifymap(;cameramatrix::AbstractMatrix{<:Real}, distcoeffs::Union{AbstractVector{<:Real},Nothing}, newcameramatrix::AbstractMatrix{<:Real}, imgsize::Tuple{Integer,Integer})
    fx0 = cameramatrix[1,1]
    fy0 = cameramatrix[2,2]
    cx0 = cameramatrix[1,3]
    cy0 = cameramatrix[2,3]
    ifx0 = 1 / fx0
    ify0 = 1 / fy0

    fx1 = newcameramatrix[1,1]
    fy1 = newcameramatrix[2,2]
    cx1 = newcameramatrix[1,3]
    cy1 = newcameramatrix[2,3]
    ifx1 = 1 / fx1
    ify1 = 1 / fy1

    k1 = distcoeffs[1]
    k2 = distcoeffs[2]
    k3 = distcoeffs[5]
    p1 = distcoeffs[3]
    p2 = distcoeffs[4]

    xmap = zeros(Real, imgsize[2], imgsize[1])
    ymap = zeros(Real, imgsize[2], imgsize[1])

    for i = 0:(imgsize[1] - 1)
        x = (i - cx1) * ifx1
        x2 = x * x
        for j = 0:(imgsize[2] - 1)
            y = (j - cy1) * ify1
            y2 = y * y

            r2 = x2 + y2
            r4 = r2 * r2
            r6 = r2 * r4
            rad = (1 + k1 * r2 + k2 * r4 + k3 * r6)

            xp = x * rad + 2p1 * x * y + p2 * (r2 + 2x2)
            yp = y * rad + 2p2 * x * y + p1 * (r2 + 2y2)

            xp = xp * fx0 + cx0
            yp = yp * fy0 + cy0

            xmap[j + 1, i + 1] = xp
            ymap[j + 1, i + 1]  = yp
        end
    end
    return xmap, ymap
end

"""

remap(src::AbstractArray{T,2}, xmap::AbstractArray{<:Real,2}, ymap::AbstractArray{<:Real,2}) where T

Applies a generic geometrical transformation to an image.

# Arguments
- `src::AbstractMatrix{T,2}`: source image
- `xmap::AbstractArray{<:Real,2}: map taking (u,v) in destination to x in src
- `yxmap::AbstractArray{<:Real,2}: map taking (u,v) in destination to y in src

"""
function remap(src::AbstractArray{T,2}, xmap::AbstractArray{<:Real,2}, ymap::AbstractArray{<:Real,2}) where T
    dst = zeros(T, size(xmap))
    for j = 1:size(xmap)[1]
        for i = 1:size(xmap)[2]
            dst[j,i] = bilinear_interpolation(src, ymap[j, i], xmap[j, i])
        end
    end
    dst
end


text:
# Defines various functions for "moving" around
# Pascal's triangle.

_upval(a::Entry{<: Integer}) = a.val*(a.n - a.k)÷a.n
_upval(a::Entry{<: AbstractFloat}) = a.val*(a.n - a.k)/a.n
function up!(a::Entry)
    if isatright(a) || isfirst(a)
        throw(OutOfBoundsError("no entry above"))
    end
    a.val = _upval(a)
    a.n -= 1
    return a
end

up(a::Entry) = up!(Entry(a))

up!(r::Row) = prev!(r)
up(r::Row) = prev(r)

_downval(a::Entry{<: Integer}) = a.val*(a.n + 1)÷(a.n - a.k + 1)
_downval(a::Entry{<: AbstractFloat}) = a.val*(a.n + 1)/(a.n - a.k + 1)
function down!(a::Entry)
    a.val = _downval(a)
    a.n += 1
    return a
end

down(a::Entry) = down!(Entry(a))

down!(r::Row) = next!(r)
down(r::Row) = next(r)

_leftval(a::Entry{<: Integer}) = a.val*a.k÷(a.n - a.k + 1)
_leftval(a::Entry{<: AbstractFloat}) = a.val*a.k/(a.n - a.k + 1)
function _left!(a)
    a.val = _leftval(a)
    a.k -= 1
    return a
end

function left!(a::Entry)
    isatleft(a) && throw(OutOfBoundsError("no entry to the left"))
    return _left!(a)
end
left(a::Entry) = left!(Entry(a))

left!(c::Column) = prev!(c)
left(c::Column) = prev(c)
left!(c::LazyColumn) = prev!(c)
left(c::LazyColumn) = prev(c)

_rightval(a::Entry{<: Integer}) = a.val*(a.n - a.k)÷(a.k + 1)
_rightval(a::Entry{<: AbstractFloat}) = a.val*(a.n - a.k)/(a.k + 1)
function _right!(a::Entry)
    a.val = _rightval(a)
    a.k += 1
    return a
end

function right!(a::Entry)
    isatright(a) && throw(OutOfBoundsError("no entry to the right"))
    return _right!(a)
end
right(a::Entry) = right!(Entry(a))

right!(c::Column) = next!(c)
right(c::Column) = next(c)
right!(c::LazyColumn) = next!(c)
right(c::LazyColumn) = next(c)

function prev!(a::Entry)
    isfirst(a) && throw(OutOfBoundsError("no previous entry"))
    if !isatleft(a)
        return _left!(a)
    else
        a.n -= 1
        a.k = a.n
        a.val = one(a.val)
        return a
    end
end
prev(a::Entry) = prev!(Entry(a))

function prev!(r::Row)
    isfirst(r) && throw(OutOfBoundsError("no previous row"))
    r.rownum -= 1
    datalength = numelements(r.rownum)
    if datalength ≥ 1
        r.data[1] -= r.rownum
        for i ∈ 2:datalength
            r.data[i] -= r.data[i-1]
        end
    end
    return r
end

function prev(r::Row)
    newrow = Row(r)
    return prev!(newrow)
end

function prev!(c::Column)
    isfirst(c) && throw(OutOfBoundsError("no previous column"))
    c.colnum -= 1
    for i ∈ length(c.data):-1:2
        c.data[i] -= c.data[i-1]
    end
    return c
end

function prev(c::Column)
    newcol = Column(c)
    return prev!(newcol)
end

function prev!(c::LazyColumn)
    isfirst(c) && throw(OutOfBoundsError("no previous column"))
    c.colnum -= 1
    for (r, val) ∈ c.data
        c.data[r] = value(up!(Entry(r + c.colnum, r - 1, val)))
    end
    return c
end

function prev(c::LazyColumn)
    newcol = LazyColumn(c)
    return prev!(newcol)
end

function next!(a::Entry)
    if !isatright(a)
        return _right!(a)
    else
        a.n += 1
        a.k = zero(a.k)
        a.val = one(a.val)
        return a
    end
end
next(a::Entry) = next!(Entry(a))

function next!(r::Row)
    datalength = numelements(r.rownum)
    if isodd(r.rownum) && r.rownum ≥ 3
        if r.rownum == 3
            newdata = one(eltype(r.data))*6
        else
            newdata = r.data[datalength]
        end
        if length(r.data) > datalength
            r.data[datalength+1] = newdata
        else
            push!(r.data, newdata)
        end
        datalength += 1
    end
    if datalength ≥ 1
        for i ∈ datalength:-1:2
            r.data[i] += r.data[i-1]
        end
        r.data[1] += r.rownum
    end
    r.rownum += 1
    return r
end

function next(r::Row)
    newrow = Row(r)
    return next!(newrow)
end

function next!(c::Column)
    c.colnum += 1
    for i ∈ 2:length(c.data)
        c.data[i] += c.data[i-1]
    end
    return c
end

function next(c::Column)
    newcol = Column(c)
    return next!(newcol)
end

function _downright!(e::Entry)
    if isatright(e)
        e.n += 1
        e.k += 1
        return e
    end
    a = right(e)
    e += a
    return e
end

function next!(c::LazyColumn)
    for (r, val) ∈ c.data
        e = Entry(r + c.colnum - 1, c.colnum, val)
        _downright!(e)
        c.data[r] = value(e)
    end
    c.colnum += 1
    return c
end

function next(c::LazyColumn)
    newcol = LazyColumn(c)
    return next!(newcol)
end

ext: jl


text:
###
### This file contains abstract types, interfaces and traits.
###
"""
    AbstractHamiltonian{T}

Supertype that provides an interface for linear operators over a linear space with scalar
type `T` that are suitable for FCIQMC (with [`lomc!`](@ref)). Indexing is done with addresses (typically not
integers) from an address space that may be large (and will not need to be completely
generated).

`AbstractHamiltonian` instances operate on vectors of type [`AbstractDVec`](@ref) from the
module `DictVectors` and work well with addresses of type [`AbstractFockAddress`](@ref) from
the module `BitStringAddresses`. The type works well with the external package
[KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl).

For available implementations see [`Hamiltonians`](@ref).

# Interface

Methods that need to be implemented:

* [`num_offdiagonals(::AbstractHamiltonian, address)`](@ref)
* [`get_offdiagonal(::AbstractHamiltonian, address, chosen::Integer)`](@ref)
* [`diagonal_element(::AbstractHamiltonian, address)`](@ref)
* [`starting_address(::AbstractHamiltonian)`](@ref)

Optional methods to implement:

* [`LOStructure(::Type{typeof(lo)})`](@ref): defaults to `AdjointUnknown`
* [`dimension(::Type{T}, ::AbstractHamiltonian)`](@ref): defaults to dimension of address
  space
* [`momentum(::AbstractHamiltonian)`](@ref): no default

Provides:

* [`offdiagonals`](@ref): iterator over reachable off-diagonal matrix elements
* [`random_offdiagonal`](@ref): function to generate random off-diagonal matrix element
* `*(H, v)`: deterministic matrix-vector multiply (allocating)
* `H(v)`: equivalent to `H * v`.
* `mul!(w, H, v)`: mutating matrix-vector multiply.
* [`dot(x, H, v)`](@ref): compute `x⋅(H*v)` minimizing allocations.
* `H[address1, address2]`: indexing with `getindex()` - mostly for testing purposes (slow!)
* [`BasisSetRep`](@ref): construct a basis set repesentation
* [`sparse`](@ref), [`Matrix`](@ref): construct a (sparse) matrix representation

See also [`Hamiltonians`](@ref), [`Interfaces`](@ref).
"""
abstract type AbstractHamiltonian{T} end

Base.eltype(::AbstractHamiltonian{T}) where {T} = T

"""
    diagonal_element(ham, add)

Compute the diagonal matrix element of the linear operator `ham` at
address `add`.

# Example

```jldoctest
julia> addr = BoseFS((3, 2, 1));

julia> H = HubbardMom1D(addr);

julia> diagonal_element(H, addr)
8.666666666666664
```
"""
diagonal_element(m::AbstractMatrix, i) = m[i, i]

"""
    num_offdiagonals(ham, add)

Compute the number of number of reachable configurations from address `add`.

# Example

```jldoctest
julia> addr = BoseFS((3, 2, 1));

julia> H = HubbardMom1D(addr);

julia> num_offdiagonals(H, addr)
10
```
"""
num_offdiagonals(m::AbstractMatrix, i) = length(offdiagonals(m, i))

"""
    newadd, me = get_offdiagonal(ham, add, chosen)

Compute value `me` and new address `newadd` of a single (off-diagonal) matrix element in a
Hamiltonian `ham`. The off-diagonal element is in the same column as address `add` and is
indexed by integer index `chosen`.

# Example

```jldoctest
julia> addr = BoseFS((3, 2, 1));

julia> H = HubbardMom1D(addr);

julia> get_offdiagonal(H, addr, 3)
(BoseFS{6,3}((2, 1, 3)), 1.0)
```
"""
get_offdiagonal(m::AbstractMatrix, i, n) = offdiagonals(m, i)[n]

"""
    starting_address(h)

Return the starting address for Hamiltonian `h`. Part of the [`AbstractHamiltonian`](@ref)
interface. When called on an `AbstractMatrix` return the index of the lowest diagonal
element.

# Example

```jldoctest
julia> addr = BoseFS((3, 2, 1));

julia> H = HubbardMom1D(addr);

julia> addr == starting_address(H)
true
```
"""
starting_address(m::AbstractMatrix) = findmin(real.(diag(m)))[2]

"""
    offdiagonals(h::AbstractHamiltonian, address)

Return an iterator over reachable off-diagonal matrix elements of type
`<:AbstractOffdiagonals`. Defaults to returning `Offdiagonals(h, a)`

# See also

* [`Offdiagonals`](@ref)
* [`AbstractOffdiagonals`](@ref)

```jldoctest
julia> addr = BoseFS((3,2,1));

julia> H = HubbardReal1D(addr);

julia> h = offdiagonals(H, addr)
6-element Rimu.Hamiltonians.Offdiagonals{BoseFS{6, 3, BitString{8, 1, UInt8}}, Float64, HubbardReal1D{Float64, BoseFS{6, 3, BitString{8, 1, UInt8}}, 1.0, 1.0}}:
 (BoseFS{6,3}((2, 3, 1)), -3.0)
 (BoseFS{6,3}((2, 2, 2)), -2.449489742783178)
 (BoseFS{6,3}((3, 1, 2)), -2.0)
 (BoseFS{6,3}((4, 1, 1)), -2.8284271247461903)
 (BoseFS{6,3}((4, 2, 0)), -2.0)
 (BoseFS{6,3}((3, 3, 0)), -1.7320508075688772)
```
"""
function offdiagonals(m::AbstractMatrix, i)
    pairs = map(=>, axes(m, 1), view(m, :, i))
    return filter!(pairs) do ((k, v))
        k ≠ i && v ≠ 0
    end
end

"""
    random_offdiagonal(offdiagonals::AbstractOffdiagonals)
    random_offdiagonal(ham::AbstractHamiltonian, add)

Generate a single random excitation, i.e. choose from one of the accessible off-diagonal
elements in the column corresponding to address `add` of the Hamiltonian matrix represented
by `ham`. Alternatively, pass as argument an iterator over the accessible matrix elements.

"""
function random_offdiagonal(offdiagonals::AbstractVector)
    nl = length(offdiagonals) # check how many sites we could get_offdiagonal to
    chosen = cRand(1:nl) # choose one of them
    naddress, melem = offdiagonals[chosen]
    return naddress, 1.0/nl, melem
end

function random_offdiagonal(ham, add)
    return random_offdiagonal(offdiagonals(ham, add))
end

"""
    LOStructure(op::AbstractHamiltonian)
    LOStructure(typeof(op))

`LOStructure` speficies properties of the linear operator `op`. If a special structure is
known this can speed up calculations. Implemented structures are:

* `IsHermitian`: The operator is complex and Hermitian or real and symmetric.
* `AdjointKnown`: The operator is not Hermitian, but its [`adjoint`](@ref) is implemented.
* `AdjointUnknown`: [`adjoint`](@ref) for this operator is not implemented.

In order to define this trait for a new linear operator type, define a method for
`LOStructure(::Type{<:MyNewLOType}) = …`.

"""
abstract type LOStructure end

struct IsHermitian <: LOStructure end
struct AdjointKnown <: LOStructure end
struct AdjointUnknown <: LOStructure end

# defaults
LOStructure(op) = LOStructure(typeof(op))
LOStructure(::Type) = AdjointUnknown()
LOStructure(::AbstractMatrix) = AdjointKnown()

"""
    has_adjoint(op)

Return true if `adjoint` is defined on `op`.
"""
has_adjoint(op) = has_adjoint(LOStructure(op))
has_adjoint(::AdjointUnknown) = false
has_adjoint(::LOStructure) = true


text:
struct Cell{T<:Shape}
    row::Int
    col::Int
    connections::Set{Cell{T}}
end

function Cell(T::Type{S}, row::Int, col::Int) where {S<:Shape}
    return Cell(row, col, Set{Cell{T}}())
end

function Cell(T::Type{S}; row::Int, col::Int) where {S<:Shape}
    return Cell(T, row, col)
end

function Base.show(io::IO, cell::Cell)
    print(io, "Cell($(cell.row),$(cell.col))")
    return nothing
end

function connect!(x::Cell, y::Cell; bidirection=true)
    push!(x.connections, y)
    if bidirection
        push!(y.connections, x)
    end
end

function disconnect!(x::Cell, y::Cell, bidirection=true)
    delete!(x.connections, y)
    if bidirection
        delete!(y.connections, x)
    end
end

function is_connected(x, y; bidirection=true)
    if bidirection
        return (y ∈ x.connections) & (x ∈ y.connections)
    else
        return y ∈ x.connections
    end
end


text:
cd(@__DIR__); include("setups/path128.jl")
gr(dpi = 200)

##
f = 1.5 * 𝚽[:, 11] + 𝚽[:, 31] .* characteristic(1:32, N) +
    𝚽[:, 61] .* characteristic(1:64, N) + 0.5 * 𝚽[:, 111]

## (a)
plot(f, c = :black, grid = false, legend = false, frame = :box, lw = 1.5)
xticks!([1; 16:16:128], vcat(string(1), [string(k) for k in 16:16:128]))
plt = plot!(xlab = latexstring("x"), ylab = latexstring("f"),
    size = (600, 400), guidefontsize = 18, xlim = [1, N])
savefig(plt, "../figs/Path128_Spectrogram_f.png")

## (b)
plot(0:N-1, 𝚽' * f, c = :black, grid = false, legend = false, frame = :box, lw = 1.5)
xticks!([0; 15:16:127], vcat(string("DC"), [string(k) for k in 15:16:127]))
plt = plot!(xlab = latexstring("l"), ylab = latexstring("g"),
    size = (600, 400), guidefontsize = 18, xlim = [1, N])
savefig(plt, "../figs/Path128_Spectrogram_g.png")


text:
module AWSTools

using AWSCore
using AWSCore.Services: sts
using AWSS3
using Mocking
using Random
using Dates: Period, unix2datetime
export assume_role

get_caller_identity() = sts("GetCallerIdentity")
account_id() = (@mock get_caller_identity())["Account"]

"""
    assume_role(role_arn, [role_session_name]) -> AWSConfig

Generate a new `AWSConfig` by assuming a new role. In order to use the assumed role you need
to use this config in the various AWS calls you perform.

# Arguments
- `role_arn::AbstractString`: The ARN of the role to assume.
- `role_session_name::AbstractString`: An optional string which is the unique identifier for
    the session name.

# Keywords
- `config::AWSConfig`: The AWS configuration to use when assuming the role.
"""
function assume_role(
    role_arn::AbstractString,
    role_session_name::AbstractString=randstring(16);
    config::AWSConfig=aws_config(),
)
    function get_role_creds(role_arn, role_session_name, config)
        response = @mock sts(
            config,
            "AssumeRole",
            RoleArn=role_arn,
            RoleSessionName=role_session_name,
        )
        credentials = response["Credentials"]
        AWSCredentials(
            credentials["AccessKeyId"],
            credentials["SecretAccessKey"],
            credentials["SessionToken"],
            expiry = unix2datetime(credentials["Expiration"])
        )
    end

    renew = () -> get_role_creds(role_arn, role_session_name, config)
    creds = renew()
    creds.renew = renew

    return aws_config(creds=creds)
end

include("timeout.jl")
include("CloudFormation.jl")
include("EC2.jl")
include("ECR.jl")
include("Docker.jl")
include("S3.jl")

end  # AWSTools


text:
# Unit tests for implementations/nbodysimulator/nbodysimulator_result.jl

@testset "nbodysimulator_result.jl" begin
    particles = [
        AtomsBase.Atom(:Ar, (@SVector [7, 7, 7])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"; meta=:data),
        AtomsBase.Atom(:Ar, (@SVector [7, 7, 21])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"; hello="world"),
        AtomsBase.Atom(:Ar, (@SVector [7, 21, 7])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [7, 21, 21])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21, 7, 7])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21, 7, 21])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21, 21, 7])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au"),
        AtomsBase.Atom(:Ar, (@SVector [21, 21, 21])u"bohr", 6e-5(@SVector randn(3))u"bohr * hartree / ħ_au")
    ]
    box = [[28.0, 0.0, 0.0], [0.0, 28.0, 0.0], [0.0, 0.0, 28.0]]u"bohr"
    system = periodic_system(particles, box)

    simulator1 = NBSimulator(400, 10, t₀=1000, thermostat=NBodySimulator.AndersenThermostat(2e-4, 2e-4))
    simulator2 = NBSimulator(400, 10, t₀=1000)

    potential = InteratomicPotentials.LennardJones(1.657e-21u"J", 0.34u"nm", 0.765u"nm", [:Ar])

    result = simulate(system, simulator1, potential)
    result2 = simulate(system, simulator2, potential)

    @test get_time_range(result) == (1000:400:5000)u"ħ_au / hartree"
    @test get_num_bodies(result) == 8
    @test get_bounding_box(result) == box
    @test get_boundary_conditions(result) == [Periodic(), Periodic(), Periodic()]

    @test reference_temperature(result) == 2e-4u"hartree / k_au"
    @test ismissing(reference_temperature(result2))

    @test get_positions(result) isa AbstractVector{<:StaticVector{3,<:Unitful.Length}}
    @test get_velocities(result) isa AbstractVector{<:StaticVector{3,<:Unitful.Velocity}}
    @test get_particles(result) isa AbstractVector{<:AtomsBase.Atom}

    @test all(get_positions(result, t) isa AbstractVector{<:StaticVector{3,<:Unitful.Length}} for t ∈ 1:11)
    @test all(get_velocities(result, t) isa AbstractVector{<:StaticVector{3,<:Unitful.Velocity}} for t ∈ 1:11)
    @test all(get_particles(result, t) isa AbstractVector{<:AtomsBase.Atom} for t ∈ 1:11)

    @test all(all(all(0.0u"bohr" ≤ c < 28.0u"bohr" for c ∈ p) for p ∈ get_positions(result, t)) for t ∈ 1:11)
    @test all(get_particles(result, t)[1].meta == :data && get_particles(result, t)[2].hello == "world" for t ∈ 1:11)

    @test Atomistic.temperature(result) isa Unitful.Temperature
    @test Atomistic.kinetic_energy(result) isa Unitful.Energy
    @test Atomistic.potential_energy(result) isa Unitful.Energy

    @test all(isapprox(Atomistic.total_energy(result2, t), Atomistic.total_energy(result2), rtol=0.1) for t ∈ 1:11)
end


text:
// Calling functions which take zero parameters

int main() {
 int x = foo();
 printInt(x);
 return 0 ;

}

int foo() {
 return 10;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
module PointCloud

using LinearAlgebra
# using PyCall

import Base:
    eltype, length, minimum, take!

import ChainRulesCore:
    rrule, NoTangent

include("queue.jl")
using .PriorityQueue

include("distance.jl")
using .Distances

export embed, upper_tri
export neighborhood, geodesics, mds, isomap, scaling

# ------------------------------------------------------------------------
# globals

const ∞ = Inf

# ------------------------------------------------------------------------
# utility functions

"""
    upper_tri(x)

Returns the upper triangular portion of matrix `x`.
"""
upper_tri(x) = [ x[i.I[1], i.I[2]] for i ∈ CartesianIndices(x) if i.I[1] < i.I[2] ]
function rrule(::typeof(upper_tri), m)
	x = upper_tri(m)
	return x, (∇) -> begin
		∇m = zeros(size(m))
		n = 1
		for i ∈ CartesianIndices(∇m)
			if i.I[1] >= i.I[2]
				continue
			end

			∇m[i.I[1],i.I[2]] = ∇[n]
			∇m[i.I[2],i.I[1]] = ∇[n]

			n += 1
		end

        (NoTangent(), ∇m)
    end
end

function embed(x, dₒ; σ=0.00)
	y = if dₒ > size(x,1)
			vcat(x, zeros(dₒ-size(x,1), size(x,2)))
		elseif dₒ == size(x,1)
		    x
		else
			error("cannot embed into smaller dimension")
		end

	return y .+ σ.*randn(size(y)...)
end

# ------------------------------------------------------------------------
# types for neighborhood graph

"""
    struct Vertex{T <: Real}
        position :: Array{T}
    end

Represents a single cell within a larger point cloud.
The embedding space is normalized gene expression.
"""
struct Vertex{T <: Real}
    position :: Array{T}
end
Vertex(x) = Vertex{eltype(x)}(x)

eltype(v::Vertex{T}) where T <: Real = T

"""
    struct Edge{T <: Real}
        verts    :: Tuple{Int, Int}
        distance :: T
    end

Connects two neighboring `verts` by an edge of length `distance`.
"""
struct Edge{T <: Real}
    verts    :: Tuple{Int, Int}
    distance :: T
end
Edge(verts, distance) = Edge{typeof(distance)}(verts, distance)

eltype(e::Edge{T}) where T <: Real = T

"""
    struct Graph{T <: Real}
        verts :: Array{Vertex{T},1}
        edges :: Array{Edge{T},1}
    end

A generic graph data structure containing vertices (points in space) stored within `verts` connected by `edges`.
"""
struct Graph{T <: Real}
    verts :: Array{Vertex{T}, 1}
    edges :: Array{Edge{T}, 1}
end
Graph(verts :: Array{Vertex{T},1}) where T <: Real = Graph{T}(verts, [])

length(G :: Graph) = length(G.verts)

# ------------------------------------------------------------------------
# operations

"""
    neighborhood(x, k :: T; D=missing, accept=(d)->true) where T <: Integer

Constructs a neighborhood graph of the `k` nearest neighbor for each point of cloud `x`.
If `D` is given, it is assumed to be a dense matrix of pairwise distances.
"""
function neighborhood(x, k :: T; D=missing, accept=(d)->true) where T <: Integer
    D = ismissing(D) ? euclidean(x) : D
    G = Graph([Vertex(x[:,i]) for i ∈ 1:size(x,2)])
    for i ∈ 1:size(D,1)
        neighbor = sortperm(D[i,:])[2:end]
        append!(G.edges, [Edge((i,j), D[i,j]) for j ∈ neighbor[1:k] if accept(D[i,j])])
    end

    return G
end

"""
    neighborhood(x, k :: T; D=missing, accept=(d)->true) where T <: AbstractFloat

Constructs a neighborhood graph of all neighbors within euclidean distance `k` for each point of cloud `x`.
If `D` is given, it is assumed to be a dense matrix of pairwise distances.
"""
function neighborhood(x, k :: T; D=missing, accept=(d)->true) where T <: AbstractFloat
    D = ismissing(D) ? euclidean(x) : D
    G = Graph([Vertex(x[:,i]) for i ∈ 1:size(x,2)])
    for i ∈ 1:size(D,1)
        neighbor = first.(findall(0 .< D[:,i] .< k))
        append!(G.edges, [Edge((i,j), D[i,j]) for j ∈ neighbor if accept(D[i,j])])
    end

    return G
end

"""
    adjacency_list(G :: Graph)

Return the flattened adjacency list for graph `G`.
"""
function adjacency_list(G :: Graph)
    adj = [ Tuple{Int, Float64}[] for v ∈ 1:length(G.verts) ]
    for e ∈ G.edges
        v₁, v₂ = e.verts
        push!(adj[v₁], (v₂, e.distance))
        push!(adj[v₂], (v₁, e.distance))
    end

    return adj
end

"""
    dijkstra!(dist, adj, src)

Compute the shortest path from `src` to all other points given adjacency list `adj` and distances `dist` using Dijkstra's algorithm.
"""
function dijkstra!(dist, adj, src)
    dist      .= ∞
    dist[src]  = 0

    Q = RankedQueue((src, 0.0))
    sizehint!(Q, length(dist))

    while length(Q) > 0
        u, d₀ = take!(Q)
        for (v, d₂) ∈ adj[u]
            d₀₂ = d₀ + d₂
            if d₀₂ < dist[v]
                dist[v] = d₀₂
                if v ∈ Q
                    update!(Q, v, d₀₂)
                else
                    insert!(Q, v, d₀₂)
                end
            end
        end
    end
end

"""
    floyd_warshall(G :: Graph)

Compute the shortest path from all vertices to all other vertices within graph `G`.
"""
function floyd_warshall(G :: Graph)
    V = length(G.verts)
    D = fill(∞, (V,V))
    # remove diagonal
    for ij ∈ CartesianIndices(D)
        if ij.I[1] == ij.I[2]
            D[ij] = 0
        end
    end

    # all length 1 paths
    for e ∈ G.edges
        D[e.verts[1],e.verts[2]] = e.distance
        D[e.verts[2],e.verts[1]] = e.distance
    end

    # naive V³ paths
    for k ∈ 1:V
        for i ∈ 1:V
            for j ∈ 1:V
                if D[i,j] > D[i,k] + D[k,j]
                    D[i,j] = D[i,k] + D[k,j]
                    D[j,i] = D[i,j]
                end
            end
        end
    end

    return D
end

"""
    geodesics(G :: Graph; sparse=true)

Compute the matrix of pairwise distances, given a neighborhood graph `G`, weighted by local Euclidean distance.
If sparse is true, it will utilize Dijkstra's algorithm, individually for each point.
If sparse is false, it will utilize the Floyd Warshall algorithm.
"""
function geodesics(G :: Graph; sparse=true)
    if sparse
        adj  = adjacency_list(G)
        dist = zeros(length(G), length(G))

        Threads.@threads for v ∈ 1:length(G)
            dijkstra!(view(dist,:,v), adj, v)
        end

        return dist
    else
        return floyd_warshall(G)
    end
end

"""
    geodesics(x, k; D=missing, accept=(d)->true, sparse=true)

Compute the matrix of pairwise distances, given a pointcloud `x` and neighborhood cutoff `k`, from the resultant neighborhood graph.
If sparse is true, it will utilize Dijkstra's algorithm, individually for each point.
If sparse is false, it will utilize the Floyd Warshall algorithm.
"""
geodesics(x, k; D=missing, accept=(d)->true, sparse=true) = geodesics(neighborhood(x, k; D=D, accept=accept); sparse=sparse)

# ------------------------------------------------------------------------
# non ml dimensional reduction

"""
    mds(D², dₒ)

Computes the lowest `dₒ` components from a _Multidimensional Scaling_ analysis given pairwise squared distances `D²`.
"""
function mds(D², dₒ)
    N = size(D²,1)
    C = I - fill(1/N, (N,N))
    B = -1/2 * C*D²*C

    eig = eigen(B)
    ι   = sortperm(real.(eig.values); rev=true)

    λ = eig.values[ι]
    ν = eig.vectors[:,ι]

    return ν[:,1:dₒ] * Diagonal(sqrt.(λ[1:dₒ]))
end

"""
    isomap(x, dₒ; k=12, sparse=true)

Compute the isometric embedding of point cloud `x` into `dₒ` dimensions.
Geodesic distances between all points are estimated by utilizing the shortest path defined by the neighborhood graph.
The embedding is computed from a multidimensional scaling analysis on the resultant geodesics.
"""
function isomap(x, dₒ; k=12, sparse=true)
    G = neighborhood(x, k)
    D = geodesics(G; sparse=sparse)
    return mds(D.^2, dₒ)
end

# ------------------------------------------------------------------------
# dimension estimation

"""
    scaling(D, N)

Estimate the Hausdorff dimension by computing how the number of points contained within balls scales with varying radius.
"""
function scaling(D, N)
	Rₘᵢₙ = minimum(D[D .> 0])
	Rₘₐₓ = maximum(D)
	Rs   = range(Rₘᵢₙ,Rₘₐₓ,length=N)

	ϕ = zeros(size(D,1), N)
	for (i,R) ∈ enumerate(Rs)
		ϕ[:,i] = sum(D .<= R, dims=1)
	end

	return ϕ, Rs
end

# ------------------------------------------------------------------------
# tests

using Statistics

function test()
    r, ξ = sphere(2000)
    D = spherical_distance(r)

    r = embed(r, size(r,2); σ=0.1)
    ks = collect(4:2:50)
    ρ  = zeros(length(ks))

    for (i,k) ∈ enumerate(ks)
        D̂ = geodesics(r, k)
        ρ[i] = cor(upper_tri(D̂), upper_tri(D))
    end

    return ks, ρ
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
#===============================================================================
  Anderson acceleration is a fixed point iteration acceleration method. A fixed
  point problem is one of finding a solution to g(x) = x or alternatively to
  find a solution to F(x) = g(x) - x = 0. To allow for easy switching between
  methods, we write the code in the F(x) = 0 form, but it really is of no impor-
  tance. This makes it a bit simpler to keep a consistent naming scheme as well,
  since convergence is measured in F(x) not g(x) directly.
===============================================================================#
struct Anderson{Ti, Tb, Td}
  delay::Ti
  memory::Ti
  beta::Tb
  droptol::Td
end
Anderson() = Anderson(0, 5, nothing, nothing)
function vv_shift!(G)
  for i = 1:length(G)-1
    G[i] = G[i+1]
  end
end

                     # args
function fixedpoint!(g, x,
                     anderson::Anderson;
                     # kwargs
                     Gx = similar(x),
                     Fx = similar(x),
                     f_abstol=sqrt(eps(eltype(x))),
                     maxiter=100)
  #==============================================================
    Do initial function iterations; default is a delay of 0,
    but we always do at least one evaluation of G, to set up AA.

    Notation: AA solves g(x) - x = 0 or F(x) = 0
  ==============================================================#
  Gx  = g(Gx, x)
  Fx .= Gx .- x
  ρF0 = norm(Fx, Inf)
  ρ2F0 = norm(Fx, 2)
  x  .= Gx
  for i = 1:anderson.delay
    Gx = g(Gx, x)
    Fx .= Gx .- x
    x .= Gx
    finite_check = isallfinite(x)
    if norm(Fx) < f_abstol || !finite_check
      return (x=x, Fx=Fx, acc_iter=0, finite=finite_check)
    end
  end

  #==============================================================
    If we got this far, then the delay was not enough to con-
    verge. However, we now hope to have moved to a region where
    everything is well-behaved, and we start the acceleration.
  ==============================================================#

  n = length(x)
  memory = min(n, anderson.memory)

  Q = x*x[1:memory]'
  R = x[1:memory]*x[1:memory]'

  #==============================================================
    Start Anderson Acceleration. We use QR updates to add new
    successive changes in G to the system, and once the memory
    is exhausted, we use QR downdates to forget the oldest chan-
    ges we have stored.
  ==============================================================#
  effective_memory = 0
  beta = nothing

  G = [copy(x) for i=1:memory]
  Δg = copy(Gx)
  Δf = copy(Fx)

  Gold = copy(Gx)
  for k = 1:maxiter
    Fold = copy(Fx)
    Gx = g(x, Gx)
	  Fx .= Gx .- x
    x .= Gx

    # is this actually needed? I think we can avoid these
    @. Δg = Gx - Gold
    @. Δf = Fx - Fold

    Gold .= Gx
    Fold .= Fx

    effective_memory += 1
    memory
    # if we've exhausted the memory, downdate
    if effective_memory > memory
      vv_shift!(G)
      qrdelete!(Q, R, memory)
      effective_memory -= 1
    end

    # Add the latest change to G
    G[effective_memory] .= Δg

    # QR update
  	qradd!(Q, R, vec(Δf), effective_memory)

  	# Create views for the system depending on the effective memory counter
  	Qv = view(Q, :, 1:effective_memory)
  	Rv = UpperTriangular(view(R, 1:effective_memory, 1:effective_memory))

  	# check condition number
    if !isa(anderson.droptol, Nothing)
        while cond(Rv) > anderson.droptol && effective_memory > 1
            qrdelete!(Q, R, effective_memory)

            effective_memory -= 1
            Qv = view(Q, :, 1:effective_memory)
            Rv = UpperTriangular(view(R, 1:effective_memory, 1:effective_memory))
        end
    end

    # solve least squares problem
    γv = zeros(effective_memory) #view(γv, 1:m_eff)
    ldiv!(Rv, mul!(γv, Qv', vec(Fx)))

    # update next iterate
    for i in 1:effective_memory
        @. x -= γv[i] * G[i]
    end
    if !isa(beta, Nothing)
        x .= x .- (1 .- beta).*(Fx .- Qv*Rv*γv) # this is suboptimal!
    end

    finite_check = isallfinite(x)
    if norm(Fx) < f_abstol || !finite_check
      return (x=x, Fx=Fx, acc_iter=0, finite=finite_check)
    end
  end
  ConvergenceInfo(method, (solution=x, best_residual=Fx, ρF0=ρF0, ρ2F0=ρ2F0, iter=iter, time=time()-t0), options)
  Gx, Gx, x
end


text:
function go()
    global i_counter = 0
    global l_counter = 0

    println("Test JULIA:")
    for i_loop1 in 0:9
        println("i_loop1: ",i_loop1)
        for i_loop2 in 0:31999
            for i_loop3 in 0:31999
                i_counter = i_counter + 1
                l_counter = l_counter + 1
                if i_counter == 50
                    i_counter = 0
                end
            end
    end
    end
    println("i_counter: ", i_counter)
    println("l_counter: ", l_counter)
end
go()

max_stars_repo_name: RaSla/lang-test

max_issues_repo_name: RaSla/lang-test

max_forks_repo_name: RaSla/lang-test


text:
using Test
using BCVUMPS
using BCVUMPS:qrpos,lqpos,leftorth,rightorth,leftenv,FLmap,rightenv,FRmap,ACenv,ACmap,Cenv,Cmap,ACCtoALAR,error,obs2x2FL,obs2x2FR,bigleftenv,BgFLmap,bigrightenv,BgFRmap
using LinearAlgebra
using OMEinsum
using Random

@testset "qr" begin
    A = rand(ComplexF64,4,4)
    Q, R = qrpos(A)
    @test Array(Q*R) ≈ Array(A)
    @test all(real.(diag(R)) .> 0)
    @test all(imag.(diag(R)) .≈ 0)
end

@testset "lq" begin
    A = rand(ComplexF64,4,4)
    L, Q = lqpos(A)
    @test Array(L*Q) ≈ Array(A)
    @test all(real.(diag(L)) .> 0)
    @test all(imag.(diag(L)) .≈ 0)
end

@testset "leftorth and rightorth" for Ni = [2,3], Nj = [2,3]
    Random.seed!(50)
    D, d = 5, 2
    A = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = rand(D, d, D)
    end
    AL, L, λ = leftorth(A)
    R, AR, λ = rightorth(A)

    for j = 1:Nj,i = 1:Ni
        M = ein"cda,cdb -> ab"(AL[i,j],conj(AL[i,j]))
        @test (Array(M) ≈ I(D))

        LA = reshape(L[i,j] * reshape(A[i,j], D, d*D), d*D, D)
        ALL = reshape(AL[i,j], d*D, D) * L[i,j] * λ[i,j]
        @test (Array(ALL) ≈ Array(LA))

        M = ein"acd,bcd -> ab"(AR[i,j],conj(AR[i,j]))
        @test (Array(M) ≈ I(D))

        AxR = reshape(reshape(A[i,j], d*D, D)*R[i,j], D, d*D)
        RAR = R[i,j] * reshape(AR[i,j], D, d*D) * λ[i,j]
        @test (Array(RAR) ≈ Array(AxR))
    end
end

@testset "leftenv and rightenv" for Ni = [2,3], Nj = [2,3]
    Random.seed!(50)
    D, d = 5, 2
    A = Array{Array,2}(undef, Ni, Nj)
    M = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = rand(D, d, D)
        M[i,j] = rand(2,2,2,2)
    end

    AL, = leftorth(A)
    λL,FL = leftenv(AL, M)
    _, AR, = rightorth(A)
    λR,FR = rightenv(AR, M)

    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        @test λL[i,j] * FL[i,j] ≈ FLmap(AL[i,:], AL[ir,:], M[i,:], FL[i,j], j)
        @test λR[i,j] * FR[i,j] ≈ FRmap(AR[i,:], AR[ir,:], M[i,:], FR[i,j], j)
    end
end

@testset "ACenv and Cenv" for Ni = [2,3], Nj = [2,3]
    Random.seed!(50)
    D, d = 5, 2
    A = Array{Array,2}(undef, Ni, Nj)
    M = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = rand(D, d, D)
        M[i,j] = rand(2,2,2,2)
    end

    AL, L = leftorth(A)
    λL,FL = leftenv(AL, M)
    R, AR, = rightorth(A)
    λR,FR = rightenv(AR, M)

    C = Array{Array,2}(undef, Ni,Nj)
    AC = Array{Array,2}(undef, Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        jr = j + 1 - (j+1>Nj) * Nj
        C[i,j] = L[i,j] * R[i,jr]
        AC[i,j] = ein"asc,cb -> asb"(AL[i,j],C[i,j])
    end
    λAC, AC = ACenv(AC, FL, M, FR)
    λC, C = Cenv(C, FL, FR)
    for j = 1:Nj, i = 1:Ni
        jr = j + 1 - Nj * (j==Nj)
        @test λAC[i,j] * AC[i,j] ≈ ACmap(AC[i,j], FL[:,j], FR[:,j], M[:,j], i)
        @test λC[i,j] * C[i,j] ≈ Cmap(C[i,j], FL[:,jr], FR[:,j], i)
    end
end

@testset "bcvumps unit test" for Ni = [2,3], Nj = [2,3]
    Random.seed!(50)
    D, d = 5, 2
    A = Array{Array,2}(undef, Ni, Nj)
    M = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = rand(D, d, D)
        M[i,j] = rand(2,2,2,2)
    end

    AL, L = leftorth(A)
    λL,FL = leftenv(AL, M)
    R, AR, = rightorth(A)
    λR,FR = rightenv(AR, M)

    C = Array{Array,2}(undef, Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        jr = j + 1 - (j+1>Nj) * Nj
        C[i,j] = L[i,j] * R[i,jr]
    end

    AL, C, AR = ACCtoALAR(AL, C, AR, M, FL, FR)
    err = error(AL,C,FL,M,FR)
    @test err !== nothing
end

@testset "obsleftenv and obsrightenv" for Ni = [2], Nj = [2]
    Random.seed!(50)
    D, d = 5, 2
    A = Array{Array,2}(undef, Ni, Nj)
    M = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = rand(D, d, D)
        M[i,j] = rand(2,2,2,2)
    end

    AL, = leftorth(A)
    λL,FL = obs2x2FL(AL, M)
    _, AR, = rightorth(A)
    λR,FR = obs2x2FR(AR, M)

    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        @test λL[i,j] * FL[i,j] ≈ FLmap(AL[i,:], AL[i,:], M[i,:], FL[i,j], j)
        @test λR[i,j] * FR[i,j] ≈ FRmap(AR[i,:], AR[i,:], M[i,:], FR[i,j], j)
    end
end

@testset "bigleftenv and bigrightenv" for Ni = [2,3], Nj = [2,3]
    Random.seed!(50)
    D, d = 5, 2
    A = Array{Array,2}(undef, Ni, Nj)
    M = Array{Array,2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = rand(D, d, D)
        M[i,j] = rand(2,2,2,2)
    end

    AL, = leftorth(A)
    λL,BgFL = bigleftenv(AL, M)
    _, AR, = rightorth(A)
    λR,BgFR = bigrightenv(AR, M)

    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        @test λL[i,j] * BgFL[i,j] ≈ BgFLmap(AL[i,:], AL[ir,:], M[i,:], M[ir,:], BgFL[i,j], j)
        @test λR[i,j] * BgFR[i,j] ≈ BgFRmap(AR[i,:], AR[ir,:], M[i,:], M[ir,:], BgFR[i,j], j)
    end
end


text:
using Diderot: Arc, Node, RestrictLowDistance
using Diderot.Knapsack: Instance, DecreasingWeight, relax

@testset "model methods" begin
    instance = Instance([4.0, 3.0, 2.0], [3, 2, 2], 4)

    @test Diderot.initial_state(instance) == 4
    @test Diderot.value_type(instance) == Float64
    @test Diderot.domain_type(instance) == Bool

    diagram = Diderot.Diagram(instance)

    @test Diderot.next_variable(instance, diagram, Diderot.InOrder()) == 1
    @test Diderot.next_variable(instance, diagram, DecreasingWeight()) == 1

    @test Diderot.transitions(instance, 2, 1) ==
        Dict(Arc(2, false, 0.0) => 2)
    @test Diderot.transitions(instance, 2, 2) ==
        Dict(Arc(2, false, 0.0) => 2,
             Arc(2, true,  3.0) => 0)

    instance2 = Instance([4.0, 3.0, 2.0], [2, 3, 2], 4)
    @test Diderot.next_variable(instance2, diagram, DecreasingWeight()) == 2
end

@testset "top-down decision diagrams" begin
    N = Node{Int,Bool,Float64}

    instance = Instance([4.0, 3.0, 2.0], [3, 2, 2], 4)
    diagram = Diderot.Diagram(instance)
    Diderot.top_down!(diagram, instance)

    # layers and nodes
    @test length(diagram.layers) == 4

    @test length(diagram.layers[1]) == 1 # root
    @test diagram.layers[1][4] == N(0.0, nothing)

    @test length(diagram.layers[2]) == 2 # 4, 1
    @test diagram.layers[2][4] == N(0.0, Arc(4, false, 0.0))
    @test diagram.layers[2][1] == N(4.0, Arc(4, true, 4.0))

    @test length(diagram.layers[3]) == 3 # 4, 2, 1
    @test diagram.layers[3][4] == N(0.0, Arc(4, false, 0.0))
    @test diagram.layers[3][2] == N(3.0, Arc(4, true, 3.0))
    @test diagram.layers[3][1] == N(4.0, Arc(1, false, 0.0))

    @test length(diagram.layers[4]) == 1 # terminal
    @test diagram.layers[4][0] == N(5.0, Arc(2, true, 2.0))
end

@testset "longest path" begin
    instance = Instance([4.0, 3.0, 2.0], [3, 2, 2], 4)
    diagram = Diderot.Diagram(instance)
    Diderot.top_down!(diagram, instance)
    solution = Diderot.longest_path(diagram)
    @test solution.decisions == [false, true, true]
    @test solution.objective ≈ 5.0

    instance2 = Instance([3.0, 4.0, 2.0], [2, 3, 2], 4)
    diagram2 = Diderot.Diagram(instance2)
    Diderot.top_down!(diagram2, instance2, variable_order=DecreasingWeight())
    solution2 = Diderot.longest_path(diagram2)
    @test solution2.decisions == [true, false, true]
    @test solution2.objective ≈ 5.0
end

@testset "restriction" begin
    instance = Instance([4.0, 3.0, 2.0], [3, 2, 2], 4)

    @testset "width 1" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=RestrictLowDistance(1))
        @test length(diagram.layers) == 4
        @test all(l -> length(l) == 1, diagram.layers)

        solution = Diderot.longest_path(diagram)
        @test solution.decisions == [true, false, false]
        @test solution.objective ≈ 4.0
    end

    @testset "width 2" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=RestrictLowDistance(2))
        @test length(diagram.layers) == 4
        @test all(l -> length(l) in (1, 2), diagram.layers)

        solution = Diderot.longest_path(diagram)
        @test solution.decisions == [false, true, true]
        @test solution.objective ≈ 5.0 # optimum!
    end
end

@testset "relaxation" begin
    instance = Instance([4.0, 3.0, 2.0], [3, 2, 2], 4)

    @testset "width 1" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=relax(1))
        @test length(diagram.layers) == 4
        @test all(l -> length(l) == 1, diagram.layers)

        solution = Diderot.longest_path(diagram)
        @test solution.decisions == [true, true, true]
        @test solution.objective ≈ 9.0
    end

    @testset "width 2" begin
        diagram = Diderot.Diagram(instance)
        Diderot.top_down!(diagram, instance, processing=relax(2))
        @test length(diagram.layers) == 4
        @test all(l -> length(l) in (1, 2), diagram.layers)

        solution = Diderot.longest_path(diagram)
        @test solution.decisions == [true, false, true]
        @test solution.objective ≈ 6.0
    end
end

@testset "branch and bound" begin
    instance = Instance([4.0, 3.0, 2.0], [3, 2, 2], 4)

    # Relaxation needs at least width 2 for branching.

    @testset "width 1" begin
        solution = Diderot.branch_and_bound(
            instance, restrict=RestrictLowDistance(1), relax=relax(2))
        @test solution.decisions == [false, true, true]
        @test solution.objective ≈ 5.0
    end

    @testset "width 2" begin
        solution = Diderot.branch_and_bound(
            instance, restrict=RestrictLowDistance(2), relax=relax(2))
        @test solution.decisions == [false, true, true]
        @test solution.objective ≈ 5.0
    end

    @testset "width 3" begin
        solution = Diderot.branch_and_bound(
            instance, restrict=RestrictLowDistance(3), relax=relax(3))
        @test solution.decisions == [false, true, true]
        @test solution.objective ≈ 5.0
    end

    @testset "n 5, width 2" begin
        instance2 = Instance([5, 3, 2, 7, 4], [2, 8, 4, 2, 5], 10)
        solution = Diderot.branch_and_bound(
            instance2, restrict=RestrictLowDistance(2), relax=relax(2))
        @test solution.decisions == [1, 0, 0, 1, 1]
        @test solution.objective ≈ 16.0
    end
end


text:
const _ComponentnameReferences = Set{Tuple{UUIDs.UUID, String}}

struct _TimeSeriesRecord
    component_names::_ComponentnameReferences
    ts::TimeSeriesData
end

function _TimeSeriesRecord(component_uuid, name, ts)
    record = _TimeSeriesRecord(_ComponentnameReferences(), ts)
    push!(record.component_names, (component_uuid, name))
    return record
end

"""
Stores all time series data in memory.
"""
struct InMemoryTimeSeriesStorage <: TimeSeriesStorage
    data::Dict{UUIDs.UUID, _TimeSeriesRecord}
end

function InMemoryTimeSeriesStorage()
    storage = InMemoryTimeSeriesStorage(Dict{UUIDs.UUID, _TimeSeriesRecord}())
    @debug "Created InMemoryTimeSeriesStorage" _group = LOG_GROUP_TIME_SERIES
    return storage
end

"""
Constructs InMemoryTimeSeriesStorage from an instance of Hdf5TimeSeriesStorage.
"""
function InMemoryTimeSeriesStorage(hdf5_storage::Hdf5TimeSeriesStorage)
    storage = InMemoryTimeSeriesStorage()
    for (component, name, time_series) in iterate_time_series(hdf5_storage)
        serialize_time_series!(storage, component, name, time_series)
    end

    return storage
end

check_read_only(storage::InMemoryTimeSeriesStorage) = nothing

get_compression_settings(storage::InMemoryTimeSeriesStorage) =
    CompressionSettings(enabled = false)

is_read_only(storage::InMemoryTimeSeriesStorage) = false

function serialize_time_series!(
    storage::InMemoryTimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
    ts::TimeSeriesData,
)
    uuid = get_uuid(ts)
    if !haskey(storage.data, uuid)
        @debug "Create new time series entry." _group = LOG_GROUP_TIME_SERIES uuid component_uuid name
        storage.data[uuid] = _TimeSeriesRecord(component_uuid, name, ts)
    else
        add_time_series_reference!(storage, component_uuid, name, uuid)
    end

    return
end

function add_time_series_reference!(
    storage::InMemoryTimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
    ts_uuid::UUIDs.UUID,
)
    @debug "Add reference to existing time series entry." _group = LOG_GROUP_TIME_SERIES ts_uuid component_uuid name
    record = storage.data[ts_uuid]
    key = (component_uuid, name)

    # It's possible that this is overly restrictive, but as of now there is not a good
    # reason for a caller to add a reference multiple times. This should be a bug.
    @assert !in(key, record.component_names) "There is already a reference to $key for time series $ts_uuid"

    push!(record.component_names, key)
end

function remove_time_series!(
    storage::InMemoryTimeSeriesStorage,
    uuid::UUIDs.UUID,
    component_uuid::UUIDs.UUID,
    name::AbstractString,
)
    if !haskey(storage.data, uuid)
        throw(ArgumentError("$uuid is not stored"))
    end

    record = storage.data[uuid]
    component_name = (component_uuid, name)
    if !(component_name in record.component_names)
        throw(ArgumentError("$component_name wasn't stored for $uuid"))
    end

    pop!(record.component_names, component_name)
    @debug "Removed $component_name from $uuid." _group = LOG_GROUP_TIME_SERIES

    if isempty(record.component_names)
        @debug "$uuid has no more references; delete it." _group = LOG_GROUP_TIME_SERIES
        pop!(storage.data, uuid)
    end
end

function deserialize_time_series(
    ::Type{T},
    storage::InMemoryTimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: StaticTimeSeries}
    uuid = get_time_series_uuid(ts_metadata)
    if !haskey(storage.data, uuid)
        throw(ArgumentError("$uuid is not stored"))
    end

    ts = storage.data[uuid].ts
    total_rows = length(ts_metadata)
    if rows.start == 1 && length(rows) == total_rows
        # No memory allocation
        return ts
    end

    # TimeArray doesn't support @view
    return T(ts, get_data(ts)[rows])
end

function deserialize_time_series(
    ::Type{T},
    storage::InMemoryTimeSeriesStorage,
    ts_metadata::TimeSeriesMetadata,
    rows::UnitRange,
    columns::UnitRange,
) where {T <: TimeSeriesData}
    uuid = get_time_series_uuid(ts_metadata)
    if !haskey(storage.data, uuid)
        throw(ArgumentError("$uuid is not stored"))
    end

    ts = storage.data[uuid].ts
    if ts isa SingleTimeSeries
        return deserialize_deterministic_from_single_time_series(
            storage,
            ts_metadata,
            rows,
            columns,
            length(ts),
        )
    end

    total_rows = length(ts_metadata)
    total_columns = get_count(ts_metadata)
    if length(rows) == total_rows && length(columns) == total_columns
        return ts
    end

    full_data = get_data(ts)
    initial_timestamp = get_initial_timestamp(ts)
    resolution = get_resolution(ts)
    interval = get_interval(ts)
    start_time = initial_timestamp + interval * (columns.start - 1)
    data = SortedDict{Dates.DateTime, eltype(typeof(full_data)).parameters[2]}()
    for initial_time in range(start_time; step = interval, length = length(columns))
        if rows.start == 1
            it = initial_time
        else
            it = initial_time + (rows.start - 1) * resolution
        end
        data[it] = @view full_data[initial_time][rows]
    end

    if T <: AbstractDeterministic
        return Deterministic(ts, data)
    else
        return T(ts, data)
    end
end

function clear_time_series!(storage::InMemoryTimeSeriesStorage)
    empty!(storage.data)
    @info "Cleared all time series."
end

function get_num_time_series(storage::InMemoryTimeSeriesStorage)
    return length(storage.data)
end

function convert_to_hdf5(storage::InMemoryTimeSeriesStorage, filename::AbstractString)
    create_file = true
    hdf5_storage = Hdf5TimeSeriesStorage(create_file; filename = filename)
    for record in values(storage.data)
        for pair in record.component_names
            serialize_time_series!(hdf5_storage, pair[1], pair[2], record.ts)
        end
    end
end

function compare_values(
    x::InMemoryTimeSeriesStorage,
    y::InMemoryTimeSeriesStorage;
    compare_uuids = false,
)
    keys_x = sort!(collect(keys(x.data)))
    keys_y = sort!(collect(keys(y.data)))
    if keys_x != keys_y
        @error "keys don't match" keys_x keys_y
        return false
    end

    for key in keys_x
        record_x = x.data[key]
        record_y = y.data[key]
        if record_x.component_names != record_y.component_names
            @error "component_names don't match" record_x.component_names record_y.component_names
            return false
        end
        if TimeSeries.timestamp(record_x.ts.data) != TimeSeries.timestamp(record_y.ts.data)
            @error "timestamps don't match" record_x record_y
            return false
        end
        if TimeSeries.values(record_x.ts.data) != TimeSeries.values(record_y.ts.data)
            @error "values don't match" record_x record_y
            return false
        end
    end

    return true
end


text:
__precompile__()

module SimulatedAnnealing

    export

    program,
    CoolingScheduler,

    StabilityCounter,

    State,
    StatePair,
    ConcreteState,

    update,
    get_acceptance,
    MH_init,
    MH_body!,
    MH_algorithm,

    boltzmann_factor,
    acceptance,
    SA_algorithm

    include("cooling.jl")
    include("stable.jl")
    include("state.jl")
    include("mh.jl")
    include("sa.jl")

end # module


text:
include("../src/gridded_module.jl")

module gtest
using BenchmarkTools
using ..gridded
using ..util
using PyPlot
PyPlot.pygui(true)

function test_static_curve()

    x = [0.0,  1.0, 2.0, 3.0, 4.0]
    y = [10.0, 20,  30,  40,  50]

    z = 2 * copy(x)

    sc = gridded.StaticCurve(x, y)
    # sc.x = [10] # fails to bind a new curve to that -> good
    sc.x[1] = 10 # still, elements can still be changed

    mc = gridded.MutableCurve(x, y) # the vectors itself can be changed to other vectors == the pointers
    # mc.x[1] = 20
    mc.x = z;

    println( length(mc) )

    PyPlot.figure()
    PyPlot.plot(sc.x, sc.y)

    PyPlot.figure()
    PyPlot.plot(mc.x, mc.y)

end

function test_push()

    x = [0.0,  1.0, 2.0, 3.0, 4.0]
    y = [10.0, 20,  30,  40,  50]

    x_new = 10
    y_new = 50

    sc = gridded.StaticCurve(x, y)
    mc = gridded.MutableCurve(x, y)

    show(sc)
    push!(sc, x_new, y_new)
    show(sc)

    show(mc)
    push!(mc, x_new, y_new)
    show(mc)

end

function test_inv()
    x = [0.0,  1.0, 2.0, 3.0, 4.0]
    y = [10.0, 20,  30,  40,  50]

    mc = gridded.MutableCurve(x, y)

    mc_inv = gridded.inverse(mc)
    show(mc_inv)

end

function test_linear_interpolation()

    crv = gridded.ramp_curve(gridded.MutableCurve{Float64}, 0.0, 0.1, 5.0)
    gridded.map_y!(crv, sin)

    xq = pi/2 # -1.0 # 6.0 # 3.14
    # (yq, right_idx) = gridded.linear_search_interpolate(crv, xq)

    (yq, right_idx) = gridded.binary_search_interpolate(crv, xq)

    println(yq)
    println(right_idx)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(crv.x, crv.y, marker = :.)
    PyPlot.plot(xq, yq, marker = :o)

end

function benchmark_fcn_linear_interpolation(crv, x_vec)
    # x_query = x_vec
    x_query = sort(x_vec)
    left_idx = 1
    for kk = 1:length(x_query)
        xq = x_query[kk]
        (yq, left_idx) = gridded.linear_search_interpolate(crv, xq, init_index = left_idx)
    end
end

function benchmark_fcn_binary_interpolation(crv, x_vec)
    # x_query = x_vec
    x_query = sort(x_vec)
    left_idx = 1
    for kk = 1:length(x_query)
        xq = x_query[kk]
        (yq, left_idx) = gridded.binary_search_interpolate(crv, xq, left_init_index = left_idx)
    end
end

function benchmark_linear_interpolation(crv, x_vec)
    b = @benchmark benchmark_fcn_linear_interpolation($crv, $x_vec)

    println("=== benchmark_linear_interpolation ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_binary_interpolation(crv, x_vec)
    b = @benchmark benchmark_fcn_binary_interpolation($crv, $x_vec)

    println("=== benchmark_binary_interpolation ===")
    show(
            stdout,
            MIME("text/plain"),
            b )
    println("\n\n")
end

function benchmark_all_interpolations()

    crv = gridded.ramp_curve(gridded.MutableCurve{Float64}, 0.0, 0.001, 5.0)
    crv.y = gridded.map_y(crv, sin)

    n_vec = [50, 500, 5_000, 50_000]
    for n in n_vec
        # n = 5_0
        a = -1.0
        b = 6.0
        x_vec = (b - a) .* rand(n) .+ a # gridded.ramp(-1.0, 0.003, 6.0)

        println(" ========================================== ")
        println("N: $(n)")
        println(" ========================================== ")
        benchmark_linear_interpolation(crv, x_vec)
        benchmark_binary_interpolation(crv, x_vec)
    end

end

# =========================================================================== #
# Map
# =========================================================================== #

function test_map()

    x = [10.0, 20, 30, 40, 50]
    y = [100.0, 200, 300]
    z = [
        1.0 2   3;
        4   5   6;
        7   8   9;
        10  11  12;
        13  14  15;
        ]

    map = gridded.GriddedMap(x, y, z)
    show(map)

    gridded.map_z!(map, q -> q^2)
    show(map)

    (X, Y) = util.mesh_grid(map.x, map.y)

    PyPlot.figure()
    PyPlot.surf(X, Y, map.z)

end

function test_map_interp()
    x = [10.0, 20, 30, 40, 50]
    y = [100.0, 200, 300]
    z = [
        1.0 2   3;
        4   5   6;
        7   8   9;
        10  11  12;
        13  14  15;
        ]

    map = gridded.GriddedMap(x, y, z)

    xq = 15.0
    yq = 250.0
    (zq, x_idx, y_idx) = gridded.linear_search_interpolate(map, xq, yq)

    println(zq)

    (X, Y) = util.mesh_grid(map.x, map.y)

    PyPlot.figure()
    PyPlot.surf(X, Y, map.z, alpha = 0.5)
    PyPlot.plot3D(xq, yq, zq, marker = :., markersize = 20, color = :k)

end

function test_peaks()

    x_vec = gridded.linspace(-3.0, 3, 10)
    y_vec = gridded.linspace(-3.0, 3, 10)
    map = gridded.peaks(x_vec, y_vec)

    xq = 1.1
    yq = 2.1
    # (zq, x_idx, y_idx) = gridded.linear_search_interpolate(map, xq, yq)
    (zq, x_idx, y_idx) = gridded.binary_search_interpolate(map, xq, yq)

    println(zq)
    println(x_idx)
    println(y_idx)

    (X, Y) = util.mesh_grid(map.x, map.y)

    PyPlot.figure()
    PyPlot.surf(X, Y, map.z, alpha = 0.9)
    PyPlot.plot3D(xq, yq, zq, marker = :., markersize = 20, color = :k)

end

end # gtest

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
#
# Copyright (c) 2015, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

"""
Apply zero padding on all sides of an image. This keeps the image relatively
well centered even after zero padding, which is useful if you want to
avoiding a large tilt in the phase of the function after applying the
Fourier transform (i.e. via the Fourier shift property).  Such a tilt is
generally undesirable when using the FFT to model the Fourier transforming
property of an optical system.
"""
function fftpad(u1::AbstractArray{T,1}) where {T}
    padM = Int(ceil(size(u1, 1) / 2.0))
    temp = zeros(eltype(u1), size(u1, 1) * 2)
    temp[padM+1:padM+size(u1, 1)] .= u1
    return temp
end

function fftpad(u1::AbstractArray{T,2}) where {T}
    padM = Int(ceil(size(u1, 1) / 2.0))
    padN = Int(ceil(size(u1, 2) / 2.0))
    temp = zeros(eltype(u1), size(u1, 1) * 2, size(u1, 2) * 2)
    temp[padM+1:padM+size(u1, 1), padN+1:padN+size(u1, 2)] .= u1 # Zero pad the input
    return temp
end

function fftpad(u1::AbstractArray{T,3}) where {T}
    padM = Int(ceil(size(u1, 1) / 2.0))
    padN = Int(ceil(size(u1, 2) / 2.0))
    padK = Int(ceil(size(u1, 3) / 2.0))
    temp = zeros(eltype(u1), size(u1, 1) * 2, size(u1, 2) * 2, size(u1, 3) * 2)
    temp[padM+1:padM+size(u1, 1), padN+1:padN+size(u1, 2), padK+1:padK+size(u1, 3)] .= u1 # Zero pad the input
    return temp
end

"""
Remove symmetric zero padding.
"""
function fftunpad(u1::AbstractArray{T,1}) where {T}
    M = size(u1, 1)
    origM = div(M, 2)
    padM = Int(ceil(origM / 2.0))
    return view(u1, padM+1:padM+origM)
end

function fftunpad(u1::AbstractArray{T,2}) where {T}
    (M, N) = size(u1)
    origM = div(M, 2)
    origN = div(N, 2)
    padM = Int(ceil(origM / 2.0))
    padN = Int(ceil(origN / 2.0))
    return view(u1, padM+1:padM+origM, padN+1:padN+origN)                 # remove zero padding
end

function fftunpad(u1::AbstractArray{T,3}) where {T}
    (M, N, K) = size(u1)
    origM = div(M, 2)
    origN = div(N, 2)
    origK = div(K, 2)
    padM = Int(ceil(origM / 2.0))
    padN = Int(ceil(origN / 2.0))
    padK = Int(ceil(origK / 2.0))
    return view(u1, padM+1:padM+origM, padN+1:padN+origN, padK+1:padK+origK)                 # remove zero padding
end

function fftunpad2(u1::AbstractArray{T,3}) where {T}
    (M, N, K) = size(u1)
    origM = div(M, 2)
    origN = div(N, 2)
    padM = Int(ceil(origM / 2.0))
    padN = Int(ceil(origN / 2.0))
    return view(u1, padM+1:padM+origM, padN+1:padN+origN, 1:K)                 # remove zero padding
end

#-------------------------------------------------------------------------------------
# Centered, optimal size fft padding routines

"""
Zero pads the FFT up to at least double the size of the original array, or larger
if a nearby larger dimension will have only small prime least common denominators.
Choosing the dimensions in this way can speed up the FFT algorithm in some cases.
"""
function fftpad_optimal(u1::AbstractArray{T,2}) where {T}
    optimal_shape =
        (nextprod([2, 3, 5, 7], size(u1, 1) * 2), nextprod([2, 3, 5, 7], size(u1, 2) * 2))
    padM = Int(ceil(optimal_shape[1] / 2.0 - size(u1, 1) / 2.0))
    padN = Int(ceil(optimal_shape[2] / 2.0 - size(u1, 2) / 2.0))
    temp = zeros(eltype(u1), optimal_shape)
    temp[padM+2:padM+size(u1, 1)+1, padN+2:padN+size(u1, 2)+1] .= u1 # Zero pad the input in a fashion consistent with direct convolution
    return temp
end

"""
Removes optimal zero padding.  You must supply the original (non-zero padded)
array shape.
"""
function fftunpad_optimal(u1::AbstractArray{T,2}, orig_size::Tuple{Int64,Int64}) where {T}
    optimal_shape =
        (nextprod([2, 3, 5, 7], orig_size[1] * 2), nextprod([2, 3, 5, 7], orig_size[2] * 2))
    origM, origN = orig_size
    padM = Int(ceil(optimal_shape[1] / 2.0 - orig_size[1] / 2.0))
    padN = Int(ceil(optimal_shape[2] / 2.0 - orig_size[2] / 2.0))
    return view(u1, padM+2:padM+origM+1, padN+2:padN+origN+1)                 # remove zero padding
end


text:
# This file generates the solution used in method of manufactured solutions
using LinearAlgebra, SymPy, Printf

@syms x y z t real=true
μ = 1 // 100
γ = 7 // 5

output = open("mms_solution_generated.jl", "w")

@printf output "const γ_exact = %s\n" γ
@printf output "const μ_exact = %s\n" μ

for dim = 2:3
  if dim == 3
    ρ = cos(π * t) * sin(π * x) * cos(π * y) * cos(π * z) + 3
    U = cos(π * t) * ρ * sin(π * x) * cos(π * y) * cos(π * z)
    V = cos(π * t) * ρ * sin(π * x) * cos(π * y) * cos(π * z)
    W = cos(π * t) * ρ * sin(π * x) * cos(π * y) * sin(π * z)
    E = cos(π * t) * sin(π * x) * cos(π * y) * cos(π * z) + 100
  else
    ρ = cos(π * t) * sin(π * x) * cos(π * y) + 3
    U = cos(π * t) * ρ * sin(π * x) * cos(π * y)
    V = cos(π * t) * ρ * sin(π * x) * cos(π * y)
    W = cos(π * t) * 0
    E = cos(π * t) * sin(π * x) * cos(π * y) + 100
  end

  P = (γ-1)*(E - (U^2 + V^2 + W^2) / 2ρ)

  u, v, w = U / ρ, V / ρ, W / ρ

  dudx, dudy, dudz = diff(u, x), diff(u, y), diff(u, z)
  dvdx, dvdy, dvdz = diff(v, x), diff(v, y), diff(v, z)
  dwdx, dwdy, dwdz = diff(w, x), diff(w, y), diff(w, z)

  ϵ11 = dudx
  ϵ22 = dvdy
  ϵ33 = dwdz
  ϵ12 = (dudy + dvdx) / 2
  ϵ13 = (dudz + dwdx) / 2
  ϵ23 = (dvdz + dwdy) / 2

  τ11 = 2μ * (ϵ11 - (ϵ11 + ϵ22 + ϵ33) / 3)
  τ22 = 2μ * (ϵ22 - (ϵ11 + ϵ22 + ϵ33) / 3)
  τ33 = 2μ * (ϵ33 - (ϵ11 + ϵ22 + ϵ33) / 3)
  τ12 = τ21 = 2μ * ϵ12
  τ13 = τ31 = 2μ * ϵ13
  τ23 = τ32 = 2μ * ϵ23

  Fx_x = diff.([U;
                u * U + P - τ11;
                u * V     - τ12;
                u * W     - τ13;
                u * (E + P) - u * τ11 - v * τ12 - w * τ13], x)
  Fy_y = diff.([V;
                v * U     - τ21;
                v * V + P - τ22;
                v * W     - τ23;
                v * (E + P) - u * τ21 - v * τ22 - w * τ23], y)
  Fz_z = diff.([W;
                w * U     - τ31;
                w * V     - τ32;
                w * W + P - τ33;
                w * (E + P) - u * τ31 - v * τ32 - w * τ33], z)

  dρdt = simplify(Fx_x[1] + Fy_y[1] + Fz_z[1] + diff(ρ, t))
  dUdt = simplify(Fx_x[2] + Fy_y[2] + Fz_z[2] + diff(U, t))
  dVdt = simplify(Fx_x[3] + Fy_y[3] + Fz_z[3] + diff(V, t))
  dWdt = simplify(Fx_x[4] + Fy_y[4] + Fz_z[4] + diff(W, t))
  dEdt = simplify(Fx_x[5] + Fy_y[5] + Fz_z[5] + diff(E, t))

  @printf output "ρ_g(t, x, y, z, ::Val{%d}) = %s\n" dim ρ
  @printf output "U_g(t, x, y, z, ::Val{%d}) = %s\n" dim U
  @printf output "V_g(t, x, y, z, ::Val{%d}) = %s\n" dim V
  @printf output "W_g(t, x, y, z, ::Val{%d}) = %s\n" dim W
  @printf output "E_g(t, x, y, z, ::Val{%d}) = %s\n" dim E
  @printf output "Sρ_g(t, x, y, z, ::Val{%d}) = %s\n" dim dρdt
  @printf output "SU_g(t, x, y, z, ::Val{%d}) = %s\n" dim dUdt
  @printf output "SV_g(t, x, y, z, ::Val{%d}) = %s\n" dim dVdt
  @printf output "SW_g(t, x, y, z, ::Val{%d}) = %s\n" dim dWdt
  @printf output "SE_g(t, x, y, z, ::Val{%d}) = %s\n" dim dEdt
end

close(output)

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
using SwitchTimeOpt
using Test
using SafeTestsets

@safetestset "Linear system" begin include("./linsys.jl") end
@safetestset "Noninear system" begin include("./nonlinsys.jl") end


text:
mutable struct twoxBC <: BCType
end
function (obj::twoxBC)(params::AdvectionEquationMod.ParamType, u, coords, nrm, t)
  u_bc = 2*coords[1]
  bndryflux = AdvectionEquationMod.RoeSolver(params, u, u_bc, nrm)
  return bndryflux
end

# input file used by all tests
global const test_lowlevel_inputfile = "input_vals_channel.jl"

"""
  Test that the mesh counts things correctly, and that coordinates and metrics
  are calculated correctly.  This makes it dependent on the element ordering
"""
function test_lowlevel_mesh(mesh, sbp, eqn, opts)
  @testset "--- Testing Mesh --- " begin

    @test ( mesh.numVert )== 4
    @test ( mesh.numEdge )== 5
    @test ( mesh.numEl )== 2
    @test ( mesh.order )== 1
    @test ( mesh.numDof )== 4
    @test ( mesh.numNodes )== 4
    @test ( mesh.numDofPerNode )== 1
    @test ( mesh.numBoundaryFaces )== 4
    @test ( mesh.numInterfaces )== 1
    @test ( mesh.numNodesPerElement )== 3
    @test ( mesh.numNodesPerType )== [1, 0 , 0]

    @test ( mesh.bndry_funcs[1] )== AdvectionEquationMod.x5plusy5BC()
    @test ( mesh.bndryfaces[1].element )== 1
    @test ( mesh.bndryfaces[1].face )== 3
    @test ( mesh.bndryfaces[2].element )== 1
    @test ( mesh.bndryfaces[2].face )== 1
    @test ( mesh.bndryfaces[3].element )== 2
    @test ( mesh.bndryfaces[3].face )== 2
    @test ( mesh.bndryfaces[4].element )== 2
    @test ( mesh.bndryfaces[4].face )== 1

    @test ( mesh.interfaces[1].elementL )== 1
    @test ( mesh.interfaces[1].elementR )== 2
    @test ( mesh.interfaces[1].faceL )== 2
    @test ( mesh.interfaces[1].faceR )== 3

    jac_fac = 0.25
    fac = 2
    @test isapprox( mesh.coords[:, :, 2], [4 4 0; 0 4 4.0])
    @test isapprox( mesh.coords[:, :, 1], [0.0 4 0; 0 0 4])

    @test isapprox( mesh.dxidx[:, :, 1, 2], fac*[1 1; -1 0.0]) atol=1e-14
    @test isapprox( mesh.dxidx[:, :, 2, 2], fac*[1 1; -1 0.0]) atol=1e-14
    @test isapprox( mesh.dxidx[:, :, 3, 2], fac*[1 1; -1 0.0]) atol=1e-14

    @test isapprox( mesh.dxidx[:, :, 1, 1], fac*[1 0; 0 1.0]) atol=1e-14
    @test isapprox( mesh.dxidx[:, :, 2, 1], fac*[1 0; 0 1.0]) atol=1e-14
    @test isapprox( mesh.dxidx[:, :, 3, 1], fac*[1 0; 0 1.0]) atol=1e-14

    @test isapprox( mesh.jac, jac_fac*ones(3,2))

  end  # end fact block

  return nothing
end

#test_lowlevel_mesh(mesh, sbp, eqn, opts)
add_func2!(AdvectionTests, test_lowlevel_mesh, test_lowlevel_inputfile, [TAG_SHORTTEST])

"""
  Test some basic functionality: dissassmble solution, mass matrix, etc.
"""
function test_lowlevel_core(mesh, sbp, eqn, opts)
  @testset "--- Testing Functions Within AdvectionData_--- " begin
    Tsol = Float64
    u_vec = Tsol[1,2,3,4]
    u  = zeros(Tsol, 1, 3, 2)

    # checking array1DTo3D
    array1DTo3D(mesh, sbp, eqn, opts, u_vec, u)
    @test isapprox( u[1,1,2], 2.0)
    @test isapprox( u[1,2,2], 4.0)
    @test isapprox( u[1,3,2], 3.0)
    @test isapprox( u[1,1,1], 1.0)
    @test isapprox( u[1,2,1], 2.0)
    @test isapprox( u[1,3,1], 3.0)

    #checking array3DTo1D
    fill!(u_vec, 0.0)
    array3DTo1D(mesh, sbp, eqn, opts, u, u_vec)
    @test isapprox( u_vec, [1.0,4.0,6.0,4.0])

    # check mass matrix
    # just for testing, make jac != 1
    w_val = sbp.w[1]  # all of sbp.w entries are the same
    jac_val = 0.25
    M1 = w_val/jac_val
    M2 = 2*w_val/jac_val
    M3 = 2*w_val/jac_val
    M4 = w_val/jac_val
    M_test = [M1, M2, M3, M4]
    M_code = AdvectionEquationMod.calcMassMatrix(mesh, sbp, eqn)
    @test isapprox( M_code, M_test) atol=1e-13

    Minv_test = 1./M_test
    Minv_test = AdvectionEquationMod.calcMassMatrixInverse(mesh, sbp, eqn)
    @test isapprox( Minv_test, Minv_test) atol=1e-13

    arr = rand(1, 3, 2)
    arr_orig = copy(arr)
    AdvectionEquationMod.matVecA0inv(mesh, sbp, eqn, opts, arr)
    @test isapprox( arr, arr_orig) atol=1e-14

    AdvectionEquationMod.matVecA0(mesh, sbp, eqn, opts, arr)
    @test isapprox( arr, arr_orig) atol=1e-14
  end

  return nothing
end

#test_lowlevel_core(mesh, sbp, eqn, opts)
add_func2!(AdvectionTests, test_lowlevel_core, test_lowlevel_inputfile, [TAG_SHORTTEST])

"""
  Test things from common_funcs.jl
"""
function test_lowlevel_common(mesh, sbp, eqn, opts)
  @testset "--- Testing common functions ---" begin

    x = 1.
    y = 2.
    coords = [x, y]
    eqn.params.alpha_x = 0.0
    eqn.params.alpha_y = 0.0
    t = 0.0
    val = AdvectionEquationMod.calc_x5plusy5(eqn.params, coords, t)
    @test isapprox( val, x^5 + y^5) atol=1e-14

    val = AdvectionEquationMod.calc_exp_xplusy(eqn.params, coords, t)
    @test isapprox( val, exp(x + y)) atol=1e-14
  end  # end facts block

  return nothing
end

#test_lowlevel_common(mesh, sbp, eqn, opts)
add_func2!(AdvectionTests, test_lowlevel_common, test_lowlevel_inputfile, [TAG_SHORTTEST])

"""
  Test Roe solver used for boundary conditions and some hand calculated
  boundary flux values.
"""
function test_lowlevel_bc(mesh, sbp, eqn, opts)
  @testset "--- Testing Boundary Function ---" begin

    # testing choice of u or u_bc
    u = 5.0
    u_bc = 2.5
    eqn.params.alpha_x = 1.5
    eqn.params.alpha_y = 0.5
    dxidx = [1. 0; 0 1]
    nrm = [1., 0]
    nrm2 = zeros(nrm)
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)

    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)

    @test isapprox( val, u*eqn.params.alpha_x) atol=1e-14

    nrm = [-1.0, 0]
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)
    @test isapprox( val, -u_bc*eqn.params.alpha_x) atol=1e-14

    nrm = [0, 1.0]
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)
    @test isapprox( val, u*eqn.params.alpha_y) atol=1e-14

    nrm = [0, -1.0]
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)
    @test isapprox( val, -u_bc*eqn.params.alpha_y) atol=1e-14

    # now test rotation using dxidx

    function get_rotation_matrix(theta)
      return [cos(theta) -sin(theta); sin(theta) cos(theta)]
    end

    theta = 30*pi/180  # angle between x axis and xi axis
    flow_direction = 25*pi/180
    alpha_mag = 2.0
    eqn.params.alpha_x = alpha_mag*cos(flow_direction)
    eqn.params.alpha_y = alpha_mag*sin(flow_direction)

    dxidx = get_rotation_matrix( theta)

    # nrm is in the xi-eta space
    # the xi axis is 30 degrees clockwise of the x axis, so the alpha vector
    # which is already at 25 degrees from the x axis needs to go 30 more
    # to be in the xi-eta coordinate system
    # this is still an outflow for nrm = [1, 0]
    nrm = [1., 0]  # flow is in the xi direction
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    angle_diff = theta + flow_direction
    alpha_eff = alpha_mag*cos(angle_diff)  # effective alpha in the wall normal
                                           # direction
    val_exp = alpha_eff*u
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)

    @test isapprox( val, val_exp) atol=1e-14

    # now check eta direction
    nrm = [0, 1.]
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    alpha_eff = alpha_mag*sin(angle_diff)
    val_exp = u*alpha_eff
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)

    @test isapprox( val, val_exp) atol=1e-14

    # now rotate the coordinate system so much that this becomes an inflow
    theta = 120*pi/180
    nrm = [1., 0]
    angle_diff = theta + flow_direction
    alpha_eff = alpha_mag*cos(angle_diff)
    val_exp = alpha_eff*u_bc

    dxidx = get_rotation_matrix( theta)
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)
    @test isapprox( val, val_exp) atol=1e-14

    # check eta direction
    nrm = [0, 1]
    calcBCNormal(eqn.params, dxidx, nrm, nrm2)
    angle_diff = theta + flow_direction
    alpha_eff = alpha_mag*sin(angle_diff)
    val_exp = alpha_eff*u
    val = AdvectionEquationMod.RoeSolver(eqn.params, u, u_bc, nrm2)
    @test isapprox( val, val_exp) atol=1e-14
  end  # end facts block

  return nothing
end

#test_lowlevel_bc(mesh, sbp, eqn, opts)
add_func2!(AdvectionTests, test_lowlevel_bc, test_lowlevel_inputfile, [TAG_BC, TAG_SHORTTEST])

"""
  Test computing volume integrals over entire mesh
"""
function test_lowlevel_volumeintegrals()
  @testset "--- Testing Volume Integrals ---" begin

    # use the 8 element mesh
    fname = "input_vals_8el.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    println("precompute_boundary_flux = ", opts["precompute_boundary_flux"])
    Tmsh = eltype(mesh.dxidx)

    fill!(eqn.q, 0.0)
    eqn.params.alpha_x = 1.0
    eqn.params.alpha_y = 1.0
#    eqn.params.alpha_x = ones(Float64 ,1, mesh.numNodesPerElement, mesh.numEl)
    alpha_y = zero(eqn.params.alpha_x)

    fill!(eqn.res, 0.0)
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    @test isapprox( eqn.res, zeros(1, mesh.numNodesPerElement, mesh.numEl)) atol=1e-12

    # check that Qx.'*q = 0 when q = ones(3)
    fill!(eqn.q, 1.0)
    fill!(eqn.res, 0.0)
    dxidx1 = mesh.dxidx[:, :, :, 1]
    jac = mesh.jac[:, 1]
    dxidx_true = zeros(2, 2, 3)
    for i=1:3
      dxidx_true[:, :, i] = dxidx1[:, :, i]*jac[i]
    end

    Qx = sbp.Q[:, :, 1]*dxidx1[1, 1, 1] + sbp.Q[:, :, 2]*dxidx1[2, 1, 1]
    q = [eqn.q[1, 1, 1], eqn.q[1, 2, 1], eqn.q[1, 3, 1]]  # extract q values
    # sum reduces the vector to the value of the integral
    val_test = sum(Qx.'*q)
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    # extract residual values
    vec_code =  [eqn.res[1, 1, 1], eqn.res[1, 2, 1], eqn.res[1, 3, 1]]
    val_code = sum(vec_code)
    xi_flux = dxidx1[1,1,1]*1*q
    xi_component = sbp.Q[:, :, 1].'*xi_flux
    eta_flux = dxidx1[2, 1]*1*q
    eta_component = sbp.Q[:, :, 2].'*eta_flux
    val2_test = sum(xi_component + eta_component)
    @test isapprox( val_test, val2_test) atol=1e-14
    @test isapprox( val_test, val_code) atol=1e-14

    # test that the integral of qx dOmega^e when q = 2*x works
    # dq/dx = 2, thus the integral is = 1

    println("----- Checking q=2*x case -----")
    bc_func = twoxBC()
    mesh.bndry_funcs[1] = bc_func

    x1 = zeros(Tmsh, 3)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        eqn.q[1, j, i] = 2*x
        if i == 1
          x1[j] = x
        end

      end
    end

    fill!(eqn.res, 0.0)
    eqn.params.alpha_x = 1.0
    eqn.params.alpha_y = 0.0
    # check the boundry contribution
    AdvectionEquationMod.evalBoundaryIntegrals(mesh, sbp, eqn, opts)

    @test isapprox( sum(eqn.res[:, :, 1]), -2.0) atol=1e-14
    @test isapprox( sum(eqn.res[:, :, 3]), -2.0) atol=1e-14

    @test isapprox( sum(eqn.res[:, :, 6]), -2.0) atol=1e-14
    @test isapprox( sum(eqn.res[:, :, 8]), -2.0) atol=1e-14

    fill!(eqn.res, 0.0)
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    for i=1:mesh.numEl
      Qx_i = sbp.Q[:, :, 1]*mesh.dxidx[1, 1, 1, i] + sbp.Q[:, :, 2]*mesh.dxidx[2, 1, 1, i]
      q_i = reshape(eqn.q[1, :, i], 3)
      val_test = Qx_i.'*q_i
      val_code = reshape(eqn.res[:, :, i], 3)
      @test isapprox( sum(val_code), 0.0) atol=1e-14# proven by hand calc
      @test isapprox( val_code, val_test) atol=1e-14
    end

    println("----- Checking q=2*x^2 + 5 case -----")

    x1 = zeros(Tmsh, 3)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        eqn.q[1, j, i] = 2*x*x + 5
      end
    end

    fill!(eqn.res, 0.0)
    eqn.params.alpha_x = 1.0
    eqn.params.alpha_y = 0.0
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    for i=1:mesh.numEl
      Qx_i = sbp.Q[:, :, 1]*mesh.dxidx[1, 1, 1, i] + sbp.Q[:, :, 2]*mesh.dxidx[2, 1, 1, i]
      q_i = reshape(eqn.q[1, :, i], 3)
      val_test = Qx_i.'*q_i
      val_code = reshape(eqn.res[:, :, i], 3)
      @test isapprox( val_code, val_test) atol=1e-14
    end

    fname = "input_vals_8el_large.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    println("----- Checking q=2*x^2 + 5 case on large grid -----")

    x1 = zeros(Tmsh, 3)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        eqn.q[1, j, i] = 2*x*x + 5
      end
    end

    fill!(eqn.res, 0.0)
    eqn.params.alpha_x = 1.0
    eqn.params.alpha_y = 0.0
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    for i=1:mesh.numEl
      Qx_i = sbp.Q[:, :, 1]*mesh.dxidx[1, 1, 1, i] + sbp.Q[:, :, 2]*mesh.dxidx[2, 1, 1, i]
      q_i = reshape(eqn.q[1, :, i], 3)
      val_test = Qx_i.'*q_i
      val_code = reshape(eqn.res[:, :, i], 3)
      @test isapprox( val_code, val_test) atol=1e-14
    end

    # back to the original mesh
    fname = "input_vals_8el.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    println("----- Checking sinwave case -----")
    x1 = zeros(Tmsh, 3)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        alpha_x = eqn.params.alpha_x
        alpha_y = eqn.params.alpha_y
        eqn.q[1, j, i] = AdvectionEquationMod.calc_sinwave(eqn.params, mesh.coords[:, j, i], 0.25)
      end
    end

    mesh.bndry_funcs[1] = AdvectionEquationMod.sinwave_BC()
    fill!(eqn.res, 0.0)
    eqn.params.alpha_x = 1.0
    eqn.params.alpha_y = 0.0
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    for i=1:mesh.numEl
      Qx_i = sbp.Q[:, :, 1]*mesh.dxidx[1, 1, 1, i] + sbp.Q[:, :, 2]*mesh.dxidx[2, 1, 1, i]
      q_i = reshape(eqn.q[1, :, i], 3)
      val_test = Qx_i.'*q_i
      val_code = reshape(eqn.res[:, :, i], 3)
      @test isapprox( val_code, val_test) atol=1e-14
    end

    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
#TODO: uncomment when SBP boundaryintegrate is fixed
#=
    fname = "input_vals_channel_verylarge.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    println("----- Checking sinwave case on very large mesh -----")
    x1 = zeros(Tmsh, 3)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        alpha_x = eqn.params.alpha_x
        alpha_y, = eqn.params.alpha_y

        eqn.q[1, j, i] = AdvectionEquationMod.calc_sinwave(mesh.coords[:, j, i], eqn.params, 0.25)
      end
    end
    println("finished calculating sinwave")
    mesh.bndry_funcs[1] = AdvectionEquationMod.sinwave_BC()
    fill!(eqn.res, 0.0)
    eqn.params.alpha_x = 1.0
    eqn.params.alpha_y = 0.0
    AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)
    println("called eval SCResidual")
    for i=1:mesh.numEl
      Qx_i = sbp.Q[:, :, 1]*mesh.dxidx[1, 1, 1, i] + sbp.Q[:, :, 2]*mesh.dxidx[2, 1, 1, i]
      q_i = reshape(eqn.q[1, :, i], 3)
      val_test = Qx_i.'*q_i
      val_code = reshape(eqn.res[:, :, i], 3)
      @test isapprox( val_code, val_test) atol=1e-14
    end
  =#
  end # end facts block

  return nothing
end

#test_lowlevel_volumeintegrals()
add_func1!(AdvectionTests, test_lowlevel_volumeintegrals, [TAG_VOLUMEINTEGRALS, TAG_SHORTTEST])
  #=
  @testset "--- Testing evalInteriorFlux ---" begin

    fill!(eqn.res, 0.0)
    fill!(eqn.params.alpha_x, 1.0)
    fill!(eqn.params.alpha_y, 1.0)

    nbrnodeindex = Array(mesh.numNodesPerFace:-1:1)

  end
  =#

#=
 @testset "--- Testing dataPrep ---" begin

  end

  @testset "--- Testing evalBoundaryIntegrals ---" begin
  end

  @testset "--- Testing evalEuler --- " begin

end # end facts block
=#


text:
# skip

using Bukdu

function mock_action(endpoint, path, headers=[], query=nothing)
    if query !== nothing
        path = string(merge(HTTP.URI(path), query=query))
    end
    Router.call(get, path, headers).resp
end

merge!(Plug.Loggers.config, Dict(:action_pad => 25, :path_pad => 50))


text:
#= Sysimage for Developers. Sysimage.so must be written on .gitignore because it is so large file! =#
using Pkg
Pkg.add("PackageCompiler")
Pkg.add("DelimitedFiles")
Pkg.add("TimerOutputs")
Pkg.add("Test")
Pkg.add("Plots")
using PackageCompiler

PackageCompiler.create_sysimage([:CSV, :DataFrames, :Plots, :Luxor, :BenchmarkTools, :TimerOutputs, :Test]; sysimage_path="Sysimage.so")
# >>  You can use packages with 'import and using'
# $ julia -JSysimage.so

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
# findpeaks.jl
# Compare Matlab and Sound versions of findpeaks.

using Sound: findpeaks
using MATLAB
using Test: @testset, @test

ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"

@testset "findpeaks" begin
    N = 32
    i = [7, 12, 25]
    M = length(i)
    x = zeros(N)
    x[i] .= [2, 5, 3]
    pj = findpeaks(x)

    mat"pm = findPeaks4($x, 50, 0.005, 0)"
    @mget pm
    pm == Int.(pm)
    @test pj == pm
end