text:
const Height = Int8
const TreeNode{T, I, D} = Tuple{T, I, I, Height, D}

mutable struct Node{T, D}
    value::T
    left::Union{Nothing, Node{T, D}}
    right::Union{Nothing, Node{T, D}}
    height::Height
    data::D
end
const Tree{T, D} = Union{Nothing, Node{T, D}}

@inline nil() = nothing
@inline node(t::T, data::D) where {T, D} =
    Node{T, D}(t, nothing, nothing, Height(1), data)
@inline node(t::T, left, right, height, data::D) where {T, D} =
    Node{T, D}(t, left, right, height, data)

@inline nil(alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}} = 0
@inline node(t::T, data::D, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}} =
    allocateend!(alloc, (t, 0, 0, Height(1), data))
@inline node(t::T, left, right, height, data::D, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}} =
    allocateend!(alloc, (t, left, right, height, data))

@inline deallocate(i::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}} =
    deallocate!(alloc, i)

# @inline value(tree::Nothing) = @assert false
@inline value(tree::Node{T, D}) where {T, D} = tree.value

# @inline left(tree::Nothing) = @assert false
@inline left(tree::Node{T, D}) where {T, D} = tree.left

# @inline right(tree::Nothing) = @assert false
@inline right(tree::Node{T, D}) where {T, D} = tree.right

@inline height(tree::Nothing) = Height(0)
@inline height(tree::Node{T, D}) where {T, D} = tree.height

# @inline data(tree::Nothing) = @assert false
@inline data(tree::Node{T, D}) where {T, D} = tree.data

# @inline value!(tree::Nothing, left::Tree{T}) where {T, D} = @assert false
@inline function value!(tree::Node{T, D}, value::T) where {T, D}
    tree.value = value
end

# @inline left!(tree::Nothing, left::Tree{T}) where {T, D} = @assert false
@inline function left!(tree::Node{T, D}, left::Tree{T}) where {T, D}
    tree.left = left
end

# @inline right!(tree::Nothing, right::Tree{T}) where {T, D} = @assert false
@inline function right!(tree::Node{T, D}, right::Tree{T}) where {T, D}
    tree.right = right
end

# @inline height!(tree::Nothing, height::Height) = @assert false
@inline function height!(tree::Node{T, D}, height::Height) where {T, D}
    tree.height = height
end

# @inline height!(tree::Nothing, height_left::Height, height_right::Height) = @assert false
@inline function height!(tree::Node{T, D}, height_left::Height, height_right::Height) where {T, D}
    tree.height = max(height_left, height_right) + Height(1)
end

# @inline height!(tree::Nothing) = @assert false
@inline function height!(tree::Node{T, D}) where {T, D}
    _left = left(tree)
    _right = right(tree)
    height!(tree, height(_left), height(_right))
end

# @inline data!(tree::Nothing, data::D) where {T, D} = @assert false
@inline function data!(tree::Node{T, D}, data::D) where {T, D}
    tree.data = data
end

@inline function value(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    alloc[tree][1]
end
@inline function value(tree::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 1]
end

@inline function left(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    alloc[tree][2]
end
@inline function left(tree::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 2]
end

@inline function right(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    alloc[tree][3]
end
@inline function right(tree::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 3]
end

@inline function height(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    tree === 0 ? Height(0) : alloc[tree][4]
end
@inline function height(tree::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    tree === 0 ? Height(0) : alloc[tree, 4]
end

@inline function data(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    alloc[tree][5]
end
@inline function data(tree::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 5]
end

@inline function value!(tree::I, _value::T, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    # if tree !== 0
        node = alloc[tree]
        # if _value != node[1]
            t = (_value, node[2], node[3],
                 node[4], node[5])
            alloc[tree] = t
        # end
    # else
    #     @assert false
    # end
end
@inline function value!(tree::I, _value::T, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 1] = _value
end

@inline function left!(tree::I, _left::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    # if tree !== 0
        node = alloc[tree]
        # if _left != node[2]
            t = (node[1], _left, node[3],
                 node[4], node[5])
            alloc[tree] = t
        # end
    # else
    #     @assert false
    # end
end
@inline function left!(tree::I, _left::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 2] = _left
end

@inline function right!(tree::I, _right::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    # if tree !== 0
        node = alloc[tree]
        # if _right != node[3]
            t = (node[1], node[2], _right,
                 node[4], node[5])
            alloc[tree] = t
        # end
    # else
    #     @assert false
    # end
end
@inline function right!(tree::I, _right::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 3] = _right
end

@inline function height!(tree::I, _height::Height, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    # if tree !== 0
        node = alloc[tree]
        # if _height != node[5]
            t = (node[1], node[2], node[3],
                 _height, node[5])
            alloc[tree] = t
        # end
    # else
    #     @assert false
    # end
end
@inline function height!(tree::I, _height::I, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 4] = _height
end

@inline function height!(tree::I, height_left::Height, height_right::Height, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    _height = max(height_left, height_right) + Height(1)
    height!(tree, _height, alloc)
end

@inline function height!(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    # if tree !== 0
        _left = left(tree, alloc)
        _right = right(tree, alloc)
        height_left = height(_left, alloc)
        height_right = height(_right, alloc)
        height!(tree, height_left, height_right, alloc)
    # else
    #     @assert false
    # end
end

@inline function data!(tree::I, _data::D, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    # if tree !== 0
        node = alloc[tree]
        # if _data != node[5]
            t = (node[1], node[2], node[3],
                 node[4], _data)
            alloc[tree] = t
        # end
    # else
    #     @assert false
    # end
end
@inline function data!(tree::I, _data::D, alloc::A) where {T, I, D, S, A <: AbstractFreeListSOAllocator{TreeNode{T, I, D}, I, S}}
    alloc[tree, 5] = _data
end

invariantvalue(tree::Nothing) = true
function invariantvalue(tree::Node{T, D}) where {T, D}
    invariant = true
    _value = value(tree)
    _left = left(tree)
    invariant &= _left === nothing || value(_left) < _value
    _right = right(tree)
    invariant &= _right === nothing || value(_right) > _value
    @assert invariant
    invariant
end

invariantheight(tree::Nothing) = true
function invariantheight(tree::Node{T, D}) where {T, D}
    _left = left(tree)
    _right = right(tree)
    _height = height(tree)
    height_left = height(_left)
    height_right = height(_right)
    @assert _height == max(height_left, height_right) + 1
    invariant = abs(height_left - height_right) < 2
    @assert invariant
    invariant
end

function invariantvalue(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    if tree === 0
        return true
    else
        invariant = true
        _value = value(tree, alloc)
        _left = left(tree, alloc)
        invariant &= _left === 0 || value(_left, alloc) < _value
        _right = right(tree, alloc)
        invariant &= _right === 0 || value(_right, alloc) > _value
        @assert invariant
        invariant
    end
end

function invariantheight(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    if tree === 0
        return true
    else
        _left = left(tree, alloc)
        _right = right(tree, alloc)
        _height = height(tree, alloc)
        height_left = height(_left, alloc)
        height_right = height(_right, alloc)
        @assert _height == max(height_left, height_right) + 1
        invariant = abs(height_left - height_right) < 2
        @assert invariant
        invariant
    end
end

@inline function rotateright(tree::Node{T, D}) where {T, D}
    _left = left(tree)
    left_right = right(_left)
    right!(_left, tree)
    left!(tree, left_right)
    # height!(tree)
    _height = height!(tree, height(left_right), height(right(tree)))
    # height!(_left)
    height!(_left, height(left(_left)), _height)
    _left
end
@inline function rotateleft(tree::Node{T, D}) where {T, D}
    _right = right(tree)
    right_left = left(_right)
    @assert _right !== nothing
    left!(_right, tree)
    right!(tree, right_left)
    # height!(tree)
    _height = height!(tree, height(left(tree)), height(right_left))
    # height!(_right)
    height!(_right, _height, height(right(_right)))
    _right
end

@inline function rotateright(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    _left = left(tree, alloc)
    left_right = right(_left, alloc)
    right!(_left, tree, alloc)
    left!(tree, left_right, alloc)
    # height!(tree, alloc)
    _right = right(tree, alloc)
    height!(tree, height(left_right, alloc), height(_right, alloc), alloc)
    # height!(_left, alloc)
    left_left = left(_left, alloc)
    height!(_left, height(left_left, alloc), height(tree, alloc), alloc)
    _left
end
@inline function rotateleft(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    _right = right(tree, alloc)
    right_left = left(_right, alloc)
    left!(_right, tree, alloc)
    right!(tree, right_left, alloc)
    # height!(tree, alloc)
    _left = left(tree, alloc)
    height!(tree, height(_left, alloc), height(right_left, alloc), alloc)
    # height!(_right, alloc)
    right_right = right(_right, alloc)
    height!(_right, height(tree, alloc), height(right_right, alloc), alloc)
    _right
end

# @inline balance(tree::Nothing) = nothing
@inline function balance(tree::Node{T, D}) where {T, D}
    _left = left(tree)
    _right = right(tree)
    height_left = height(_left)
    height_right = height(_right)
    if height_left < height_right - 1
        right_left = left(_right)
        right_right = right(_right)
        height_right_left = height(right_left)
        height_right_right = height(right_right)
        if height_right_left > height_right_right
            _right = rotateright(_right)
            right!(tree, _right)
        end
        tree = rotateleft(tree)
    elseif height_left > height_right + 1
        left_left = left(_left)
        left_right = right(_left)
        height_left_left = height(left_left)
        height_left_right = height(left_right)
        if height_left_left < height_left_right
            _left = rotateleft(_left)
            left!(tree, _left)
        end
        tree = rotateright(tree)
    else
        height!(tree, height_left, height_right)
    end
    tree
end

@inline function balance(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    _left = left(tree, alloc)
    _right = right(tree, alloc)
    height_left = height(_left, alloc)
    height_right = height(_right, alloc)
    if height_left < height_right - 1
        right_left = left(_right, alloc)
        right_right = right(_right, alloc)
        height_right_left = height(right_left, alloc)
        height_right_right = height(right_right, alloc)
        if height_right_left > height_right_right
            _right = rotateright(_right, alloc)
            right!(tree, _right, alloc)
        end
        tree = rotateleft(tree, alloc)
    elseif height_left > height_right + 1
        left_left = left(_left, alloc)
        left_right = right(_left, alloc)
        height_left_left = height(left_left, alloc)
        height_left_right = height(left_right, alloc)
        if height_left_left < height_left_right
            _left = rotateleft(_left, alloc)
            left!(tree, _left, alloc)
        end
        tree = rotateright(tree, alloc)
    else
        height!(tree, height_left, height_right, alloc)
    end
    tree
end

@noinline insert(tree::Nothing, t::T, data::D) where {T, D} = node(t, data)
@noinline function insert(tree::Node{T}, t::T, data::D) where {T, D}
    _value = value(tree)
    if t < _value
        _left = left(tree)
        _height = height(_left)
        _left = insert(_left, t, data)
        left!(tree, _left)
        if height(_left) != _height
            tree = balance(tree)
        end
    elseif t > _value
        _right = right(tree)
        _height = height(_right)
        _right = insert(_right, t, data)
        right!(tree, _right)
        if height(_right) != _height
            tree = balance(tree)
        end
    else
        data!(tree, data)
    end
    tree
end

@noinline function insert(tree::I, t::T, data::D, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    if tree !== 0
        _value = value(tree, alloc)
        if t < _value
            _left = left(tree, alloc)
            _height = height(_left, alloc)
            _left = insert(_left, t, data, alloc)
            left!(tree, _left, alloc)
            if height(_left, alloc) != _height
                tree = balance(tree, alloc)
            end
        elseif t > _value
            _right = right(tree, alloc)
            _height = height(_right, alloc)
            _right = insert(_right, t, data, alloc)
            right!(tree, _right, alloc)
            if height(_right, alloc) != _height
                tree = balance(tree, alloc)
            end
        else
            data!(tree, data, alloc)
        end
    else
        tree = node(t, data, alloc)
    end
    tree
end

# @inline Base.haskey(tree::Nothing, t::T) where {T, D} = false
@inline function Base.haskey(tree::Tree{T, D}, t::T) where {T, D}
    while tree !== nothing
        _value = value(tree)
        if t < _value
            tree = left(tree)
        elseif t > _value
            tree = right(tree)
        else
            return true
        end
    end
    return false
end

@inline function Base.haskey(tree::I, t::T, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    while tree !== 0
        _value = value(tree, alloc)
        if t < _value
            tree = left(tree, alloc)
        elseif t > _value
            tree = right(tree, alloc)
        else
            return true
        end
    end
    return false
end

@inline Base.getindex(tree::Nothing, t::T) where {T, D} = @assert false
@inline function Base.getindex(tree::Node{T, D}, t::T) where {T, D}
    while tree !== nothing
        _value = value(tree)
        if t < _value
            tree = left(tree)
        elseif t > _value
            tree = right(tree)
        else
            return data(tree)
        end
    end
    @assert false
end

@inline function Base.haskey(tree::I, t::T, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    while tree !== 0
        _value = value(tree, alloc)
        if t < _value
            tree = left(tree, alloc)
        elseif t > _value
            tree = right(tree, alloc)
        else
            return true
        end
    end
    return false
end

@inline function Base.getindex(tree::I, t::T, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    while tree !== 0
        _value = value(tree, alloc)
        if t < _value
            tree = left(tree, alloc)
        elseif t > _value
            tree = right(tree, alloc)
        else
            return data(tree, alloc)
        end
    end
    @assert false
end

@inline getleftmost(::Nothing) = @assert false
@inline function getleftmost(tree::Node{T, D}) where {T, D}
    _left = left(tree)
    while _left !== nothing
        tree = _left
        _left = left(tree)
    end
    value(tree)
end

@inline getrightmost(::Nothing) = @assert false
@inline function getrightmost(tree::Node{T, D}) where {T, D}
    _right = right(tree)
    while _right !== nothing
        tree = _right
        _right = right(tree)
    end
    value(tree)
end

@inline function getleftmost(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    _left = left(tree, alloc)
    while _left !== 0
        tree = _left
        _left = left(tree, alloc)
    end
    value(tree, alloc)
end

@inline function getrightmost(tree::I, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    _right = right(tree, alloc)
    while _right !== 0
        tree = _right
        _right = right(tree, alloc)
    end
    value(tree, alloc)
end

@noinline delete(tree::Nothing, t::T) where T = tree
@noinline function delete(tree::Node{T, D}, t::T) where {T, D}
    _value = value(tree)
    if t < _value
        _left = left(tree)
        _height = height(_left)
        _left = delete(_left, t)
        left!(tree, _left)
        if height(_left) != _height
            tree = balance(tree)
        end
    elseif t > _value
        _right = right(tree)
        _height = height(_right)
        _right = delete(_right, t)
        right!(tree, _right)
        if height(_right) != _height
            tree = balance(tree)
        end
    else
        _left = left(tree)
        _right = right(tree)
        if _left === nothing && _right === nothing
            tree = nothing
        elseif _left === nothing
            value!(tree, value(_right))
            left!(tree, left(_right))
            right!(tree, right(_right))
            # height!(tree)
            height!(tree, height(_right))
            data!(tree, data(_right))
        elseif _right === nothing
            value!(tree, value(_left))
            left!(tree, left(_left))
            right!(tree, right(_left))
            # height!(tree)
            height!(tree, height(_left))
            data!(tree, data(_left))
        else
            height_left = height(_left)
            height_right = height(_right)
            if height_left >= height_right
                _value = getrightmost(_left)
                value!(tree, _value)
                _left = delete(_left, _value)
                left!(tree, _left)
                # height!(tree)
                height!(tree, height(_left), height_right)
            else
                _value = getleftmost(_right)
                value!(tree, _value)
                _right = delete(_right, _value)
                right!(tree, _right)
                # height!(tree)
                height!(tree, height_left, height(_right))
            end
        end
    end
    tree
end

@noinline function delete(tree::I, t::T, alloc::A) where {T, I, D, A <: AbstractFreeListAllocator{TreeNode{T, I, D}, I}}
    if tree !== 0
        _value = value(tree, alloc)
        if t < _value
            _left = left(tree, alloc)
            _height = height(_left, alloc)
            _left = delete(_left, t, alloc)
            left!(tree, _left, alloc)
            if height(_left, alloc) != _height
                tree = balance(tree, alloc)
            end
        elseif t > _value
            _right = right(tree, alloc)
            _height = height(_right, alloc)
            _right = delete(_right, t, alloc)
            right!(tree, _right, alloc)
            if height(_right, alloc) != _height
                tree = balance(tree, alloc)
            end
        else
            _left = left(tree, alloc)
            _right = right(tree, alloc)
            if _left === 0 && _right === 0
                deallocate(tree, alloc)
                tree = 0
            elseif _left === 0
                alloc[tree] = alloc[_right]
                deallocate(_right, alloc)
            elseif _right === 0
                alloc[tree] = alloc[_left]
                deallocate(_left, alloc)
            else
                height_left = height(_left, alloc)
                height_right = height(_right, alloc)
                if height_left >= height_right
                    _value = getrightmost(_left, alloc)
                    value!(tree, _value, alloc)
                    _left = delete(_left, _value, alloc)
                    left!(tree, _left, alloc)
                    # height!(tree, alloc)
           		    height!(tree, height(_left, alloc), height_right, alloc)
                else
                    _value = getleftmost(_right, alloc)
                    value!(tree, _value, alloc)
                    _right = delete(_right, _value, alloc)
                    right!(tree, _right, alloc)
                    # height!(tree, alloc)
	                height!(tree, height_left, height(_right, alloc), alloc)
                end
            end
        end
    end
    tree
end


text:
# ---
# title: 478. Generate Random Point in a Circle
# id: problem478
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Random, Rejection Sampling
# link: <https://leetcode.com/problems/generate-random-point-in-a-circle/description/>
# hidden: true
# ---
#
# Given the radius and x-y positions of the center of a circle, write a function
# `randPoint` which generates a uniform random point in the circle.
#
# Note:
#
#   1. input and output values are in [floating-point](https://www.webopedia.com/TERM/F/floating_point_number.html).
#   2. radius and x-y position of the center of the circle is passed into the class constructor.
#   3. a point on the circumference of the circle is considered to be in the circle.
#   4. `randPoint` returns a size 2 array containing x-position and y-position of the random point, in that order.
#
# **Example 1:**
#
#
#
#     Input: ["Solution","randPoint","randPoint","randPoint"]
#     [[1,0,0],[],[],[]]
#     Output: [null,[-0.72939,-0.65505],[-0.78502,-0.28626],[-0.83119,-0.19803]]
#
#
# **Example 2:**
#
#
#
#     Input: ["Solution","randPoint","randPoint","randPoint"]
#     [[10,5,-7.5],[],[],[]]
#     Output: [null,[11.52438,-8.33273],[2.46992,-16.21705],[11.13430,-12.42337]]
#
# **Explanation of Input Syntax:**
#
# The input is two lists: the subroutines called and their arguments.
# `Solution`'s constructor has three arguments, the radius, x-position of the
# center, and y-position of the center of the circle. `randPoint` has no
# arguments. Arguments are always wrapped with a list, even if there aren't any.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
using ADCME
using PyCall
using LinearAlgebra
using SparseArrays
using PyPlot
using Random
using Test
Random.seed!(233)

sparse_mat_mul = load_op("build/libSparseMatMul", "sparse_sparse_mat_mul")
diag_sparse_mat_mul = load_op("build/libSparseMatMul", "diag_sparse_mat_mul")
sparse_diag_mat_mul = load_op("build/libSparseMatMul", "sparse_diag_mat_mul")

function Base.:*(s1::SparseTensor, s2::SparseTensor)
    ii1, jj1, vv1 = find(s1)
    ii2, jj2, vv2 = find(s2)
    m, n = size(s1)
    n_, k = size(s2)
    if n!=n_
        error("IGACS: matrix size mismatch: ($m, $n) vs ($n_, $k)")
    end
    mat_mul_fn = sparse_mat_mul
    if s1._diag
        mat_mul_fn = diag_sparse_mat_mul
    elseif s2._diag
        mat_mul_fn = sparse_diag_mat_mul
    end
    ii3, jj3, vv3 = mat_mul_fn(ii1-1,jj1-1,vv1,ii2-1,jj2-1,vv2,m,n,k)
    SparseTensor(ii3, jj3, vv3, m, k)
end
# sparse_mat_mul = py"sparse_mat_mul"
################## End Load Operator ##################
A = sprand(10,5,0.3)
B = sprand(5,20,0.3)
C = A*B
CC = SparseTensor(A)*SparseTensor(B)
# TODO: specify your input parameters
sess = Session()
init(sess)
C_ = run(sess, CC)
@test C_≈C

A = spdiagm(0=>[1.;2.;3;4;5])
B = sprand(5,20,0.3)
C = A*B
CC = SparseTensor(A)*SparseTensor(B)
# TODO: specify your input parameters
sess = Session()
init(sess)
C_ = run(sess, CC)
@test C_≈C

A = sprand(10,5,0.5)
B = spdiagm(0=>[1.;2.;3;4;5])
C = A*B
CC = SparseTensor(A)*SparseTensor(B)
# TODO: specify your input parameters
sess = Session()
init(sess)
C_ = run(sess, CC)
@test C_≈C


text:
using VDPTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter

frames = Frames(MIME("image/png"), fps=2)

dmdp = DiscreteVDPTagMDP()
mdp = cproblem(dmdp)

hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(1))
policy = translate_policy(ToNextML(mdp), mdp, dmdp, dmdp)
hist = simulate(hr, dmdp, policy)

cstates = [convert_s(TagState, s, dmdp) for s in state_hist(hist)]
gr()
@showprogress "Creating gif..." for s in cstates
    push!(frames, plot(mdp, s))
end

filename = string(tempname(), "_vdprun.gif")
write(filename, frames)
println(filename)
run(`setsid gifview $filename`)


text:
# Part 1 of final project: Diffusion equation

@views function init()
    nx   = 32             # number of global grid points
    ny   = 32             # number of global grid points
    nz   = 32             # number of global grid points
    xc   = zeros(nx, 1) # global coord vector
    H    = zeros(nx, ny, nz) # global solution as obtained by implicitGlobalGrid's `gather!()`
    inds = Int.(ceil.(LinRange(1, length(xc), 12)))
    xc[inds]          .= [0.15625, 1.09375, 2.03125, 2.96875, 3.90625, 4.843750000000001, 5.468749999999999, 6.40625, 7.34375, 8.28125, 9.21875, 9.84375]
    H[inds, inds, 15] .= [6.711882385110737e-21 2.4514117276084958e-17 1.5437604567532982e-14 1.6762365118986786e-12 3.138212957188894e-11 1.0130275516790658e-10 8.332937341028532e-11 1.4367785881450301e-11 4.2714301545963265e-13 2.1895179117406884e-15 1.9351535299019154e-18 6.711882385110737e-21; 2.4514117276084958e-17 0.0004814874277834956 0.0014445717212427402 0.0032493045299438012 0.005557594231157323 0.006946731379668997 0.006690697700021068 0.004803038520579223 0.0025518750281671872 0.0010411002857436305 0.00029141318932587976 2.4514117276084958e-17; 1.5437604567532982e-14 0.0014445717212427402 0.004624363091812729 0.011159479624478753 0.020196376510907457 0.025913346399889084 0.02484567081606447 0.017173580755691007 0.008566902048912462 0.003255743161934395 0.0008608124023284637 1.5437604567532982e-14; 1.6762365118986786e-12 0.0032493045299438012 0.011159479624478753 0.029095533074715414 0.056125785858625855 0.07421776052483857 0.0707876802440615 0.04685098387699614 0.02177896579933202 0.007659758643077384 0.0019037581454307546 1.6762365118986786e-12; 3.138212957188894e-11 0.005557594231157323 0.020196376510907457 0.056125785858625855 0.11433141121531959 0.15522789812463358 0.1473718066389615 0.0939249428097469 0.04113536620718002 0.013581240422684215 0.0032124500970325837 3.138212957188894e-11; 1.0130275516790658e-10 0.006946731379668997 0.025913346399889084 0.07421776052483857 0.15522789812463358 0.21350368660148342 0.20223719296375614 0.12653170317321735 0.053851078219800314 0.017260273812497083 0.003990509147237947 1.0130275516790658e-10; 8.332937341028532e-11 0.006690697700021068 0.02484567081606447 0.0707876802440615 0.1473718066389615 0.20223719296375614 0.19164191533936806 0.12029222684071515 0.051452364828806635 0.01657647018690786 0.0038475611129833603 8.332937341028532e-11; 1.4367785881450301e-11 0.004803038520579223 0.017173580755691007 0.04685098387699614 0.0939249428097469 0.12653170317321735 0.12029222684071515 0.07752458731237559 0.03454914825401289 0.011616853176862237 0.0027870615866269714 1.4367785881450301e-11; 4.2714301545963265e-13 0.0025518750281671872 0.008566902048912462 0.02177896579933202 0.04113536620718002 0.053851078219800314 0.051452364828806635 0.03454914825401289 0.016437844426928104 0.005928856172534874 0.001503341539029276 4.2714301545963265e-13; 2.1895179117406884e-15 0.0010411002857436305 0.003255743161934395 0.007659758643077384 0.013581240422684215 0.017260273812497083 0.01657647018690786 0.011616853176862237 0.005928856172534874 0.0023115783364469583 0.0006238849765532687 2.1895179117406884e-15; 1.9351535299019154e-18 0.00029141318932587976 0.0008608124023284637 0.0019037581454307546 0.0032124500970325837 0.003990509147237947 0.0038475611129833603 0.0027870615866269714 0.001503341539029276 0.0006238849765532687 0.00017702766887648104 1.9351535299019154e-18; 6.711882385110737e-21 2.4514117276084958e-17 1.5437604567532982e-14 1.6762365118986786e-12 3.138212957188894e-11 1.0130275516790658e-10 8.332937341028532e-11 1.4367785881450301e-11 4.2714301545963265e-13 2.1895179117406884e-15 1.9351535299019154e-18 6.711882385110737e-21]
    return xc, H
end

xc, H = init();


text:
abstract type FilterBank end

    struct Bark <: FilterBank end

    abstract type MelFilterBank <: FilterBank end
        struct Mel <: MelFilterBank end

        abstract type HTKMelFilterBank <: MelFilterBank end
            struct HTKMel <: HTKMelFilterBank end
            struct FCMel <: HTKMelFilterBank end

fb_type(::Bark) = "bark"
fb_type(::Mel) = "mel"
fb_type(::HTKMel) = "htkmel"
fb_type(::FCMel) = "fcmel"


text:
# Functions for getting parasitic assist from using other people's results.

if !isdefined(Main, :Response)
include("core.jl")
end

"""
    readResponses(filename)

Read responses from a file. Each line should be one of:
* a letter-coded response (b=black, y=yellow, g=green),
* a comment line beginning with #.
* a blank line
"""
function readResponses(filename::AbstractString)
    result = Response[]
    for line in eachline(filename)
        m = match(r"^[byg]{5}", line)
        if m ≠ nothing
            push!(result, Response(m.match))
            continue
        end
        m = match(r"^#.*$", line)
        m ≠ nothing && continue
        m = match(r"^ *$", line)
        m ≠ nothing && continue
        error("invalid line: $line")
    end
    return result
end

"""
    hasResponse(answer, response)

True if for given answer, there exists a guess that yields the given response.
This is a low-level function that operates on an AnswerIndex.
"""
function hasResponse(answer::AnswerIndex, response::Response)
    for g in 1:length(guessWords)
        if responseMatrix[answer,g] == response
            return true
        end
    end
    return false
end

"""
    sift(answers, response)

Return answers that have at least one guess that yields the given response.
This is a low-level routine that operates on Vector{AnswerIndex}.
"""
sift(answers::Vector{AnswerIndex}, response::Response) = filter(a->hasResponse(a, response), answers)

"""
    sift(answers, response)

Return answers that have at least one guess that yields the given response.
"""
function sift(answers::Vector{String}, response::String)
    return answerWords[sift(indicesOfWords(answers, answerWords), Response(response))]
end

"""
    sift(answers, responses)

Return answers that are plausible for the given responses.
"""
function sift(answers::Vector{String}, responses::Vector{Response})
    a = indicesOfWords(answers, answerWords)
    for r in responses
        a = sift(a, r)
    end
    return answerWords[a]
end

"""
    findErrors(answer, reponses)

Return responses that are inconsistent with the answer.
"""
function findErrors(answer::String, responses::Vector{Response})
    result = Response[]
    for r in responses
        if !hasResponse(indexOfWord(answer,answerWords), r)
            push!(result, r)
        end
    end
    return result
end

nothing

max_stars_repo_name: ArchRobison/WordleAnalysis

max_issues_repo_name: ArchRobison/WordleAnalysis

max_forks_repo_name: ArchRobison/WordleAnalysis


text:
ensurenode(g::Graph, i::Integer) = addnode!(g, i - nodecount(g));

"""
    load_ngraph(path, directed::Bool, [nodeids::Type])

Loads a graph from [ngraph binary format](https://github.com/anvaka/ngraph.tobinary#linksbin-format).
It requires a folder with 3 files in it: `links.bin`, `labels.json`, and `positions.bin`.
"""
function load_ngraph(path,
                     directed::Bool,
                     nodeids::Type=UInt32)
    g = Graph(directed=directed, nodeids=nodeids)
    # Load the links
    open(joinpath(path, "links.bin")) do f
        src = zero(nodeids)
        for i = reinterpret(nodeids, read(f))
            if i < 0
                src = -i
                ensurenode(g, src)
            else
                ensurenode(g, i)
                addedge!(g, src, i)
            end
        end
    end
    # Load the labels
    open(joinpath(path, "labels.json")) do f
        l = readline(f)[2:end - 1]
        g[:, :label] = [String(strip(i, [' ', '"'])) for i = split(l, ',')]
    end
    # Load the layout
    open(joinpath(path, "positions.bin")) do f
        ints = reinterpret(Int32, read(f))
        g[:, :x] = ints[1:3:end]
        g[:, :y] = ints[2:3:end]
        g[:, :z] = ints[3:3:end]
    end
    return g
end


text:
using NamedArrays, BenchmarkTools

I = [randstring(4) for i=1:100]
J = [randstring(4) for j=1:100]

d1 = Dict((i,j) => rand() for i in I, j in J)
d2 = Dict((i,j) => rand() for i in I, j in J)

d3 = Dict((i,j) => d1[i,j]*d2[i,j] for i in I, j in J)
d4 = Dict{Tuple{String,String}, Float64}()
for i in I, j in J
    d4[i,j] = d1[i,j]*d2[i,j]
end

println("\nDict comprehension\n", @benchmark d3 = Dict((i,j) => d1[i,j]*d2[i,j] for i in I, j in J))
println("\nDict assignment\n", @benchmark begin
    d4 = Dict{Tuple{String,String}, Float64}()
    for i in I, j in J
        d4[i,j] = d1[i,j]*d2[i,j]
    end
end)

n1 = NamedArray([rand() for i in I, j in J], (I,J))
n2 = NamedArray([rand() for i in I, j in J], (I,J))

n3 = NamedArray([n1[i,j]*n2[i,j] for i in I, j in J], (I,J))
n4 = NamedArray(zeros(100,100), (I,J))
for i in I, j in J
    n4[i,j] = n1[i,j]*n2[i,j]
end
println("\nNamedArray comprehension\n", @benchmark n3 = NamedArray([n1[i,j]*n2[i,j] for i in I, j in J], (I,J)))
println("\nNamedArray assignment\n", @benchmark begin
    n4 = NamedArray(zeros(100,100), (I,J))
    for i in I, j in J
        n4[i,j] = n1[i,j]*n2[i,j]
    end
end)

b1 = [rand() for i=1:100, j=1:100]
b2 = [rand() for i=1:100, j=1:100]

b3 = [b1[i,j]*b2[i,j] for i=1:100, j=1:100]
println("\nBasic Array comprehension (no name lookup)\n", @benchmark b3 = [b1[i,j]*b2[i,j] for i=1:100, j=1:100])

i, j = I[4], J[4]
println("\nBasic Dict comprehension (no name lookup)\n", @benchmark d3[i, j])
println("\nBasic Dict comprehension (no name lookup)\n", @benchmark d4[i, j])

println("\nNamedArray assignment\n", @benchmark n3[i, j])
println("\nNamedArray assignment\n", @benchmark n4[i, j])

println("\nBasic Array comprehension (no name lookup)\n", @benchmark b3[4, 4])


text:
module EBIC

export run_ebic

using JSON
using DataStructures: SortedSet
using ProgressMeter: next!, finish!, Progress
using ArgParse: ArgParseSettings, @add_arg_table, parse_args

include("constants.jl")
include("biclusterseval.jl")
include("evolution.jl")
include("scoring.jl")
include("algorithm.jl")

using .evolution: init_population, mutate
using .scoring: score_population
using .algorithm: update_rank_list!, ReverseOrdering
using .biclusterseval: get_biclusters, initialize_input_on_gpus

function run_ebic(;
    input_path::String,
    best_bclrs_stats = true,
    verbose = true,
    max_iterations = MAX_ITERATIONS,
    max_biclusters = MAX_BICLUSTERS_NUMBER,
    overlap_threshold = OVERLAP_THRESHOLD,
    negative_trends = NEGATIVE_TRENDS_ENABLED,
    approx_trends_ratio::Float32 = APPROX_TRENDS_RATIO,
    gpus_num = GPUS_NUMBER,
    output = false,
)
    data_load_time = @elapsed begin
        d_input_data = initialize_input_on_gpus(input_path, gpus_num)
    end
    @debug "Loading data to GPU took: $(data_load_time)s"

    # used to evaluate the iteration and timing of the best bclrs finding
    prev_top_bclrs = Vector()
    last_top_bclrs_change = (0, 0)

    if verbose
        p_bar = Progress(max_iterations, barlen = 20)
    end

    @debug "Starting algorithm..."
    # algorithm initialization steps
    start_time = time_ns()
    tabu_list = Set()
    tabu_hits = 0
    top_rank_list = SortedSet(Vector(), ReverseOrdering())

    cols_number = size(d_input_data[1], 2)

    old_population = init_population(cols_number, tabu_list)

    old_scored_population =
        score_population(d_input_data, old_population, gpus_num = gpus_num)

    update_rank_list!(top_rank_list, old_scored_population, overlap_threshold)

    i = 0
    while i < max_iterations
        i += 1
        new_population = Population()

        # reset penalties
        penalties = fill(1, cols_number)

        # elitism
        for (_, chromo) in top_rank_list
            !(length(new_population) < REPRODUCTION_SIZE) && break

            push!(new_population, chromo)
        end

        # perform mutations to replenish new population
        new_population, tabu_hits =
            mutate(new_population, old_scored_population, tabu_list, penalties, cols_number)

        # check if algorithm has found new solutions
        if tabu_hits >= MAX_NUMBER_OF_TABU_HITS
            verbose && finish!(p_bar)
            break
        end

        # evaluate fitness for new population
        new_scored_population =
            score_population(d_input_data, new_population, gpus_num = gpus_num)

        # save best chromosomes
        update_rank_list!(top_rank_list, old_scored_population, overlap_threshold)

        # proceed to the next generation
        old_scored_population = new_scored_population

        if best_bclrs_stats
            new_top_bclrs = collect(top_rank_list)[1:max_biclusters]
            if !isempty(prev_top_bclrs)
                changed = new_top_bclrs != prev_top_bclrs
                if changed
                    last_top_bclrs_change = (i, time_ns() - start_time)
                    prev_top_bclrs = new_top_bclrs
                end
            else
                prev_top_bclrs = new_top_bclrs
            end
        end

        verbose && next!(p_bar)
    end

    algorithm_time = time_ns() - start_time

    biclusters = get_biclusters(
        d_input_data,
        [last(p) for p in top_rank_list],
        gpus_num,
        negative_trends,
        approx_trends_ratio,
    )

    run_summary = Dict(
        "algorithm_time" => algorithm_time / 1e9,
        "data_load_time" => data_load_time,
        "biclusters" => biclusters[1:max_biclusters],
        "performed_iters" => i,
        "last_iter_tabu_hits" => tabu_hits,
    )

    if best_bclrs_stats
        run_summary["best_bclrs_iter"] = last_top_bclrs_change[1]
        run_summary["best_bclrs_time"] = last_top_bclrs_change[2] / 1e9
    end

    if output
        open("$(basename(input_path))-res.json", "w") do f
            JSON.print(f, run_summary["biclusters"])
        end
    end

    return run_summary
end

function main()
    args = ArgParseSettings("""
    EBIC is a next-generation biclustering algorithm based on artificial intelligence (AI). EBIC is probably the first algorithm capable of discovering the most challenging patterns (i.e. row-constant, column-constant, shift, scale, shift-scale and trend-preserving) in complex and noisy data with average accuracy of over 90%. It is also one of the very few parallel biclustering algorithms that use at least one graphics processing unit (GPU) and is ready for big-data challenges.
    """)

    @add_arg_table args begin
        "--input", "-i"
        help = "The path to the input file."
        arg_type = String
        default = INPUT_PATH
        dest_name = "input_path"

        "--max_iterations", "-n"
        help = "The maximum number of iterations of the algorithm."
        arg_type = Int
        default = MAX_ITERATIONS

        "--biclusters_num", "-b"
        help = "The number of biclusters that will be returned in the end."
        arg_type = Int
        default = MAX_BICLUSTERS_NUMBER
        dest_name = "max_biclusters"

        "--overlap_threshold", "-x"
        help = "The maximum similarity level of each two chromosomes held in top rank list."
        arg_type = Float64
        default = OVERLAP_THRESHOLD

        "--negative_trends", "-t"
        help = "Enable negative trends."
        action = :store_true

        "--gpus_num", "-g"
        help = "The number of gpus the algorithm should run on."
        arg_type = Int
        default = GPUS_NUMBER

        "--approx_trends", "-a"
        help = ""
        arg_type = Float32
        default = APPROX_TRENDS_RATIO
        dest_name = "approx_trends_ratio"

        "--verbose", "-v"
        help = "Turn on the progress bar."
        action = :store_true

        "--best_bclrs_stats", "-s"
        help = "Evaluate resulting biclusters finding iteration and time. Enabled, it slightly worsens overall algorithm performance."
        action = :store_true

        "--output", "-o"
        help = "Save biclusters to a file in the JSON format. The output file name is a concatenation of the input file name and '-res.json' suffix."
        action = :store_true

    end
    args = parse_args(args, as_symbols = true)

    results = run_ebic(; args...)

    JSON.print(results)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

end


text:
```Write a function that computes the list of the first 100 Fibonacci numbers.
The first two Fibonacci numbers are 1 and 1. The n+1-st Fibonacci number can be
computed by adding the n-th and the n-1-th Fibonacci number. The first few are
therefore 1, 1, 1+1=2, 1+2=3, 2+3=5, 3+5=8.```

using ResumableFunctions

@resumable function fibonnaci(n::Int) :: Int
  a = 0
  b = 1
  for i in 1:n-1
    @yield a
    a, b = b, a+b
  end
  a
end

function FibSeq(N::Int64)
    for fib in fibonnaci(N)
      println(fib)
    end
end

FibSeq(25)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using DemoPackageSSD
using Documenter

DocMeta.setdocmeta!(DemoPackageSSD, :DocTestSetup, :(using DemoPackageSSD); recursive=true)

makedocs(;
    modules=[DemoPackageSSD],
    authors="stephdoering <stephanie@thedoerings.net> and contributors",
    repo="https://github.com/stephdoering/DemoPackageSSD.jl/blob/{commit}{path}#{line}",
    sitename="DemoPackageSSD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://stephdoering.github.io/DemoPackageSSD.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/stephdoering/DemoPackageSSD.jl",
    devbranch="master",
)


text:
#a module for a web crawler. sort out comments and stuff in a bit

module WebCrawler

export crawl, begin_crawl. setup_crawler

include("crawler.jl")

end #module

max_stars_repo_path: src/Web Crawler.jl

max_issues_repo_path: src/Web Crawler.jl

max_forks_repo_path: src/Web Crawler.jl


text:
# This program computes Feigenbaum's alpha constant using
# Feigenbaum's functional equation derived via renormalization
# group methods.
# Oct 2017 -- Stuart Brorson, sdb@cloud9.net

using ForwardDiff

# I define the helper fcns first.

#-------------------------------------------------------
# Define universal fcn g(z, a).  It is expressed as a power
# series, coefficients a.  Only even powers are needed since
# g(z) is even.
function g(z, a::Vector)
  #println("---> Entered g\n")
  s = 0;
  for i=length(a):-1:1
    # println("i = $i, s = $s\n")
    s = z*z*(s + a[i]);
  end
  #println("<--- Leaving g\n")
  return (1+s)
end

#---------------------------------------------------------
# Define relation obeyed by g(z).  Name it f(z, a)
function f(z, a::Vector)
  #println("---> Entered f\n")
  const one = big(1.0)
  alpha = one/g(one, a);
  r = g(z, a) - alpha*g(g(z/alpha, a), a);
  #println("<--- Leaving g\n")
  return r
end

#--------------------------------------------------------
# Define fcn which takes gradient of f w.r.t. a.
# Returns row vector which is gradient of f w.r.t. a,
# evaluated at z, a.
function gradf(z, a::Vector)
  #println("---> Entered gradf\n")
  # Specialize to position z.  Result is fcn of coeffs a only.
  fbeta(a) = f(z, a);

  # Compute gradient w.r.t. beta and return it
  gf1 = ForwardDiff.gradient(fbeta, a)
  #println("<--- Leaving gradf\n")
  return gf1'
end

# ------------------------------------------------------------
# The main event -- a computation of alpha using Newton's method.
function compute_alpha(Numdigs, N, betain)

  # First set up variables, including BigFloat settings
  # Setting precision to 5x Numdigs seems to keep Newton's method
  # from wandering away.  Smaller values don't.  More fidding might
  # improve things here.
  setprecision(5*Numdigs);

  # Set stopping tol for Newton's method so I get the number of digits
  # I want.
  const tol = BigFloat(10.0)^-(Numdigs);

  const zi = big.(linspace(1/N, 1, N));   # Grid of sample points zi
  fn = zeros(BigFloat, N);             # Function f vector
  Jn = zeros(BigFloat, N, N);        # Jacobian

  # Need expansion coeffs beta.  Start with the one passed in, betain.
  # N should always be larger than length(betain), so we are adding to the
  # betas we computed last time.
  betan = zeros(BigFloat, N);
  for i = 1:length(betain)
    betan[i] = betain[i];
  end

  # println("-----------------------------\n")
  sttim = now();  # Keep track of loop timing.
  # Now enter Newton loop.
  for cnt = 1:50
    nowtim = now();  # Keep track of loop timing
    println("Newton iteration = $cnt, time elapsed = $(nowtim-sttim)\n")

    # Compute new f and Jacobian upon each iteration
    for i = 1:N
      fn[i] = f(zi[i], betan);
      Jn[i,:] = gradf(zi[i], betan);
    end

    # Now compute step to take, then take it.
    sn = Jn\fn;
    betanp1 = betan - sn;

    # Check for convergence
    if (norm(sn) < tol)
      println("\nDone!  Converged after $cnt iterations.\n")
      # compute alpha
      myalpha = 1/g(1, betanp1);
      return myalpha, betanp1;
    end

    # Move values back
    betan = betanp1

  end
  println("!!!!!!!  Didn't converge!  !!!!!!!!\n")

end

#--------------------------------------------------
function iterate_alpha()
  # This is a runner -- call this fcn.  This fcn calls
  # compute_alpha asking for different number of digits
  # at different precision/tolerance levels.
  # It calls with the beta coefficients computed in the
  # last iteration as the starting point to help Newton's
  # method from wandering away.

  # Seed betan with 3 values.
  betan = zeros(BigFloat, 3);
  betan[1] = big(-1.52763883147);
  betan[2] = big(1.048327004372e-1);
  betan[3] = big(2.669121419134012e-2);

  # Step up number of pts to compute on.  Change this to get
  # more digits
  for N = 10:20:1700
    println("=================================\n")
    println("N = $N\n")
    Numdigs = Int(floor(2*N));

    # Call the function which computes alpha at this level
    # of precision.
    myalpha, betan = compute_alpha(Numdigs, N, betan)

    # Print out myalpha to N digits of precision.  Must add 3
    # to account for negative sign and "2."
    a = string(myalpha);
    println("myalpha = $(a[1:N+3])\n")

    f = open("alpha.dat", "a");
    write(f, "==========\nN = $N\n$(a[1:N+3])\n")
    close(f)

    # Compare against the version from
    # http://www.plouffe.fr/simon/constants/feigenbaum.txt
    oeisalpha = oeis_alpha(3*Numdigs);
    diff = Float64(abs(myalpha) - abs(oeisalpha));
    println("diff = $diff\n")
  end

end

#------------------------------------------------
function oeis_alpha(Numdigs)
  # Return up to 1018 digits of alpha from
  # http://www.plouffe.fr/simon/constants/feigenbaum.txt

  if (Numdigs <= 1018)
    setprecision(Numdigs);
  else
    setprecision(1018);
  end

  str = "-2."
  str = str*"5029078750958928222839028732182157863812713767271499773361920567"
  str = str*"7923546317959020670329964974643383412959523186999585472394218237"
  str = str*"7785445179272863314993372578112163594879503744781260997380598671"
  str = str*"2397117373289276654044010306698313834600094139322364490657889951"
  str = str*"2205843172507873377463087853424285351988587500042358246918740820"
  str = str*"4281700901714823051821621619413199856066129382742649709844084470"
  str = str*"1008054549677936760888126446406885181552709324007542506497157047"
  str = str*"0475419932831783645332562415378693957125097066387979492654623137"
  str = str*"6745918909813116752434221110130913127837160951158341230841503716"
  str = str*"4997020224681219644081216686527458043026245782561067150138521821"
  str = str*"6449532543349873487413352795815351016583605455763513276501810781"
  str = str*"1948369459574850237398235452625632779475397269902012891516645793"
  str = str*"9420198920248803394051699686551494477396533876979741232354061781"
  str = str*"9896112494095990353128997733611849847377946108428833293833903950"
  str = str*"9008914086351525626803381414669279913310743349705143545201344643"
  str = str*"4264752001621384610729922641994332772918977769053802596851"

  alpha = parse(BigFloat, str)

end

max_stars_repo_name: brorson/FeigenbaumConstants

max_issues_repo_name: brorson/FeigenbaumConstants

max_forks_repo_name: brorson/FeigenbaumConstants


text:
using SequentialMonteCarlo
using RNGPool
using CoupledConditionalSMC
import SMCExamples.Particles.Float64Particle

@inline function simplelM(::Int64, particle::Float64Particle,
  newParticle::Float64Particle, ::Nothing)
  x::Float64 = particle.x
  y::Float64 = newParticle.x
  t::Float64 = x - y
  return -0.5 * t * t
end

function setupSimpleModel(n::Int64, boxSize::Float64)
  @inline function lG(::Int64, particle::Float64Particle, ::Nothing)
    return abs(particle.x) < boxSize ? 0.0 : -Inf
  end
  @inline function M!(newParticle::Float64Particle, rng::RNG, p::Int64,
    particle::Float64Particle, ::Nothing)
    if p == 1
      newParticle.x = randn(rng)
    else
      newParticle.x = particle.x + randn(rng)
    end
  end
  return SMCModel(M!, lG, n, Float64Particle, Nothing), simplelM
end

model, lM = setupSimpleModel(10000, 5.0)

import Statistics: mean, std

setRNGs(12345)

m = 1000

using DataFrames
box_at_df = DataFrame(n=Int[], N=Int[], mean=Union{Float64, Missing}[],
  std=Union{Float64, Missing}[])

ns = 200:200:1600
Ns = round.(Int64, 100*(ns/100).^(1.4))

for i in 1:length(ns)
  n = ns[i]
  N = Ns[i]
  println(n, " ", N)
  vs = CoupledConditionalSMC.couplingTimes(model, N, n, m, true, true)
  println(mean(vs))
  push!(box_at_df, [n, N, mean(vs), std(vs)])
end

using JLD2
@save "paper/box_at.jld2" box_at_df

println(box_at_df)

box_all_df = DataFrame(n=Int[], N=Int[], type=String[], mean=Union{Float64, Missing}[],
  std=Union{Float64, Missing}[])

ns = 200:200:1600
Ns = round.(Int64, 128*(ns/200))

for i in 1:length(ns)
  n = ns[i]
  N = Ns[i]
  println(n, " ", N)
  vs = CoupledConditionalSMC.couplingTimes(model, N, n, m, true, true)
  println(mean(vs))
  push!(box_all_df, [n, N, "AT", mean(vs), std(vs)])
  vs = CoupledConditionalSMC.couplingTimes(model, lM, N, n, m, :AS, true, true)
  println(mean(vs))
  push!(box_all_df, [n, N, "AS", mean(vs), std(vs)])
  vs = CoupledConditionalSMC.couplingTimes(model, lM, 128, n, m, :BS, true, true)
  println(mean(vs))
  push!(box_all_df, [n, 128, "BS", mean(vs), std(vs)])
end

using JLD2
@save "paper/box_all.jld2" box_all_df

println(box_all_df)


text:
# This is a Node in the Binning Analysis tree that averages two values. There
# is one of these for each binning level. When two values should be
# compressed, this is done immediately, so that only one value needs to be saved.
# switch indicates whether value should be written to or averaging should happen.
mutable struct Compressor{T}
    value::T
    switch::Bool
end

struct LogBinner{N, T}
    # list of Compressors, one per level
    compressors::NTuple{N, Compressor{T}}

    # sum(x) for all values on a given lvl
    x_sum::Vector{T}
    # sum(x.^2) for all values on a given lvl
    x2_sum::Vector{T}
    # number of values that are summed on a given lvl
    count::Vector{Int64}
end

"""
    LogBinner([T = Float64, N = 32])

Creates a Binning Analysis which can take 2^N-1 value sof type T.
"""
LogBinner(T::Type = Float64, N::Int64 = 32) = LogBinner(zero(T), N)

# TODO
# Currently, the Binning Analysis requires a "zero" to initialize x_sum and
# x2_sum. It's not necessary for the Compressors.
# Since Arrays do not have a static size, we cannot generate a fitting zero
# automatically. So currently we let the user supply it.
# Other possibilities:
# > generate x_sum, x2_sum with #undef
# bad: requires many checks (if isassigned ... else ...)
# > force StaticArrays (and/or tuples)
# good: faster for small arrays
# bad: requires user to use another package
# > initialize is first push!
# bad: also requires frequent checks (if first push ... else ...)
# ...?
"""
    LogBinner([zero = 0.0, N = 32])

Creates a new Binning Analysis which can take 2^N-1 values of type T. The type
is inherited by the given zero. Returns a Binning Analysis object.

Values can be added using `push!(LogBinner, value)`.
"""
function LogBinner(_zero::T = zero(Float64), N::Int64 = 32) where {T}
    # heuristic to set sum type (#2)
    D = ndims(_zero)
    S = if eltype(T)<:Real
        D > 0 ? Array{Float64, D} : Float64
    else
        D > 0 ? Array{ComplexF64, D} : ComplexF64
    end

    LogBinner{N, S}(
        tuple([Compressor{S}(copy(_zero), false) for i in 1:N]...),
        [copy(_zero) for _ in 1:N],
        [copy(_zero) for _ in 1:N],
        zeros(Int64, N)
    )
end

# TODO typing?
"""
    append!(LogBinner, values)

Adds an array of values to the Binning Analysis by applying push! to each
element.
"""
function append!(B::LogBinner, values::AbstractArray)
    for value in values
        push!(B, value)
    end
    nothing
end

"""
    push!(LogBinner, value)

Pushes a new value into the Binning Analysis.
"""
function push!(B::LogBinner{N, T}, value::S) where {N, T, S}
    ndims(T) == ndims(S) || throw(DimensionMismatch("Expected $(ndims(T)) dimensions but got $(ndims(S))."))

    _push!(B, 1, value)
end

_square(x) = x^2
_square(x::Complex) = Complex(real(x)^2, imag(x)^2)
_square(x::AbstractArray) = _square.(x)

# recursion, back-end function
function _push!(B::LogBinner{N, T}, lvl::Int64, value::S) where {N, T <: Number, S}
    C = B.compressors[lvl]

    # any value propagating through this function is new to lvl. Therefore we
    # add it to the sums. Note that values pushed to the output arrays are not
    # added here until the array drops to the next level. (New compressors are
    # added)
    B.x_sum[lvl] += value
    B.x2_sum[lvl] += _square(value)
    B.count[lvl] += 1

    if !C.switch
        # Compressor has space -> save value
        C.value = value
        C.switch = true
        return nothing
    else
        # Do averaging
        if lvl == N
            # No more propagation possible -> throw error
            throw(OverflowError("The Binning Analysis ha exceeddd its maximum capacity."))
        else
            # propagate to next lvl
            C.switch = false
            _push!(B, lvl+1, 0.5 * (C.value + value))
            return nothing
        end
    end
    return nothing
end

function _push!(
        B::LogBinner{N, T},
        lvl::Int64,
        value::S
    ) where {N, T <: AbstractArray, S}

    C = B.compressors[lvl]
    B.x_sum[lvl] .+= value
    B.x2_sum[lvl] .+= _square(value)
    B.count[lvl] += 1

    if !C.switch
        C.value = value
        C.switch = true
        return nothing
    else
        if lvl == N
            throw(OverflowError("The Binning Analysis ha exceeddd its maximum capacity."))
        else
            C.switch = false
            _push!(B, lvl+1, 0.5 * (C.value .+ value))
            return nothing
        end
    end
    return nothing
end


text:
function Example_4_3(; T = 200, seed = 1)

    # Set random number generator

    rg = MersenneTwister(seed)

    # state x is four dimensional, comprises car position and velocities

    Δt  = 1 / 10
    σ₁  = σ₂  = 1 / 2
    qᶜ₁ = qᶜ₂ = 1

    # Matrices in the dynamic model are:

    A = [1  0  Δt  0;
         0  1   0  Δt;
         0  0   1  0;
         0  0   0  1]

    Q = [(qᶜ₁*Δt^3)/3         0        (qᶜ₁*Δt^2)/2       0        ;
              0         (qᶜ₂*Δt^3)/3         0        (qᶜ₁*Δt^2)/2 ;
         (qᶜ₁*Δt^2)/2         0            qᶜ₁*Δt         0        ;
              0         (qᶜ₂*Δt^2)/2         0          qᶜ₂*Δt     ]

    # Measurement matrices

    H = [1 0 0 0;
         0 1 0 0]

    R = [σ₁^2  0;
         0    σ₂^2]

    # initialise

    m₀, P₀ = zeros(4), 0.1*Matrix(I, 4, 4)

    x₀ = rand(MvNormal(m₀, P₀))

    # store measurements and states

    yclean = Array{Vector{Float64},1}(undef, T)

    y = Array{Vector{Float64},1}(undef, T)

    x = Array{Vector{Float64},1}(undef, T)

    # first state

    x[1] = A * x₀

    y[1] = rand(rg, MvNormal(H * x[1], R))

    for k in 2:T

        x[k] = rand(rg, MvNormal(A * x[k-1], Q)) # (4.32)

        yclean[k] = H * x[k]

        y[k] = rand(rg, MvNormal(H * x[k], R)) # (4.33)

    end

    return yclean, y, x, A, H, Q, R, m₀, P₀

end

function runExample_4_3(seed = 1)

    # Simulate data

    yclean, y, x, A, H, Q, R, m₀, P₀ = Example_4_3(seed=seed)

    # Inference

    μfilter, Σfilter = filteringrecursion(y; A = A, H = H, Q = Q, R = R, m₀ = m₀, P₀ = P₀)

    μsmooth, Σsmooth = smoothingrecursion(y; A = A, H = H, Q = Q, R = R, m = μfilter, P = Σfilter)

    # collect positions

    truepositions = reduce(hcat, [xₖ[1:2] for xₖ in x])

    filtpositions = reduce(hcat, [xₖ[1:2] for xₖ in μfilter])

    smoopositions = reduce(hcat, [xₖ[1:2] for xₖ in μsmooth])

    # plot positions

    figure(1)

    subplot(211); cla()

    plot(truepositions[1,:], truepositions[2,:], "o", label = "true positions")

    plot(filtpositions[1,:], filtpositions[2,:], ".", label = "filter positions", alpha=0.5)

    legend()

    subplot(212); cla()

    plot(truepositions[1,:], truepositions[2,:], "o", label = "true positions")

    plot(smoopositions[1,:], smoopositions[2,:], ".", label = "smooth positions", alpha=0.5)

    legend()

    nothing

end


text:
using Distributions
using Cytof5
using Flux
using BSON

# Plotting
using PyCall
matplotlib = pyimport("matplotlib")
plt = pyimport("matplotlib.pyplot")
matplotlib.use("Agg")
# Load current dir
pushfirst!(PyVector(pyimport("sys")."path"), "../vb")
plot_yz = pyimport("plot_yz")
blue2red = pyimport("blue2red")

# Directory containing CB results
data_dir = "../../sims/vb/results/vb-cb-paper/test/"

# Path to mm1 and best output
path_to_vb_cb_output = "$(data_dir)/output.bson"

# Get output for miss-mech-0 and miss-mech-1
output = BSON.load(path_to_vb_cb_output)

# Posterior Samples
samps = [Cytof5.VB.rsample(output[:state])[2] for _ in 1:100]
sig2s = [s.sig2 for s in samps]
Zs = [reshape(s.v, 1, length(s.v)) .> s.H for s in samps]
Ws = cat([s.W for s in samps]..., dims=3)


text:
module HackRvaJulia

include("introduction.jl")

end # module


text:
quatErr(q_err::Quat{TD}) where {TD} = (0 <= q_err.w) ? SVector( q_err.x,  q_err.y,  q_err.z) : SVector(-q_err.x, -q_err.y, -q_err.z)

function quatErr(q1::Quat{TD}, qRef::Quat{TD}) where {TD}
  q_err = q1 * inv(qRef)
  return quatErr(q_err)
end

cheapRV(q::Quat) = 2 * quatErr(q)
cheapRV(spq::SPQuat) = cheapRV(Quat(spq))

components(q::Quat{T}) where {T} = SVector{4,T}(q.w, q.x, q.y, q.z)
components(spq::SPQuat{T}) where {T} = SVector{3,T}(spq.x, spq.y, spq.z)


text:
#=
day17:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-17
=#

using Test, AoC

function launch(v, target)
    p = [0,0]
    maxHeight = p[2]
    while (p[1] <= target[1,2]) && (p[2] >= target[2,1])
        if p[2] > maxHeight
            maxHeight = p[2]
        end
        if (target[1,1] <= p[1] <= target[1,2]) && (target[2,1] <= p[2] <= target[2,2])
            return maxHeight
        end
        p = p + v
        v[2] -= 1
        v[1] -= sign(v[1])
    end
    return false
end
@test launch([17,-4], [[20, -10] [30, -5]]) == false
@test launch([7,2], [[20, -10] [30, -5]]) == 3
@test launch([6,3], [[20, -10] [30, -5]]) == 6
@test launch([6,9], [[20, -10] [30, -5]]) == 45
@test launch([23,-10], [[20, -10] [30, -5]]) == 0

exampleLines(17,1) .|> line -> @test launch(parse.(Int, split(line, ",")), [[20, -10] [30, -5]]) !== false

function part1(target)
    maxHeight = 0
    for v_x ∈ 1:target[1,2]
        for v_y ∈ 1:abs(target[2,1])+1
            height = launch([v_x,v_y], target)
            if height > maxHeight
                maxHeight = height
            end
        end
    end
    return maxHeight
end

function part2(target)
    hits = 0
    for v_x ∈ 1:target[1,2]
        for v_y ∈ -abs(target[2,1])-1:abs(target[2,1])+1
            height = launch([v_x,v_y], target)
            if height !== false
                hits += 1
            end
        end
    end
    return hits
end

@test part1([[20, -10] [30, -5]]) == 45
@test part2([[20, -10] [30, -5]]) == 112

@show @time part1([[96, -144] [125,-98]])
@show @time part2([[96, -144] [125,-98]])


text:
# See also test_pngsuite.jl file

@testset "Images with background" begin
    img_no_bg = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=false)
    img_file_bg = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=true)
    img_index_bg = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=UInt8(46), expand_paletted=true)
    img_gray_bg = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=Gray{N0f8}(1.0))
    img_rgb_bg = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=RGB{N0f8}(1.0))
    img_gray_bg16 = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=Gray{N0f16}(1.0))
    img_rgb_bg16 = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=RGB{N0f16}(1.0))

    @testset "load" begin
        @testset "Background is not ignored" begin
            @test img_file_bg[1, 1] == RGBA{N0f8}(1.0, 1.0, 1.0, 1.0)
            @test img_no_bg[1, 1] == RGBA{N0f8}(0.0, 0.0, 0.0, 0.0)
        end

        @testset "Background as index into palette" begin
            @test img_file_bg ≈ img_index_bg
        end
        @testset "Background as RGB 8 bit" begin
            @test img_file_bg ≈ img_rgb_bg
        end
        @testset "Background as Gray 8 bit" begin
            @test img_file_bg ≈ img_gray_bg
        end
        @testset "Background as RGB 16 bit" begin
            @test img_file_bg ≈ img_rgb_bg16
        end
        @testset "Background as Gray 16 bit" begin
            @test img_file_bg ≈ img_gray_bg16
        end
    end

    img_no_bg_no_gamma = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=false, gamma=1.0)
    img_file_bg_no_gamma = PNGFiles.load("test_images/paletted_implicit_transparency.png", background=true, gamma=1.0)
    new_path = "test_images/paletted_implicit_transparency_bg_new.png"
    @testset "Paletted background rountripping" begin
        PNGFiles.save(new_path, img_no_bg_no_gamma, background=UInt8(46), file_gamma=1.0)
        new_img = PNGFiles.load(new_path, background=true, gamma=1.0)
        @test img_file_bg_no_gamma ≈ new_img
    end
end


text:
function localSearch2(estimator::SBMEstimator, dataset::Dataset)::Tuple{SBM, Matrix{Int}, OptResults}
    """ Local Search 2
    Initialization with random assignments x.
    This local search heuristic works in
    a solution representation-decoder scheme
    The local search (i.e. the search for improving moves) is done
    in the space of assignments x. However, whenever a solution
    is to be evaluated the complete solution is first retrieved,
    by finding the optimal value of w analytically.

    Parameters
    ----------
    estimator : SBMEstimator
        Optimization method specifications.
    dataset : Dataset
        Dataset representing an observed graph.

    Returns
    -------
    sbm : SBM
        Stochastic Block Model
    x : Array{Int64,2}
        Matrix of assignments of nodes to groups
    opt_results : OptResults
        Results of the optimization process
    """
    n = dataset.n
    q = dataset.n_communities

    timeLimit = estimator.time_limit

    start = time(); # Start counting time

    # Initialize with random assignments
    x = randomAssignments(dataset, seed=estimator.seed)
    w = optimalProbMatrix(dataset, x)

    # Calculate objective value
    obj_value = calculateObjective(dataset, w, x)
    best_solution = (copy(w), copy(x))
    best_obj = obj_value
    best_move = nothing

    status = nothing
    iterations = 0
    improved = true
    while true
        ## Base cases
        # Break if there is no improvement
        if ~improved
            status = :LocalOptimum_LS2
            break
        end
        improved = false

        # Break if time limit is exceeded
        availableTime = timeLimit - (time() - start)
        if availableTime <= 0
            status = :UserLimit
            break
        end
        iterations += 1

        # Find improving relocation move
        improved, best_obj, w, x = findImprovingRelocation(dataset, w, x, true, best_imp=true, time_limit=availableTime)
    end
    solvetime = time() - start;

    # SBM
    sbm = SBM(w, "poisson") # assumes a poisson distribution

    # Build OptResults
    obj_lb = -Inf
    obj_ub = best_obj
    nodecount = nothing
    lazycount = nothing
    opt_results = OptResults(obj_lb, obj_ub, status, solvetime, iterations, nodecount, lazycount)
    if estimator.verbose
        displayResults(opt_results)
    end

    return sbm, x, opt_results
end


text:
function render(text::String, cursor::TextCursor, font::GLFont, displayableAray::Shape)
	for char in text
        if char == '\t'
        	cursor.x += font.properties.advance * 3
        elseif char == '\r' || char == '\n'
        	cursor.y -= font.properties.lineHeight
         	cursor.x = cursor.intend
        else
            if (cursor.x <= displayableAray.x + displayableAray.w && cursor.y  <= displayableAray.y + displayableAray.h - font.properties.lineHeight) &&
            	(cursor.y  >= displayableAray.y )
                render("charOffset", [cursor.x, cursor.y], font.gl.program.id)
                #in the GLBuffer for the text triangle data, there are six verts per char
    			glDrawArrays(GL_TRIANGLES, int(char) * 6, 6)
            end
            cursor.x += font.properties.advance
        end
    end
end

global selectionStyle = Dict{ASCIIString, Any}(["textColor" => Float32[0,0,0,0], "backgroundColor" => Float32[0.02,0.8,0.01,0.2]])

function render(t::TextField, font::GLFont)
	displayableAray = t.area
	render(font.gl)
	#render selection
	if(t.hasFocus)
		startLine, line = findline(t.newLineIndexes, first(t.selection))
		xPosition = first(t.selection) - first(line)
		cursorX = t.x + (xPosition * font.properties.advance)
		cursorY = t.y - ((startLine-1) * font.properties.lineHeight)
		selectionStart =  TextCursor(cursorX, cursorY, cursorX)
		render(selectionStyle, font.gl.program.id)
		if length(t.selection) < 1
			render("|", TextCursor(cursorX + (font.properties.advance / 2f0), cursorY, t.x), font, displayableAray)
		else
			#naive approach, just render selection as a blank string with backgroundcolor
			a = chr2ind(t.text, max(first(t.selection), 1))
			b = chr2ind(t.text, min(last(t.selection), length(t.text)))
			selectionString = map(x -> (x == '\r' || x == '\n') ? x : ' ', t.text[a:b])
			render(selectionString, TextCursor(cursorX, cursorY, t.x), font, displayableAray)
		end
	end

	#render text
	startCursor = TextCursor(t.x, t.y, t.x)
	if !isempty(t.text)
		for elem in t.styles

				# assert range is in text range
				segment = intersect(elem.segment, 1:(length(t.text) + 1))
				if first(segment) > last(segment)
					error("segment not in textrange! segment: ", segment, " length(text): ", length(t.text), " stylesegment: ", elem.segment)
				end

				render(elem.style, font.gl.program.id)
				a = chr2ind(t.text, first(segment))
				b = chr2ind(t.text, last(segment))
				render(t.text[a:b], startCursor, font, displayableAray)
		end
	end
end

function render_to_stencil_buffer(mask)

	glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
	glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE)
	glEnable(GL_STENCIL_TEST)
	glStencilFunc(GL_ALWAYS, 1, 1)
	glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
	glDisable(GL_DEPTH_TEST)
	glDisable(GL_TEXTURE_2D)

	render(mask)

	glEnable(GL_TEXTURE_2D)
	glEnable(GL_DEPTH_TEST)
	glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)
	glStencilFunc(GL_EQUAL, 1, 1)
	glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP)

end


text:
function sample_points(N_domain, N_boundary)
    # interior nodes
    X_domain = rand(Float64,(N_domain, 2))
    X_boundary = zeros((N_boundary, 2))

    # generate random boundary points on the faces of the domain
    N_bd_each=convert(Int64, N_boundary/4)
    # bottom face
    X_boundary[1:N_bd_each, 1] = rand(Float64,(N_bd_each,1))
    # right face
    X_boundary[N_bd_each+1:2*N_bd_each, 1] .= 1
    X_boundary[N_bd_each+1:2*N_bd_each, 2] = rand(Float64,(N_bd_each,1))
    # top face
    X_boundary[2*N_bd_each+1:3*N_bd_each, 1] = rand(Float64,(N_bd_each,1))
    X_boundary[2*N_bd_each+1:3*N_bd_each, 2] .= 1
    # left face
    X_boundary[3*N_bd_each+1:N_boundary, 2] = rand(Float64,(N_bd_each,1))
    return X_domain, X_boundary
end


text:
import ViscousFlow: Fields
using FFTW

struct Not{T}
  idx::T
end
import Base: to_indices, uncolon, tail, _maybetail

@inline to_indices(A, inds, I::Tuple{Not, Vararg{Any}}) =
   (setdiff(uncolon(inds, (:, tail(I)...)), I[1].idx), to_indices(A, _maybetail(inds), tail(I))...)

@testset "Grid Routines" begin

  # size
  nx = 12; ny = 12

  # sample point
  i = 5; j = 7
  cellzero = Nodes(Dual,(nx,ny))
  nodezero = Nodes(Primal,cellzero)
  facezero = Edges(Primal,cellzero)
  dualfacezero = Edges(Dual,cellzero)

  @test typeof(cellzero) <: ScalarGridData

  @test typeof(nodezero) <: ScalarGridData

  @test typeof(facezero) <: VectorGridData

  @test typeof(dualfacezero) <: VectorGridData

  @test typeof(facezero.u) <: XEdges{Primal}
  @test typeof(facezero.v) <: YEdges{Primal}

  @test typeof(dualfacezero.u) <: XEdges{Dual}
  @test typeof(dualfacezero.v) <: YEdges{Dual}

  cellunit = deepcopy(cellzero)
  cellunit[i,j] = 1.0

  nodeunit = deepcopy(nodezero)
  nodeunit[i,j] = 1.0

  facexunit = deepcopy(facezero)
  facexunit.u[i,j] = 1.0

  faceyunit = deepcopy(facezero)
  faceyunit.v[i,j] = 1.0

  dualfacexunit = deepcopy(dualfacezero)
  dualfacexunit.u[i,j] = 1.0

  dualfaceyunit = deepcopy(dualfacezero)
  dualfaceyunit.v[i,j] = 1.0

  @testset "Basic array operations" begin
    w = zero(cellunit)
    w .= cellunit
    @test w[i,j] == 1.0
    q = similar(facexunit)
    q .= facexunit
    @test q.u[i,j] == 1.0
    @test iszero(q.v)
  end

  @testset "Inner products and norms" begin
    w = zero(cellunit)
    i0, j0 = rand(2:nx-1), rand(2:ny-1)
    w[i0,j0] = 1.0
    @test Fields.norm(w)*sqrt((nx-2)*(ny-2)) == 1.0
    w .= 1.0
    @test Fields.norm(w) == 1.0

    @test 2*w == w*2

    p = Nodes(Primal,w)
    p .= 1.0
    @test Fields.norm(p) == 1.0
    p2 = deepcopy(p)
    @test Fields.dot(p,p2) == 1.0
    @test Fields.norm(p-p2) == 0.0

    q = Edges(Dual,w)
    q.u .= 1.0
    q2 = deepcopy(q)
    @test Fields.dot(q,q2) == 1.0

    q = Edges(Primal,w)
    q.u .= 1.0
    q2 = deepcopy(q)
    @test Fields.dot(q,q2) == 1.0

    @test Fields.integrate(w) == 1.0

    @test Fields.integrate(p) == 1.0

    q .= 1
    @test Fields.norm(2*q) == sqrt(8)

  end

  @testset "Dual cell center data Laplacian" begin
    lapcell = laplacian(cellunit)
    @test lapcell[i,j] == -4.0
    lapcell[i,j] = 0.0
    @test lapcell[i+1,j] == lapcell[i-1,j] == lapcell[i,j-1] == lapcell[i,j+1] == 1.0
    lapcell[i+1,j] = lapcell[i-1,j] = lapcell[i,j-1] = lapcell[i,j+1] = 0.0
    @test iszero(lapcell)
  end

  @testset "Dual cell center data curl" begin
    q = curl(cellunit)
    @test q.u[i,j-1] == 1.0 && q.u[i,j] == -1.0
    q.u[i,j-1] = q.u[i,j] = 0.0
    @test iszero(q.u)
    @test q.v[i-1,j] == -1.0 && q.v[i,j] == 1.0
    q.v[i-1,j] = q.v[i,j] = 0.0
    @test iszero(q.v)
  end

  @testset "Dual cell node gradient" begin
    q = grad(nodeunit)
    @test q.u[i,j] == 1.0 && q.u[i+1,j] == -1.0
    q.u[i,j] = q.u[i+1,j] = 0.0
    @test iszero(q.u)
    @test q.v[i,j] == 1.0 && q.v[i,j+1] == -1.0
    q.v[i,j] = q.v[i,j+1] = 0.0
    @test iszero(q.v)
  end

  @testset "Face data curl" begin
    cellcurl = curl(facexunit)
    @test cellcurl[i,j] == -1.0 && cellcurl[i,j+1] == 1.0
    cellcurl[i,j] = cellcurl[i,j+1] = 0.0
    @test iszero(cellcurl)
    cellcurl = curl(faceyunit)
    @test cellcurl[i,j] == 1.0 && cellcurl[i+1,j] == -1.0
    cellcurl[i,j] = cellcurl[i+1,j] = 0.0
    @test iszero(cellcurl)
  end

  @testset "Face data divergence" begin
    nodediv = divergence(facexunit)
    @test nodediv[i,j] == -1.0 && nodediv[i-1,j] == 1.0
    nodediv[i,j] = nodediv[i-1,j] = 0.0
    @test iszero(nodediv)
    nodediv = divergence(faceyunit)
    @test nodediv[i,j] == -1.0 && nodediv[i,j-1] == 1.0
    nodediv[i,j] = nodediv[i,j-1] = 0.0
    @test iszero(nodediv)
  end

  @testset "Face data Laplacian" begin
    lap = laplacian(facexunit)
    @test lap.u[i,j] == -4.0
    lap.u[i,j] = 0.0
    @test lap.u[i+1,j] == lap.u[i-1,j] == lap.u[i,j-1] == lap.u[i,j+1] == 1.0
    lap.u[i+1,j] = lap.u[i-1,j] = lap.u[i,j-1] = lap.u[i,j+1] = 0.0
    @test iszero(lap.u)
    @test iszero(lap.v)

    lap = laplacian(faceyunit)
    @test lap.v[i,j] == -4.0
    lap.v[i,j] = 0.0
    @test lap.v[i+1,j] == lap.v[i-1,j] == lap.v[i,j-1] == lap.v[i,j+1] == 1.0
    lap.v[i+1,j] = lap.v[i-1,j] = lap.v[i,j-1] = lap.v[i,j+1] = 0.0
    @test iszero(lap.u)
    @test iszero(lap.v)
  end

  @testset "Dual face data divergence" begin
    celldiv = divergence(dualfacexunit)
    @test celldiv[i,j] == 1.0 && celldiv[i+1,j] == -1.0
    celldiv[i,j] = celldiv[i+1,j] = 0.0
    @test iszero(celldiv)
    celldiv = divergence(dualfaceyunit)
    @test celldiv[i,j] == 1.0 && celldiv[i,j+1] == -1.0
    celldiv[i,j] = celldiv[i,j+1] = 0.0
    @test iszero(celldiv)
  end

  @testset "Face data shift to dual face" begin
    shiftx = Edges(Dual,cellzero)
    interpolate!(shiftx,facexunit)
    @test shiftx.u[i,j] == shiftx.u[i-1,j] == shiftx.u[i,j+1] == shiftx.u[i-1,j+1] == 0.25
    shiftx.u[i,j] = shiftx.u[i-1,j] = shiftx.u[i,j+1] = shiftx.u[i-1,j+1] = 0.0
    @test iszero(shiftx.u)
    @test iszero(shiftx.v)
    shifty = Edges(Dual,cellzero)
    interpolate!(shifty,faceyunit)
    @test shifty.v[i,j] == shifty.v[i,j-1] == shifty.v[i+1,j] == shifty.v[i+1,j-1] == 0.25
    shifty.v[i,j] = shifty.v[i,j-1] = shifty.v[i+1,j] = shifty.v[i+1,j-1] = 0.0
    @test iszero(shifty.u)
    @test iszero(shifty.v)
  end

  @testset "Dual cell center data shift to dual face" begin
    w = Edges(Dual,cellzero)
    interpolate!(w,cellunit)
    @test w.u[i,j] == w.u[i-1,j] == 0.5
    w.u[i,j] = w.u[i-1,j] = 0.0
    @test iszero(w.u)
    @test w.v[i,j] == w.v[i,j-1] == 0.5
    w.v[i,j] = w.v[i,j-1] = 0.0
    @test iszero(w.v)
  end

  @testset "Face data shift to dual cell center" begin
    cellx = Nodes(Dual,cellzero)
    celly = Nodes(Dual,cellzero)
    interpolate!((cellx,celly),facexunit)
    @test cellx[i,j] == 0.5 && cellx[i,j+1] == 0.5
    cellx[i,j] = cellx[i,j+1] = 0.0
    @test iszero(cellx)
    @test iszero(celly)

    cellx = Nodes(Dual,cellzero)
    celly = Nodes(Dual,cellzero)
    interpolate!((cellx,celly),faceyunit)
    @test celly[i,j] == 0.5 && celly[i+1,j] == 0.5
    celly[i,j] = celly[i+1,j] = 0.0
    @test iszero(cellx)
    @test iszero(celly)
  end

  @testset "div curl" begin
    @test iszero(divergence(curl(cellunit)))
    @test iszero(curl(grad(nodeunit)))
  end

  L = plan_laplacian(nx,ny;with_inverse=true)

  @testset "Laplacian of the LGF" begin
    ψ = L\cellunit
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())
  end

  @testset "LGF for Helmholtz equation" begin
    alpha = 0.02
    LH(i,j,f::Function,α) = im*α*f(i,j,α)-(f(i-1,j,α)+f(i+1,j,α)+f(i,j+1,α)+f(i,j-1,α)-4*f(i,j,α))

    i0, j0 = rand(1:100), rand(0:100)
    @test abs(LH(i0,j0,Fields.lgf_helmholtz,alpha)) < 100.0*eps()

    @test isapprox(real(LH(0,0,Fields.lgf_helmholtz,alpha)),1.0;atol=100.0*eps())

  end

end

##### COMPLEX GRID #####

@testset "Complex Grid Routines" begin

# size
nx = 12; ny = 12

# sample point
i = 5; j = 7

cellzero = Nodes(Dual,(nx,ny),dtype=ComplexF64)
nodezero = Nodes(Primal,cellzero)
facezero = Edges(Primal,cellzero)
dualfacezero = Edges(Dual,cellzero)

a = 1.0+2.0im

cellunit = deepcopy(cellzero)
cellunit[i,j] = a

nodeunit = deepcopy(nodezero)
nodeunit[i,j] = a

facexunit = deepcopy(facezero)
facexunit.u[i,j] = a

faceyunit = deepcopy(facezero)
faceyunit.v[i,j] = a

dualfacexunit = deepcopy(dualfacezero)
dualfacexunit.u[i,j] = a

dualfaceyunit = deepcopy(dualfacezero)
dualfaceyunit.v[i,j] = a

@testset "Basic array operations" begin
  w = zero(cellunit)
  w .= cellunit
  @test w[i,j] == a
  q = similar(facexunit)
  q .= facexunit
  @test q.u[i,j] == a
  @test iszero(q.v)
end

@testset "Basic complex operations" begin
  w = conj(cellunit)
  @test typeof(w) == typeof(cellunit)
  @test w[i,j] == conj(a)

  w = real(cellunit)
  @test typeof(w) <: Nodes{Dual,nx,ny,Float64}
  @test real(w[i,j]) == real(a) && imag(w[i,j]) == 0.0

  w = imag(cellunit)
  @test typeof(w) <: Nodes{Dual,nx,ny,Float64}
  @test real(w[i,j]) == imag(a) && imag(w[i,j]) == 0.0

end

@testset "Inner products and norms" begin
  w = zero(cellunit)
  i0, j0 = rand(2:nx-1), rand(2:ny-1)
  w[i0,j0] = 1.0im
  @test Fields.norm(w)*sqrt((nx-2)*(ny-2)) == 1.0
  w .= 1.0im
  @test Fields.norm(w) == 1.0

  @test 2*w == w*2

  p = Nodes(Primal,w)
  p .= 1.0im
  @test Fields.norm(p) == 1.0
  p2 = deepcopy(p)
  @test Fields.dot(p,p2) == 1.0
  @test Fields.norm(p-p2) == 0.0

  q = Edges(Dual,w)
  q.u .= 1.0im
  q2 = deepcopy(q)
  @test Fields.dot(q,q2) == 1.0

  q = Edges(Primal,w)
  q.u .= 1.0im
  q2 = deepcopy(q)
  @test Fields.dot(q,q2) == 1.0

  @test Fields.integrate(w) == 1.0im

  @test Fields.integrate(p) == 1.0im

  q .= 1.0im
  @test Fields.norm(2*q) == sqrt(8)

end

@testset "Dual cell center data Laplacian" begin
  lapcell = laplacian(cellunit)
  @test lapcell[i,j] == -4.0*a
  lapcell[i,j] = 0.0
  @test lapcell[i+1,j] == lapcell[i-1,j] == lapcell[i,j-1] == lapcell[i,j+1] == a
  lapcell[i+1,j] = lapcell[i-1,j] = lapcell[i,j-1] = lapcell[i,j+1] = 0.0
  @test iszero(lapcell)
end

@testset "Dual cell center data curl" begin
  q = curl(cellunit)
  @test q.u[i,j-1] == 1.0*a && q.u[i,j] == -1.0*a
  q.u[i,j-1] = q.u[i,j] = 0.0
  @test iszero(q.u)
  @test q.v[i-1,j] == -1.0*a && q.v[i,j] == 1.0*a
  q.v[i-1,j] = q.v[i,j] = 0.0
  @test iszero(q.v)
end

@testset "Dual cell node gradient" begin
  q = grad(nodeunit)
  @test q.u[i,j] == 1.0*a && q.u[i+1,j] == -1.0*a
  q.u[i,j] = q.u[i+1,j] = 0.0
  @test iszero(q.u)
  @test q.v[i,j] == 1.0*a && q.v[i,j+1] == -1.0*a
  q.v[i,j] = q.v[i,j+1] = 0.0
  @test iszero(q.v)
end

@testset "Face data curl" begin
  cellcurl = curl(facexunit)
  @test cellcurl[i,j] == -1.0*a && cellcurl[i,j+1] == 1.0*a
  cellcurl[i,j] = cellcurl[i,j+1] = 0.0
  @test iszero(cellcurl)
  cellcurl = curl(faceyunit)
  @test cellcurl[i,j] == 1.0*a && cellcurl[i+1,j] == -1.0*a
  cellcurl[i,j] = cellcurl[i+1,j] = 0.0
  @test iszero(cellcurl)
end

@testset "Face data divergence" begin
  nodediv = divergence(facexunit)
  @test nodediv[i,j] == -1.0*a && nodediv[i-1,j] == 1.0*a
  nodediv[i,j] = nodediv[i-1,j] = 0.0
  @test iszero(nodediv)
  nodediv = divergence(faceyunit)
  @test nodediv[i,j] == -1.0*a && nodediv[i,j-1] == 1.0*a
  nodediv[i,j] = nodediv[i,j-1] = 0.0
  @test iszero(nodediv)
end

@testset "Face data Laplacian" begin
  lap = laplacian(facexunit)
  @test lap.u[i,j] == -4.0*a
  lap.u[i,j] = 0.0
  @test lap.u[i+1,j] == lap.u[i-1,j] == lap.u[i,j-1] == lap.u[i,j+1] == 1.0*a
  lap.u[i+1,j] = lap.u[i-1,j] = lap.u[i,j-1] = lap.u[i,j+1] = 0.0
  @test iszero(lap.u)
  @test iszero(lap.v)

  lap = laplacian(faceyunit)
  @test lap.v[i,j] == -4.0*a
  lap.v[i,j] = 0.0
  @test lap.v[i+1,j] == lap.v[i-1,j] == lap.v[i,j-1] == lap.v[i,j+1] == 1.0*a
  lap.v[i+1,j] = lap.v[i-1,j] = lap.v[i,j-1] = lap.v[i,j+1] = 0.0
  @test iszero(lap.u)
  @test iszero(lap.v)
end

@testset "Dual face data divergence" begin
  celldiv = divergence(dualfacexunit)
  @test celldiv[i,j] == 1.0*a && celldiv[i+1,j] == -1.0*a
  celldiv[i,j] = celldiv[i+1,j] = 0.0
  @test iszero(celldiv)
  celldiv = divergence(dualfaceyunit)
  @test celldiv[i,j] == 1.0*a && celldiv[i,j+1] == -1.0*a
  celldiv[i,j] = celldiv[i,j+1] = 0.0
  @test iszero(celldiv)
end

@testset "Face data shift to dual face" begin
  shiftx = Edges(Dual,cellzero)
  interpolate!(shiftx,facexunit)
  @test shiftx.u[i,j] == shiftx.u[i-1,j] == shiftx.u[i,j+1] == shiftx.u[i-1,j+1] == 0.25*a
  shiftx.u[i,j] = shiftx.u[i-1,j] = shiftx.u[i,j+1] = shiftx.u[i-1,j+1] = 0.0
  @test iszero(shiftx.u)
  @test iszero(shiftx.v)
  shifty = Edges(Dual,cellzero)
  interpolate!(shifty,faceyunit)
  @test shifty.v[i,j] == shifty.v[i,j-1] == shifty.v[i+1,j] == shifty.v[i+1,j-1] == 0.25*a
  shifty.v[i,j] = shifty.v[i,j-1] = shifty.v[i+1,j] = shifty.v[i+1,j-1] = 0.0
  @test iszero(shifty.u)
  @test iszero(shifty.v)
end

@testset "Dual cell center data shift to dual face" begin
  w = Edges(Dual,cellzero)
  interpolate!(w,cellunit)
  @test w.u[i,j] == w.u[i-1,j] == 0.5*a
  w.u[i,j] = w.u[i-1,j] = 0.0
  @test iszero(w.u)
  @test w.v[i,j] == w.v[i,j-1] == 0.5*a
  w.v[i,j] = w.v[i,j-1] = 0.0
  @test iszero(w.v)
end

@testset "Face data shift to dual cell center" begin
  cellx = Nodes(Dual,cellzero)
  celly = Nodes(Dual,cellzero)
  interpolate!((cellx,celly),facexunit)
  @test cellx[i,j] == 0.5*a && cellx[i,j+1] == 0.5*a
  cellx[i,j] = cellx[i,j+1] = 0.0
  @test iszero(cellx)
  @test iszero(celly)

  cellx = Nodes(Dual,cellzero)
  celly = Nodes(Dual,cellzero)
  interpolate!((cellx,celly),faceyunit)
  @test celly[i,j] == 0.5*a && celly[i+1,j] == 0.5*a
  celly[i,j] = celly[i+1,j] = 0.0
  @test iszero(cellx)
  @test iszero(celly)
end

@testset "div curl" begin
  @test iszero(divergence(curl(cellunit)))
  @test iszero(curl(grad(nodeunit)))
end

L = plan_laplacian(nx,ny;with_inverse=true,dtype=ComplexF64)

@testset "Laplacian of the LGF" begin
  ψ = L\cellunit
  lapψ = L*ψ
  @test lapψ[i,j]≈1.0*a
  @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=100.0*eps()) &&
          isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=100.0*eps())
end

α = 0.07
LH = plan_helmholtz(nx,ny,α;with_inverse=true)

@testset "Helmholtz of the LGF" begin
  ψ = LH\cellunit
  helmψ = LH*ψ
  @test helmψ[i,j]≈1.0*a
  @test isapprox(maximum(abs.(helmψ[Not(i),:])),0.0;atol=100.0*eps()) &&
          isapprox(maximum(abs.(helmψ[:,Not(j)])),0.0;atol=100.0*eps())
end

end

@testset "Fields" begin
    @testset "Hadamard Product" begin
        edges_p  = Edges(Primal,(30,40))
        edges_p.u .= rand(Float64,size(edges_p.u))

        # Should be safe for the output to be the same as the input
        edges_p2 = edges_p ∘ edges_p
        product!(edges_p, edges_p, edges_p)
        @test edges_p2.u == edges_p.u
        @test edges_p2.v == edges_p.v

        edges_d  = Edges{Dual, 30, 40, Float64}()
        @test_throws MethodError (edges_p ∘ edges_d)
    end

    @testset "Discrete Laplacian" begin
        s = Nodes(Dual,(30,40))
        s[3:end-2, 3:end-2] .= rand(26, 36)

        L = plan_laplacian(30, 40)

        @test L*s ≈ -curl(curl(s))

        @test_throws MethodError (L \ s)

        L = plan_laplacian(30, 40, with_inverse = true, fftw_flags = FFTW.PATIENT)
        @test L \ (L*s) ≈ s

        L! = plan_laplacian!(30, 40, with_inverse = true, fftw_flags = FFTW.PATIENT)
        sold = deepcopy(s)
        L! \ (L!*s)
        @test s ≈ sold
    end

    @testset "Integrating factor" begin
        s = Nodes(Dual,(30,40))
        s[15,15] = 1.0

        E1 = plan_intfact(1,s)
        E2 = plan_intfact(2,s)

        @test E1*(E1*s) ≈ E2*s

        E! = plan_intfact!(2,s)
        s2 = deepcopy(s)
        E! * s

        @test s ≈ E2*s2

    end

    @testset "Discrete Divergence" begin
        s = Nodes{Dual, 5, 4, Float64}()
        s .= rand(5, 4)

        @test iszero(divergence(curl(s)))

        s = Nodes{Primal, 5, 4, Float64}()
        q′ = Edges{Primal, 5, 4, Float64}()
        q′.u .= reshape(1:15, 5, 3)
        q′.v .= reshape(1:16, 4, 4)

        # Not sure if this is the behavior we want yet
        # Currently, the ghost cells are not affected
        # by the divergence operator
        #s .= 1.0
        divergence!(s, q′)
        @test s == [ 5.0  5.0  5.0
                     5.0  5.0  5.0
                     5.0  5.0  5.0
                     5.0  5.0  5.0 ]
    end

    @testset "Discrete Curl" begin
        s = Nodes{Dual, 5, 4, Float64}()
        s .= reshape(1:20, 4, 5)'

        q = curl(s)

        @test q.u == [ 1.0  1.0  1.0
                       1.0  1.0  1.0
                       1.0  1.0  1.0
                       1.0  1.0  1.0
                       1.0  1.0  1.0 ]

        @test q.v == [ -4.0  -4.0  -4.0  -4.0
                       -4.0  -4.0  -4.0  -4.0
                       -4.0  -4.0  -4.0  -4.0
                       -4.0  -4.0  -4.0  -4.0 ]
    end

    @testset "Shifting Primal Edges to Dual Edges" begin

        q = Edges{Primal, 5, 4, Float64}()
        q.u .= reshape(1:15, 5, 3)
        q.v .= reshape(1:16, 4, 4)
        Qq = Edges{Dual, 5, 4, Float64}()

        interpolate!(Qq,q)
        @test Qq.u == [ 0.0  4.0  9.0  0.0
                        0.0  5.0  10.0 0.0
                        0.0  6.0  11.0 0.0
                        0.0  7.0  12.0 0.0 ]

        @test Qq.v == [ 0.0  0.0   0.0
                        3.5  7.5  11.5
                        4.5  8.5  12.5
                        5.5  9.5  13.5
                        0.0  0.0   0.0 ]

    end

    @testset "Shifting Dual Edges to Primal Edges" begin

        q = Edges{Dual, 5, 4, Float64}()
        q.u .= reshape(1:16, 4, 4)
        q.v .= reshape(1:15, 5, 3)
        v = Edges{Primal, 5, 4, Float64}()
        interpolate!(v,q)

        @test v.u == [ 0.0  0.0   0.0
                       3.5  7.5  11.5
                       4.5  8.5  12.5
                       5.5  9.5  13.5
                       0.0  0.0   0.0 ]

        @test v.v == [ 0.0  4.0   9.0  0.0
                        0.0  5.0  10.0  0.0
                        0.0  6.0  11.0  0.0
                        0.0  7.0  12.0  0.0 ]

    end

    @testset "Shifting Dual Nodes to Dual Edges" begin

        w = Nodes{Dual, 5, 4, Float64}()
        w .= reshape(1:20, 5, 4)

        Ww = Edges{Dual, 5, 4, Float64}()
        interpolate!(Ww,w)

        @test Ww.u == [ 0.0  6.5  11.5  0.0
                        0.0  7.5  12.5  0.0
                        0.0  8.5  13.5  0.0
                        0.0  9.5  14.5  0.0 ]
        @test Ww.v == [ 0.0   0.0   0.0
                        4.5   9.5  14.5
                        5.5  10.5  15.5
                        6.5  11.5  16.5
                        0.0   0.0   0.0 ]
    end

    @testset "Physical grid" begin

        g = PhysicalGrid((-1.0,3.0),(-2.0,3.0),0.02)
        @test size(g) == (208,252)
        @test size(g,1) == 208
        @test size(g,2) == 252
        @test length(g) == 208*252
        @test origin(g) == (54,101)
        @test cellsize(g) == 0.02
        @test limits(g,1) == (-1.06,3.06)
        @test limits(g,2) == (-2.0,3.0)

    end
end


text:
abstract type Edge{N,T} end
Base.broadcastable(e::Edge) = Ref(e)
Base.@propagate_inbounds function Base.getindex(e::Edge, i::Int)
    getfield(e, :points)[i]
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
function particle_in_box(b, n::Int)
    L = get_length(b)
    N = sqrt(2/L)
    a = n * π / L
    #return N .* sin.(a .* (b .+ 0.5*L) )
    return N .* sin.(a .* (b .- b[begin]) )
end

function initial_orbitals(b::AbstractBasis)
    l = length(b)
    ψ = zeros(l,l)
    for j in 1:l-2 # end points are zero so -2
        ψ[:,j] = particle_in_box(b, j)
    end
    ψ[end,:] .= 0  # Improve accuracy at the end
    # Add nonzero end functions
    ψ[1,end-1] = 1
    ψ[end,end] = 2
    w = get_weight(b)
    return gram_schmit(ψ, w)
end

function gram_schmit(orbitals::AbstractMatrix, w::AbstractVector)
    out = similar(orbitals)
    norm = sqrt( sum( w .* orbitals[:,1].^2 ) )
    out[:,1] = orbitals[:,1] ./ norm
    for i in 2:size(out,2)
        tmp = orbitals[:,i]
        for j in 1:i-1
           tmp -= sum( out[:,j] .* w .* orbitals[:,i] ) .* out[:,j]
        end
        norm = sqrt( sum( w .* tmp.^2 ) )
        out[:,i] = tmp ./ norm
    end
    return out
end

"""
    solve_hartree_fock(b::AbstractBasis, Vn; Ve=x->exp(-x^2), rtol=1E-9, max_iter=100)

Solve Hartree Fock equations for two electrons for given nuclear potential `Vn`.

Electron-electron repulsion is fixed to exp(-Δx^2).

# Returns
- `E`         :  Energy of electronic system
- `Eorb`      :  Orbital energy
- `orbitals`  :  Orbitals in collumns
"""
function solve_hartree_fock(b::AbstractBasis, Vn ; Ve=x->exp(-x^2), rtol=1E-9, max_iter=100)
    # Solve FC = SCε for two electrons
    orbitals = initial_orbitals(b)
    F = fock_matrix(b, orbitals; Vn=Vn)
    h₁ = one_electron_operator(b, Vn)

    # We have orthogonal, but not orthonormal basis
    X = Diagonal( 1.0 ./ sqrt.( get_weight(b) ) )

    # Energy is h₁ + F (=2h₁+J) for two electrons
    E₀ = orbitals[:,1]' * ( h₁ + F ) * orbitals[:,1]

    @info "Energy at start = $(E₀)"
    eo, orbitals = eigen( X * F * X )
    orbitals = X * orbitals
    F = fock_matrix(b, orbitals; Vn=Vn)
    for i in 1:max_iter
        E₁ = orbitals[:,1]' * ( h₁ + F ) * orbitals[:,1]
        @info "Energy at iteration $i = $(E₁)"
        if abs(E₁-E₀) < rtol
            break
        else
            E₀ = E₁
            eo, orbitals = eigen( X * F * X )
            orbitals = X * orbitals
            F = fock_matrix(b, orbitals; Vn=Vn)
        end
    end
    E = orbitals[:,1]' * ( h₁ + F ) * orbitals[:,1]
    @info "Final energy = $(E)"
    return E, eo, orbitals
end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation