text:
# Code was originally written in the Jupyter Notebook, so it is recommended to run this script from there instead.
# Some code may not work in console since displaying output images is not supported in console.

using Luxor

function generate_tree(treex, treey, size, hue, snow = false)
#     Set color and opacity
    setopacity(1)
    sethue(hue)
#     Decide number of triangles to draw for tree depending on size
    for i in 1:round(size/4, digits = 0)+ 3
#         Decide position of each triangle
        pointa = Point(treex, treey + (size+1)*i)
#         Draw triangle
        ngonside(pointa, size*i, 3, pi/6, :fill)
    end
#     If tree is snow-covered, then render the two top triangles as white.
#     This is placed after the loop since other triangles draw over the first 2 in the loop, and make it look weird.
    if (snow)
        setopacity(1)
        sethue("snow")
        ngonside(Point(treex, treey + (size+1)*1), size, 3, pi/6, :fill)
        ngonside(Point(treex, treey + (size+1)*2), size*2, 3, pi/6, :fill)
    end
end

function generate_house(x, y)

#     Setup
    setopacity(1)

#     Draw polygon on the left (left side of house, below the snowy roof)
    sethue(0.7451, 0.5882, 0.549)
    poly([Point(x-157, y+60), Point(x-86, y+50), Point(x-87, y+99), Point(x-157, y+103)], close = true, :fill)
    sethue(0.6588, 0.4745, 0.4431)
    poly([Point(x-157, y+60), Point(x-86, y+50), Point(x-87, y+99), Point(x-157, y+103)], close = true, :stroke)

#     Draw the front of the house (Pentagon shaped face)
    sethue(0.7451, 0.5882, 0.549)
    poly([Point(x-86, y+50), Point(x-87, y+99), Point(x+44, y+103), Point(x+42, y+53), Point(x-11, y-41)], close = true, :fill)
    sethue(0.6588, 0.4745, 0.4431)
    poly([Point(x-86, y+50), Point(x-87, y+99), Point(x+44, y+103), Point(x+42, y+53), Point(x-11, y-41)], close = true, :stroke)

#     Draw the snowy roof
    sethue("snow")
    poly([Point(x-157, y+60), Point(x-86, y+50), Point(x-11, y-41), Point(x-116, y+11)], close = true, :fill)
    sethue("snow2")
    poly([Point(x-157, y+60), Point(x-86, y+50), Point(x-11, y-41), Point(x-116, y+11)], close = true, :stroke)

#     Draw the chimney
    sethue(0.7451, 0.5882, 0.549)
    poly([Point(x-114, y+20), Point(x-90, y+20), Point(x-89, y-45), Point(x-106, y-40)], close = true, :fill)
    sethue(0.6588, 0.4745, 0.4431)
    poly([Point(x-114, y+20), Point(x-90, y+20), Point(x-89, y-45), Point(x-106, y-40)], close = true, :stroke)

#     Generate the bezier for the smoke and draw it
    sethue("whitesmoke")
    setline(20)
    setopacity(0.2)
    drawbezierpath(makebezierpath([Point(x-100, y-45), Point(x-80, y-50), Point(x-130, y-80), Point(x-30, y-132), Point(x-100, y-200)]), close = false, :stroke)

#     Reset to original stroke properties after drawing smoke
    setopacity(1)
    setline(1)

#     Draw Windows
#     Draw left most window
    sethue("gold")
    poly([Point(x-150, y+87), Point(x-130, y+85), Point(x-131, y+61), Point(x-150, y+64)], close = true, :fill)
    sethue("tan4")
    poly([Point(x-150, y+87), Point(x-130, y+85), Point(x-131, y+61), Point(x-150, y+64)], close = true, :stroke)
    line(Point(x-140, y+85), Point(x-141, y+63), :stroke)
    line(Point(x-130.5, y+73), Point(x-150, y+76), :stroke)

#     Draw the window to the immediate right of the leftmost window
    sethue("gold")
    poly([Point(x-117, y+82), Point(x-97, y+80), Point(x-97, y+56), Point(x-117, y+59)], close = true, :fill)
    sethue("tan4")
    poly([Point(x-117, y+82), Point(x-97, y+80), Point(x-97, y+56), Point(x-117, y+59)], close = true, :stroke)
    line(Point(x-107, y+80), Point(x-108, y+57), :stroke)
    line(Point(x-117, y+71), Point(x-97, y+69), :stroke)

#     Draw the window next to the door
    sethue("gold")
    poly([Point(x, y+90), Point(x+30, y+90), Point(x+30, y+55), Point(x, y+53)], close = true, :fill)
    sethue("tan4")
    poly([Point(x, y+90), Point(x+30, y+90), Point(x+30, y+55), Point(x, y+53)], close = true, :stroke)
    line(Point(x, y+72), Point(x+30, y+73), :stroke)
    line(Point(x+15, y+54), Point(x+15, y+90), :stroke)

#     Draw the circular window above the door
    sethue("gold")
    circle(x-15, y+10, 12, :fill)
    sethue("tan4")
    circle(x-15, y+10, 12, :stroke)
    line(Point(x-27, y+10), Point(x-3, y+10), :stroke)
    line(Point(x-15, y+22), Point(x-15, y-2), :stroke)

#     Draw the door
    setopacity(0.3)
    sethue("tan3")
    poly([Point(x-61, y+100), Point(x-30, y+100.8),  Point(x-29, y+48), Point(x-60, y+48)], close = true, :fill)
    setopacity(0.4)
    sethue("tan4")
    poly([Point(x-61, y+100), Point(x-30, y+100.8),  Point(x-29, y+48), Point(x-60, y+48)], close = true, :stroke)
end

# Setup drawing
Drawing(600, 314, "Scene.png")
# Set current point to origin
origin()

# Set up the sky
nightsky = blend(Point(0, -300), Point(0, 230), "midnightblue", "steelblue2")
setblend(nightsky)
box(O, 600, 314, :fill)

# Draw the stars
sethue("grey100")
circle(-136, -78, 1, :fill)
circle(-26, -88, 1, :fill)
circle(120, -117, 1, :fill)

# Draw moon
sethue("white")
setopacity(0.8)
circle(200, -100, 20, :fill)

# Generate size noise for trees
size = noise.(range(0, 5, length=10), detail = 10, persistence = 0)

# Generate new noise for ground on which house lies
yfore = noise.(range(4, 8, length=10), detail = 10, persistence = 0)

# Process ground noise to an array of points, and get array of lowest y-coordinate of noise (we draw the rect using this)
y_fore = .-(20, 70*yfore)
# Generate the Points, and put them in an array (WARNING: LOTS of Arbitrary constants)
gpointa = Point.(-370 .+ 70*collect(1:10), y_fore[collect(1:10)] .+ (yfore[collect(1:10)]*4 .+ 12)*8.2)
# Get array of y-coordinates of noise
abcd = map(x -> x.y , gpointa)

# Process tree size noise to an array of points, and get array of lowest y-coordinate of noise (we draw the ground using this)
y_pos = .-(-10, 50*size)
# Generate the Points, and put them in an array (WARNING: LOTS of arbitrary constants)
gpoint = Point.(-330 .+ 60*collect(1:10), y_pos[collect(1:10)] .+ (size[collect(1:10)]*4 .+ 10)*8.2)
# Get array of y-coordinates of noise
abc = map(x -> x.y , gpoint)

# Add drawing endpoints to make Bezier cover a large chunk of ground without needing the rect (Background terrain)
push!(gpoint, Point(400, gpoint[end].y), Point(301, 158), Point(-301, 158), Point(-350, gpoint[1].y))
# Add drawing endpoints to make Bezier cover a large chunk of ground without needing the rect (House terrain)
push!(gpointa, Point(400, gpointa[end].y), Point(301, 158), Point(-301, 158), Point(-350, gpointa[1].y))

# Draw trees starting from edge of the screen to the other edge, while using size and position noise.
for treenum in 1:10
    x = -330 + 60*treenum
    generate_tree(x, -10 - 50*size[treenum], size[treenum]*4 + 10, "dodgerblue4")
end

# Create rect to cover up rest of terrain that Bezier didn't cover by using
# lowest point of bezier (minimum of noise y-coordinate), But using maximum since values are negative. (Background)
sethue("lavenderblush2")
box(Point(-301, maximum(abc)), Point(301, 158), :fill, vertices = false)
poly(polyfit(gpoint), :fill)

# Create rect to cover up rest of terrain that Bezier didn't cover by using
# lowest point of bezier (minimum of noise y-coordinate), But using maximum since values are negative. (House Ground)
sethue("lavenderblush1")
setopacity(1)
box(Point(-301, maximum(abcd)), Point(301, 158), :fill, vertices = false)
poly(polyfit(gpointa), :fill)

# Draw the house with slightly adjusted coordinates
generate_house(50, -10)

# Draw snow-covered trees next to house without using noise
generate_tree(250, 20, 11, "forestgreen", true)
generate_tree(210, 70, 9, "forestgreen", true)
generate_tree(150, 30, 12, "forestgreen", true)

# Draw terrain closest to viewer using a large triangle partially offscreen
sethue("snow")
setopacity(0.9)
drawbezierpath(makebezierpath([Point(-400, 100), Point(200, 200), Point(-200, 350)]), :fill)

# Draw our final tree in foreground
generate_tree(-250, -40, 16, "forestgreen", true)

finish()
preview()

max_stars_repo_name: soumitradev/LuxorScene

max_issues_repo_name: soumitradev/LuxorScene

max_forks_repo_name: soumitradev/LuxorScene


text:
function showarray(arr)
    println(summary(arr))
    Base.print_matrix(IOContext(stdout, :limit => true), arr)
    println()
end

function strided(carr::Array{Float64,1}, scale::Int64, offset=-1)
    offset = offset==-1 ? scale : offset
    stridedarr = [carr[n:n+scale-1] for n=1:offset:(length(carr)-scale+1)]
end

function fbmArea(z)
    nx, ny = size(z)
    area = 0.0
    for i=1:nx-1
        for j=1:ny-1
            v1 = [0; 1; z[i,j+1] - z[i,j]]
            v2 = [1; 0; z[i+1,j] - z[i,j]]
            v3 = [0; -1; z[i+1,j] - z[i+1,j+1]]
            v4 = [-1; 0; z[i,j+1] - z[i+1,j+1]]
            area += (norm(v2 × v1)+norm(v4 × v3))/2
        end
    end
    return area/(nx*ny)
end

ext: jl


text:
@testset "cube" begin

    x = reshape(1:24, 2, 3, 4)
    @test is_transformable(x)

    @testset "apply" begin

        @testset "OneToOne" begin
            T = FakeOneToOneTransform()
            @test FeatureTransforms.apply(x, T) == ones(2, 3, 4)
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == ones(2)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[2, 3]) == ones(2, 2, 4)
        end

        @testset "OneToMany" begin
            T = FakeOneToManyTransform()
            @test FeatureTransforms.apply(x, T) == ones(2, 6, 4)
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == ones(2, 2)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[2, 3]) == ones(2, 4, 4)
        end

        @testset "ManyToOne" begin
            T = FakeManyToOneTransform()
            @test FeatureTransforms.apply(x, T; dims=1) == ones(3, 4)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[1, 2]) == ones(2, 4)
            @test FeatureTransforms.apply(x, T; dims=3, inds=[1, 3]) == ones(2, 3)
        end

        @testset "ManyToMany" begin
            T = FakeManyToManyTransform()
            @test FeatureTransforms.apply(x, T) == ones(2, 6, 4)
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == ones(2, 2)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[2, 3]) == ones(2, 4, 4)
        end

    end

    @testset "apply!" begin

        T = FakeOneToOneTransform()

        _x = copy(x)
        FeatureTransforms.apply!(_x, T)
        @test _x == ones(2, 3, 4)

        # https://github.com/invenia/FeatureTransforms.jl/issues/68
        _x = copy(x)
        @test_broken FeatureTransforms.apply!(_x, T; dims=2, inds=[1, 2])
        @test_broken _x == cat(
            [1 1 5; 1 1 6], [1 1 11; 1 1 12], [1 1 17; 1 1 18], [1 1 23; 1 1 24]; dims=3
        )
    end

    @testset "apply_append" begin

        @testset "OneToOne" begin

            T = FakeOneToOneTransform()

            @test FeatureTransforms.apply_append(x, T; append_dim=1) == vcat(x, ones(2, 3, 4))
            @test FeatureTransforms.apply_append(x, T; append_dim=2) == hcat(x, ones(2, 3, 4))
            @test FeatureTransforms.apply_append(x, T; append_dim=3) == cat(x, ones(2, 3, 4), dims=3)

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=1, inds=[1], append_dim=1),
                vcat(x, ones(1, 3, 4))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=2, inds=[1, 2], append_dim=2),
                hcat(x, ones(2, 2, 4))
            )
        end

        @testset "OneToMany" begin

            T = FakeOneToManyTransform()

            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=3)

            @test isequal(
                FeatureTransforms.apply_append(x, T; append_dim=2),
                hcat(x, ones(2, 6, 4))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=2, inds=[1, 2], append_dim=2),
                hcat(x, ones(2, 4, 4))
            )
        end

        @testset "ManyToOne" begin

            T = FakeManyToOneTransform()

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=1, append_dim=1),
                vcat(x, ones(1, 3, 4))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=2, append_dim=2),
                hcat(x, ones(2, 1, 4))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=3, append_dim=3),
                cat(x, ones(2, 3, 1); dims=3)
            )
        end

        @testset "ManyToMany" begin

            T = FakeManyToManyTransform()

            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=3)

            @test isequal(
                FeatureTransforms.apply_append(x, T; append_dim=2),
                hcat(x, ones(2, 6, 4))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=2, inds=[1, 2], append_dim=2),
                hcat(x, ones(2, 4, 4))
            )
        end

    end

end


text:
#=

JuliaBolt is ported from Neobolt: Neo4j Bolt connector for Python
https://github.com/neo4j-drivers/neobolt

Original notices:

# Copyright (c) 2002-2018 "Neo4j,"
# Neo4j Sweden AB [http://neo4j.com]
#
# This file is part of Neo4j.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

=#

#####################################
#             Unpacker              #
#####################################

mutable struct Unpacker
    source::Union{MessageFrame, Nothing}

    Unpacker() = new(nothing)
end

# Unpacker ------------------------ #

read(u::Unpacker, n::Integer=1) = read(u.source, n)

# Unpacker ------------------------ #

read_int(u::Unpacker) = read_int(u.source)

# Unpacker ------------------------ #

function unpack(u::Unpacker)
    marker = read_int(u)

    if marker == -1
        throw(ErrorException("Nothing to unpack"))
    end

    # Tiny Integer
    if 0x00 <= marker <= 0x7F
        return marker
    elseif 0xF0 <= marker <= 0xFF
        return Int8(marker - 256)

    # Null
    elseif marker == 0xC0
        return nothing

    # Float
    elseif marker == 0xC1
        return struct_unpack(Float64, read(u, 8))

    # Boolean
    elseif marker == 0xC2
        return false
    elseif marker == 0xC3
        return true

    # Integer
    elseif marker == 0xC8
        return struct_unpack(Int8, read(u, 1))
    elseif marker == 0xC9
        return struct_unpack(Int16, read(u, 2))
    elseif marker == 0xCA
        return struct_unpack(Int32, read(u, 4))
    elseif marker == 0xCB
        return struct_unpack(Int64, read(u, 8))

    # Bytes
    elseif marker == 0xCC
        size_uint8 = struct_unpack(UInt8, read(u, 1))
        return read(u, size_uint8)
    elseif marker == 0xCD
        size_uint16 = struct_unpack(UInt16, read(u, 2))
        return read(u, size_uint16)
    elseif marker == 0xCE
        size_uint32 = struct_unpack(UInt32, read(u, 4))
        return read(u, size_uint32)

    else
        marker_high = marker & 0xF0
        # String
        if marker_high == 0x80  # TINY_STRING
            return String(read(u, marker & 0x0F))
        elseif marker == 0xD0  # STRING_8
            size_uint8 = struct_unpack(UInt8, read(u, 1))
            return String(read(u, size_uint8))
        elseif marker == 0xD1  # STRING_16
            size_uint16 = struct_unpack(UInt16, read(u, 2))
            return String(read(u, size_uint16))
        elseif marker == 0xD2  # STRING_32
            size_uint32 = struct_unpack(UInt32, read(u, 4))
            return String(read(u, size_uint32))

        # List
        elseif 0x90 <= marker <= 0x9F || 0xD4 <= marker <= 0xD7
            return unpack_list(u, marker)

        # Map
        elseif 0xA0 <= marker <= 0xAF || 0xD8 <= marker <= 0xDB
            return unpack_map(u, marker)

        # Structure
        elseif 0xB0 <= marker <= 0xBF || 0xDC <= marker <= 0xDD
            (size, tag) = unpack_structure_header(u, marker)
            value = Structure(tag, Array{Any}(nothing, size))
            for i in 1:length(value.fields)
                value.fields[i] = unpack(u)
            end
            return value

        elseif marker == 0xDF  # END_OF_STREAM
            return :EndOfStream

        else
            throw(ErrorException("Unknown PackStream marker $(marker)"))

        end
    end
end

# Unpacker ------------------------ #

function unpack_list(u::Unpacker)
    marker = read_int(u)
    return unpack_list(u, marker)
end

# Unpacker ------------------------ #

function unpack_list(u::Unpacker, marker::UInt8)
    marker_high = marker & 0xF0
    if marker_high ==  0x90
        size_ = marker & 0x0F
        if size_ == 0
            return []
        elseif size_ == 1
            return [unpack(u)]
        else
            return [unpack(u) for _ in 1:size_]
        end
    elseif marker == 0xD4  # LIST_8
        size_uint8 = struct_unpack(UInt8, read(u, 1))
        return [unpack(u) for _ in 1:size_uint8]
    elseif marker == 0xD5  # LIST_16
        size_uint16 = struct_unpack(UInt16, read(u, 2))
        return [unpack(u) for _ in 1:size_uint16]
    elseif marker == 0xD6  # LIST_32
        size_uint32 = struct_unpack(UInt32, read(u, 4))
        return [unpack(u) for _ in 1:size_uint32]
    elseif marker == 0xD7  # LIST_STREAM
        value = []
        item = nothing
        while item != :EndOfStream
            item = unpack(u)
            if item != :EndOfStream
                push!(value, item)
            end
        end
        return value
    else
        return nothing
    end
end

# Unpacker ------------------------ #

function unpack_map(u::Unpacker)
    marker = read_int(u)
    if 0 <= marker <= 255
        return unpack_map(u, marker)
    else
        return nothing
    end
end

# Unpacker ------------------------ #

function unpack_map(u::Unpacker, marker::UInt8)
    marker_high = marker & 0xF0
    if marker_high == 0xA0
        size = marker & 0x0F
        value = Dict()
        for _ in 1:size
            key = unpack(u)
            value[key] = unpack(u)
        end
        return value
    elseif marker == 0xD8  # MAP_8
        size = struct_unpack(UInt8, read(u, 1))
        value = Dict()
        for _ in 1:size
            key = unpack(u)
            value[key] = unpack(u)
        end
        return value
    elseif marker == 0xD9  # MAP_16
        size, = struct_unpack(UInt16, read(u, 2))
        value = Dict()
        for _ in 1:size
            key = unpack(u)
            value[key] = unpack(u)
        end
        return value
    elseif marker == 0xDA  # MAP_32
        size, = struct_unpack(UInt32, read(u, 4))
        value = Dict()
        for _ in 1:size
            key = unpack(u)
            value[key] = unpack(u)
        end
        return value
    elseif marker == 0xDB  # MAP_STREAM
        value = Dict()
        key = nothing
        while key != :EndOfStream
            key = unpack(u)
            if key != :EndOfStream
                value[key] = unpack(u)
            end
        end
        return value
    else
        return nothing
    end
end

# Unpacker ------------------------ #

function unpack_structure_header(u::Unpacker)
    marker = read_int(u)
    if marker == -1
        return nothing, nothing
    else
        return unpack_structure_header(u, marker)
    end
end

# Unpacker ------------------------ #

function unpack_structure_header(u::Unpacker, marker::UInt8)
    marker_high = marker & 0xF0
    if marker_high == 0xB0  # TINY_STRUCT
        signature = read(u, 1)[1]
        return marker & 0x0F, signature
    elseif marker == 0xDC  # STRUCT_8
        size = struct_unpack(UInt8, read(u, 1))
        signature = read(u, 1)[1]
        return size, signature
    elseif marker == 0xDD  # STRUCT_16
        size = struct_unpack(UInt16, read(u, 2))
        signature = read(u, 1)[1]
        return size, signature
    else
        throw(ErrorException("Expected structure, found marker $(marker)"))
    end
end

# Unpacker ------------------------ #
# # # # # # # # # # # # # # # # # # #

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
module ExampleMOP
import ..VectorMathOptInterface
MOI = VectorMathOptInterface

export MOP

Base.@kwdef struct MOP <: MOI.ModelLike
	name :: String = ""

	variables :: Vector{MOI.VariableIndex} = []

	lower_bounds :: Dict{ MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan}, MOI.GreaterThan} = Dict()

	upper_bounds :: Dict{ MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan}, MOI.LessThan} = Dict()

	scalar_affine_lt_constraints :: Dict{MOI.ConstraintIndex, Tuple{MOI.ScalarAffineFunction, MOI.LessThan}} = Dict()
end

function MOI.is_empty( mop :: MOP )
	_is_empty = true
	for fn in [:variables, :lower_bounds, :upper_bounds, :scalar_affine_eq_constraints]
		_is_empty &= isempty( getfield(mop, fn ))
	end
	return _is_empty
end

function MOI.empty!( mop :: MOP )
	for fn in [:variables, :lower_bounds, :upper_bounds, :scalar_affine_eq_constraints]
		isempty!( getfield(mop, fn ))
	end
	return nothing
end

MOI.supports_incremental_interface( :: MOP ) = true

MOI.supports( ::MOP, ::MOI.Name ) = true
MOI.get(mop::MOP, ::MOI.Name ) = mop.name
function MOI.set( mop :: MOP, ::MOI.Name, new_name)
	mop.name = new_name
end

function rand_scalar_affine_func( vars )
	n_vars = length(vars)
	coeff = rand(n_vars)
	terms = [ MOI.ScalarAffineTerm( c, x ) for (c,x) in zip( coeff, vars ) ]
	return MOI.ScalarAffineFunction( terms, rand() )
end

# VARIABLES
next_id( V ) = isempty( V ) ? 1 : maximum( v.value for v in V )
function MOI.add_variable( mop :: MOP )
	new_var_id = next_id( mop.variables )
	new_var = MOI.VariableIndex(new_var_id)
	push!(mop.variables, new_var )
	return new_var
end
function MOI.delete( mop :: MOP, vi :: MOI.VariableIndex )
	vi_pos = findfirst( isequal(vi), mop.variables )
	if !isnothing(vi_pos)
		deleteat!(mop.variables, vi_pos)
	end
	return nothing
end
function MOI.is_valid( mop :: MOP, vi :: MOI.VariableIndex )
	return vi in mop.variables
end

MOI.get( mop :: MOP, :: MOI.NumberOfVariables ) = length( mop.variables )
MOI.get( mop :: MOP, :: MOI.ListOfVariableIndices ) = mop.variables
MOI.get( mop :: MOP, :: MOI.ListOfVariableAttributesSet ) = MOI.AbstractVariableAttribute[]

#%% BOUND CONSTRAINTS
const LB_CI = MOI.ConstraintIndex{ MOI.VariableIndex, MOI.GreaterThan }
const UB_CI = MOI.ConstraintIndex{ MOI.VariableIndex, MOI.LessThan }

for (T, bound_type) in zip( [ LB_CI, UB_CI ], ["lower_bounds", "upper_bounds"])
	@eval begin

		S = $(T).parameters[2]
		err = $T == LB_CI ? MOI.LowerBoundAlreadySet : MOI.UpperBoundAlreadySet
		fn = Symbol($bound_type)

		MOI.supports_constraint( ::MOP, ::Type{<:MOI.VariableIndex}, ::Type{<:S} ) = true

		function MOI.add_constraint( mop :: MOP, vi :: MOI.VariableIndex, b :: S )
			for con in keys( getfield( mop, fn ) )
				if vi.value == con.value
					throw(err)
				end
			end

			if MOI.is_valid( mop, vi )
				next_const = T( vi.value )
				mop.lower_bounds[ next_const ] = b

				return next_const
			else
				throw( MOI.InvalidIndex(vi) )
			end
		end

		function MOI.delete( mop :: MOP, ci :: $T )
			delete!( getfield(mop, fn), ci )
			return nothing
		end

		function MOI.is_valid( mop :: MOP, ci :: $T )
			return ci in keys( getfield(mop, fn) )
		end

	end
end

# Scalar, Affine Linear LT constraints

function _list_of_all_vars( f :: MOI.ScalarAffineFunction )
	return [ term.variable for term in f.terms ]
end

function MOI.supports_constraint( ::MOP,
	::Type{<:MOI.ScalarAffineFunction}, ::Type{<:MOI.LessThan})
	return true
end

function MOI.add_constraint( mop :: MOP, func :: MOI.ScalarAffineFunction, ub ::MOI.LessThan)
	for var_ind in _list_of_all_vars( func )
		if !MOI.is_valid( mop, var_ind )
			throw( MOI.InvalidIndex )
		end
	end

	constr_ind = MOI.ConstraintIndex{MOI.ScalarAffineFunction, MOI.LessThan}( next_id(keys(mop.scalar_affine_lt_constraints) ) )
	mop.scalar_affine_lt_constraints[ constr_ind ] = ub
	return constr_ind
end

function MOI.get( mop :: MOP, ::MOI.NumberOfConstraints )
	return sum( length( getfield(mop, fn) ) for
		fn in [:lower_bounds, :upper_bounds, :scalar_affine_lt_constraints])
end
end#module


text:
mutable struct NaiveForecaster <: Forecaster
    seasonal::Bool
end

function predict(forecaster::NaiveForecaster, train_dict, info_dict)
    forecast_dict = SeriesDict()

    for (id, y) in train_dict
        info = info_dict[id]
        horizon = info.horizon
        freq = forecaster.seasonal ? info.freq : 1

        y_hat = copy(y)

        for i in 1:horizon
            push!(y_hat, y_hat[end - freq + 1])
        end

        forecast_dict[id] = y_hat[end - horizon + 1:end]
    end

    forecast_dict
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7bd669c8-443f-11ec-045d-bd3b1c74db9b
begin
    import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())
    using Omega, Distributions, UnicodePlots
end

# ╔═╡ a5d02f2a-3310-4f7c-a20f-af43332e2306
md"# Bell's Inequality"

# ╔═╡ 88dda8ef-d7c5-44d4-8e6b-13e20dcd2f03
md"### Quantum Model
First let's set up a quantum model"

# ╔═╡ 112d8745-d12e-4e0e-8af5-d069014db356
md"### Hidden Variable Model
In the hidden variable model, both particles contain hidden state that determines their spin when measures."

# ╔═╡ 03441d7c-3421-42eb-bdc4-474d2ab1c70d

# ╔═╡ Cell order:
# ╠═7bd669c8-443f-11ec-045d-bd3b1c74db9b
# ╠═a5d02f2a-3310-4f7c-a20f-af43332e2306
# ╠═88dda8ef-d7c5-44d4-8e6b-13e20dcd2f03
# ╠═112d8745-d12e-4e0e-8af5-d069014db356
# ╠═03441d7c-3421-42eb-bdc4-474d2ab1c70d

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
export iauPnm00a
"""
Form the matrix of precession-nutation for a given date (including
frame bias), equinox-based, IAU 2000A model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double     TT as a 2-part Julian Date (Note 1)

Returned:
   rbpn         double[3][3]    classical NPB matrix (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix operates in the sense V(date) = rbpn * V(GCRS), where
   the p-vector V(date) is with respect to the true equatorial triad
   of date date1+date2 and the p-vector V(GCRS) is with respect to
   the Geocentric Celestial Reference System (IAU, 2000).

   3. A faster, but slightly less accurate result (about 1 mas), can be
   obtained by using instead the iauPnm00b function.

Called:
   iauPn00a     bias/precession/nutation, IAU 2000A

Reference:

   IAU: Trans. International Astronomical Union, Vol. XXIVB;  Proc.
   24th General Assembly, Manchester, UK.  Resolutions B1.3, B1.6.
   (2000)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPnm00a(double date1, double date2, double rbpn[3][3])

function iauPnm00a(date1::Real, date2::Real)
   rbpn = zeros(Float64, 3, 3)

   ccall((:iauPnm00a, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, date1),
       convert(Float64, date2),
       rbpn)

   return SMatrix{3,3}(rbpn')
end


text:
module YaoQASM
include("grammar.jl")
include("yaoapi.jl")
include("transform.jl")

end # module


text:
export auc, RocCurve, PrecisionRecallCurve, ratio_at_threshold, fpr_at_threshold, tpr_at_threshold, auc_at_threshold, precision_at_threshold, recall_at_threshold

using LinearAlgebra
using Statistics

function auc(x,y)
    direction = 1
    dx = [x[i] - x[i-1] for i in 2:length(x)]
    if any(x-> x < 0, dx)
        if all(x-> x < 0, dx)
            direction = -1
        else
            throw(DomainError("Function is neither increasing nor decreasing"))
        end
    end
    return direction * _trapz(y,x)
end

auc(c::confusion_matrix) = auc(false_positive_rate(c; average="binary"), true_positive_rate(c; average="binary"))

struct RocCurve
    thresholds::Vector{Number}
    ratios::Vector{Tuple{Number, Number}}
    options::NamedTuple
end

function RocCurve(y_true, y_pred, thresholding_func, thresholds; ith_class = nothing, class_name = nothing, average = "macro", normalize = false, weights = nothing)
    return _calculate_for_curves(y_true, y_pred, thresholding_func, thresholds; ith_class = ith_class, class_name = class_name,
        average = average, normalize = normalize, weights = weights)
end

function RocCurve(c::Vector{confusion_matrix} , thresholds, first_func = false_positive_rate, second_func = true_positive_rate; ith_class = nothing,
        class_name = nothing, average = "macro", normalize = false, weights = nothing)
        return _calculate_for_curves_with_matrices(c, thresholding_func, thresholds; ith_class = ith_class, class_name = class_name,
            average = average, normalize = normalize, weights = weights)
end

ratio_at_threshold(r::RocCurve, threshold) = r.ratios[findfirst(x -> x == threshold, r.thresholds)]
fpr_at_threshold(r::RocCurve, threshold) = ratio_at_threshold(r, threshold)[1]
tpr_at_threshold(r::RocCurve, threshold) = ratio_at_threshold(r, threshold)[2]

function Base.show(io::IO, ::MIME"text/plain", r::RocCurve)
    println(io,"Roc Curve:")
    println(io, repeat("=", 60))
    for i in 1:length(r.thresholds)
        println(io, lpad("Threshold : ", 30), r.thresholds[i],  lpad("False Positive Rate : ", 30), r.ratios[i][1], lpad("True Positive Rate : ", 30), r.ratios[i][2])
    end
    println("These values have been calculated with the following parameters : ", r.options)
end

##

struct PrecisionRecallCurve
    thresholds::Vector{Number}
    ratios::Vector{Tuple{Number, Number}}
    options::NamedTuple
end

function PrecisionRecallCurve(y_true, y_pred, thresholding_func, thresholds; ith_class = nothing, class_name = nothing, average = "macro", normalize = false, weights = nothing)
    return _calculate_for_curves(y_true, y_pred, thresholding_func, thresholds, precision_score, recall_score; ith_class = ith_class, class_name = class_name,
        average = average, normalize = normalize, weights = weights)
end

function PrecisionRecallCurve(c::Vector{confusion_matrix} , thresholds, first_func = precision_score, second_func = recall_score; ith_class = nothing,
        class_name = nothing, average = "macro", normalize = false, weights = nothing)
        return _calculate_for_curves_with_matrices(c, thresholding_func, thresholds, precision_score, recall_score; ith_class = ith_class, class_name = class_name,
            average = average, normalize = normalize, weights = weights)
end

ratio_at_threshold(r::PrecisionRecallCurve, threshold) =  r.ratios[findfirst(x -> x == threshold, r.thresholds)]
precision_at_threshold(r::PrecisionRecallCurve, threshold) = ratio_at_threshold(r, threshold)[1]
recall_at_threshold(r::PrecisionRecallCurve, threshold) = ratio_at_threshold(r, threshold)[2]

function Base.show(io::IO, ::MIME"text/plain", r::PrecisionRecallCurve)
    println(io,"Precision Recall Curve:")
    println(io, repeat("=", 60))
    for i in 1:length(r.thresholds)
        println(io, lpad("Threshold : ", 30), r.thresholds[i],  lpad("Precision Score : ", 30), r.ratios[i][1], lpad("Recall Score : ", 30), r.ratios[i][2])
    end
    println("These values have been calculated with the following parameters : ", r.options)
end

##

#TODO
# roc_auc_score and other metrics


text:
using Mill
using SparseArrays, DataFrames

import Mill: sparsify, mapdata
import LearnBase: nobs

@testset "creating bags" begin
    k = [2, 2, 2, 1, 1, 3]
    @test Mill.bags(k).bags == [1:3,4:5,6:6]
end

metadata = fill("metadata", 4)
a = BagNode(ArrayNode(rand(3,4)),[1:4], metadata)
b = BagNode(ArrayNode(rand(3,4)),[1:2,3:4], metadata)
c = BagNode(ArrayNode(rand(3,4)),[1:1,2:2,3:4], metadata)
d = BagNode(ArrayNode(rand(3,4)),[1:4,0:-1], metadata)
wa = WeightedBagNode(ArrayNode(rand(3,4)),[1:4], rand(1:4, 4), metadata)
wb = WeightedBagNode(ArrayNode(rand(3,4)),[1:2,3:4], rand(1:4, 4), metadata)
wc = WeightedBagNode(ArrayNode(rand(3,4)),[1:1,2:2,3:4], rand(1:4, 4), metadata)
wd = WeightedBagNode(ArrayNode(rand(3,4)),[1:4,0:-1], rand(1:4, 4), metadata)
e = ArrayNode(rand(2, 2))

f = ProductNode((wb,b))
g = ProductNode([c, wc])
h = ProductNode((wc,c))
i = ProductNode((
              b,
              ProductNode((
                        b,
                        BagNode(
                                BagNode(
                                        ArrayNode(rand(2,4)),
                                        [1:1, 2:2, 3:3, 4:4]
                                       ),
                                [1:3, 4:4]
                               )
                       ))
             ))
k = ProductNode((a = wb, b = b))
l = ProductNode((a = wc, b = c))

@testset "testing nobs" begin
    @test nobs(a) == nobs(wa) == 1
    @test nobs(b) == nobs(wb) == 2
    @test nobs(c) == nobs(wc) == 3
    @test nobs(d) == nobs(wd) == 2
    @test nobs(e) == 2
    @test nobs(f) == nobs(wb) == nobs(b) == nobs(k)
    @test nobs(g) == nobs(c) == nobs(wc) == nobs(l)
    @test nobs(h) == nobs(wc) == nobs(c)
    @test nobs(i) == nobs(b)
end

@testset "testing nobs" begin
    @test nobs(a) == nobs(wa) == 1
    @test nobs(b) == nobs(wb) == 2
    @test nobs(c) == nobs(wc) == 3
    @test nobs(d) == nobs(wd) == 2
    @test nobs(e) == 2
    @test nobs(f) == nobs(wb) == nobs(b) == nobs(k)
    @test nobs(g) == nobs(c) == nobs(wc) == nobs(l)
    @test nobs(h) == nobs(wc) == nobs(c)
    @test nobs(i) == nobs(b)
end

@testset "testing ArrayNode hcat and vcat" begin
    @test all(catobs(e, e).data .== hcat(e.data, e.data) .== reduce(catobs, [e,e]).data)
    @test all(hcat(e, e).data .== hcat(e.data, e.data))
    @test all(vcat(e, e).data .== vcat(e.data, e.data))
    x = ArrayNode(randn(2,3),rand(3))
    @test catobs(x,x[0:-1]) isa ArrayNode{Array{Float64,2},Array{Float64,1}}
    @inferred catobs(x,x[0:-1]) == [1,2,3]
    @test reduce(catobs, [x, x[0:-1]]) isa ArrayNode{Array{Float64,2},Array{Float64,1}}
    @inferred reduce(catobs, [x, x[0:-1]]) == [1,2,3]
    @test all(cat(e, e, dims = ndims(e)).data .== hcat(e.data, e.data))
end

@testset "testing BagNode hcat" begin
    @test all(catobs(a, b, c).data.data .== hcat(a.data.data, b.data.data, c.data.data))
    @test all(reduce(catobs, [a, b, c]).data.data .== hcat(a.data.data, b.data.data, c.data.data))
    @test catobs(a, b, c).bags.bags == vcat(a.bags, b.bags, c.bags).bags
    @test reduce(catobs, [a, b, c]).bags.bags == vcat(a.bags, b.bags, c.bags).bags
    @test all(catobs(c, a).data.data .== hcat(c.data.data, a.data.data))
    @test all(reduce(catobs, [c, a]).data.data .== hcat(c.data.data, a.data.data))
    @test catobs(c, a).bags.bags == vcat(c.bags, a.bags).bags
    @test reduce(catobs, [c, a]).bags.bags == vcat(c.bags, a.bags).bags
    @test all(catobs(a, c).data.data .== hcat(a.data.data, c.data.data))
    @test all(reduce(catobs, [a, c]).data.data .== hcat(a.data.data, c.data.data))
    @test catobs(a, c).bags.bags == vcat(a.bags, c.bags).bags
    @test reduce(catobs, [a, c]).bags.bags == vcat(a.bags, c.bags).bags
    @test all(catobs(a, d).data.data .== hcat(a.data.data, d.data.data))
    @test all(reduce(catobs, [a, d]).data.data .== hcat(a.data.data, d.data.data))
    @test catobs(a, d).bags.bags == vcat(a.bags, d.bags).bags
    @test reduce(catobs, [a, d]).bags.bags == vcat(a.bags, d.bags).bags
    @test all(catobs(d, a).data.data .== hcat(d.data.data, a.data.data))
    @test all(reduce(catobs, [d, a]).data.data .== hcat(d.data.data, a.data.data))
    @test catobs(d, a).bags.bags == vcat(d.bags, a.bags).bags
    @test reduce(catobs, [d, a]).bags.bags == vcat(d.bags, a.bags).bags
    @test all(catobs(d).data.data .== hcat(d.data.data))
    @test all(reduce(catobs, [d]).data.data .== hcat(d.data.data))
    @test catobs(d).bags.bags == vcat(d.bags).bags
    @test reduce(catobs, [d]).bags.bags == vcat(d.bags).bags
    @test all(cat(a, b, dims = ndims(a)).data.data .== hcat(a.data.data, b.data.data))
end

@testset "testing WeightedBagNode hcat" begin
    @test all(catobs(wa, wb, wc).data.data .== hcat(wa.data.data, wb.data.data, wc.data.data))
    @test all(reduce(catobs, [wa, wb, wc]).data.data .== hcat(wa.data.data, wb.data.data, wc.data.data))
    @test catobs(wa, wb, wc).bags.bags == vcat(wa.bags, wb.bags, wc.bags).bags
    @test reduce(catobs, [wa, wb, wc]).bags.bags == vcat(wa.bags, wb.bags, wc.bags).bags
    @test all(catobs(wc, wa).data.data .== hcat(wc.data.data, wa.data.data))
    @test all(reduce(catobs, [wc, wa]).data.data .== hcat(wc.data.data, wa.data.data))
    @test catobs(wc, wa).bags.bags == vcat(wc.bags, wa.bags).bags
    @test reduce(catobs, [wc, wa]).bags.bags == vcat(wc.bags, wa.bags).bags
    @test all(catobs(wa, wc).data.data .== hcat(wa.data.data, wc.data.data))
    @test all(reduce(catobs, [wa, wc]).data.data .== hcat(wa.data.data, wc.data.data))
    @test catobs(wa, wc).bags.bags == vcat(wa.bags, wc.bags).bags
    @test reduce(catobs, [wa, wc]).bags.bags == vcat(wa.bags, wc.bags).bags
    @test all(catobs(wa, wd).data.data .== hcat(wa.data.data, wd.data.data))
    @test all(reduce(catobs, [wa, wd]).data.data .== hcat(wa.data.data, wd.data.data))
    @test catobs(wa, wd).bags.bags == vcat(wa.bags, wd.bags).bags
    @test reduce(catobs, [wa, wd]).bags.bags == vcat(wa.bags, wd.bags).bags
    @test all(catobs(wd, wa).data.data .== hcat(wd.data.data, wa.data.data))
    @test all(reduce(catobs, [wd, wa]).data.data .== hcat(wd.data.data, wa.data.data))
    @test catobs(wd, wa).bags.bags == vcat(wd.bags, wa.bags).bags
    @test reduce(catobs, [wd, wa]).bags.bags == vcat(wd.bags, wa.bags).bags
    @test all(catobs(wd).data.data .== hcat(wd.data.data))
    @test all(reduce(catobs, [wd]).data.data .== hcat(wd.data.data))
    @test catobs(wd).bags.bags == vcat(wd.bags).bags
    @test reduce(catobs, [wd]).bags.bags == vcat(wd.bags).bags
end

@testset "testing catobs with missing values" begin
    @test catobs(a, b, c).data.data == catobs(a, b, missing, c).data.data
    @test catobs(a, b, c).bags.bags == catobs(a, b, missing, c).bags.bags
    @test catobs(wa, wb, wc).data.data == catobs(wa, wb, missing, wc).data.data
    @test catobs(wa, wb, wc).bags.bags == catobs(wa, wb, missing, wc).bags.bags
    @test catobs(wa, wb, wc).weights == catobs(wa, wb, missing, wc).weights
end

@testset "testing hierarchical hcat on product nodes" begin
    @test all(catobs(f, h).data[1].data.data .== hcat(wb.data.data, wc.data.data))
    @test all(reduce(catobs, [f, h]).data[1].data.data .== hcat(wb.data.data, wc.data.data))
    @test all(catobs(f, h).data[2].data.data .== hcat(b.data.data, c.data.data))
    @test all(reduce(catobs, [f, h]).data[2].data.data .== hcat(b.data.data, c.data.data))
    @test all(catobs(f, h, f).data[1].data.data .== hcat(wb.data.data, wc.data.data, wb.data.data))
    @test all(reduce(catobs, [f, h, f]).data[1].data.data .== hcat(wb.data.data, wc.data.data, wb.data.data))

    @test all(catobs(k, l).data[1].data.data .== hcat(wb.data.data, wc.data.data))
    @test all(catobs(k, l).data[2].data.data .== hcat(b.data.data, c.data.data))
    @test nobs(catobs(k,l)) == nobs(k) + nobs(l)

    # different tuple length
    @test_skip @test_throws Exception catobs(f, i)
    @test_skip @test_throws Exception reduce(catobs, [f, i])
    @test_skip @test_throws Exception catobs(g, i)
    @test_skip @test_throws Exception reduce(catobs, [g, i])
    # different content
    @test_skip @test_throws Exception catobs(f, g)
    @test_skip @test_throws Exception reduce(catobs, [f, g])
end

@testset "testing BagNode indexing" begin
    @test all(a[1].data.data .== a.data.data)
    @test a[1].bags.bags == [1:4]
    @test all(b[1:2].data.data .== b.data.data)
    @test b[1:2].bags.bags == [1:2,3:4]
    @test all(b[2].data.data .== b.data.data[:,3:4])
    @test b[2].bags.bags == [1:2]
    @test all(b[1].data.data .== b.data.data[:,1:2])
    @test b[1].bags.bags == [1:2]
    @test all(c[1].data.data .== c.data.data[:,1:1])
    @test c[1].bags.bags == [1:1]
    @test all(c[[1,2]].data.data .== c.data.data[:,1:2])
    @test c[[1,2]].bags.bags == [1:1,2:2]
    @test all(c[[2,1]].data.data .== c.data.data[:,[2,1]])
    @test c[[2,1]].bags.bags == [1:1,2:2]
    @test all(d[[2,1]].data.data .== d.data.data)
    @test d[[2,1]].bags.bags == [0:-1,1:4]
    @test all(d[1:2].data.data .== d.data.data)
    @test d[1:2].bags.bags == [1:4,0:-1]
    @test nobs(d[2].data) == 0
    @test d[2].bags.bags == [0:-1]
    @test isempty(a[2:1].bags.bags)
    @test nobs(a[2:1].data) == 0
end

@testset "testing WeightedBagNode indexing" begin
    @test all(wa[1].data.data .== wa.data.data)
    @test wa[1].bags.bags == [1:4]
    @test all(wb[1:2].data.data .== wb.data.data)
    @test wb[1:2].bags.bags == [1:2,3:4]
    @test all(wb[2].data.data .== wb.data.data[:,3:4])
    @test wb[2].bags.bags == [1:2]
    @test all(wb[1].data.data .== wb.data.data[:,1:2])
    @test wb[1].bags.bags == [1:2]
    @test all(wc[1].data.data .== wc.data.data[:,1:1])
    @test wc[1].bags.bags == [1:1]
    @test all(wc[[1,2]].data.data .== wc.data.data[:,1:2])
    @test wc[[1,2]].bags.bags == [1:1,2:2]
    @test all(wc[[2,1]].data.data .== wc.data.data[:,[2,1]])
    @test wc[[2,1]].bags.bags == [1:1,2:2]
    @test all(wd[[2,1]].data.data .== wd.data.data)
    @test wd[[2,1]].bags.bags == [0:-1,1:4]
    @test all(wd[1:2].data.data .== wd.data.data)
    @test wd[1:2].bags.bags == [1:4,0:-1]
    @test nobs(wd[2].data) == 0
    @test wd[2].bags.bags == [0:-1]
end

@testset "testing nested ragged array" begin
    x = BagNode(ArrayNode(rand(3,10)),[1:2,3:3,0:-1,4:5,6:6,7:10])
    y = BagNode(x,[1:2,3:3,4:5,6:6])
    @test all(y[1].data.data.data .== x.data.data[:,1:3])
    @test y[1].data.bags.bags == [1:2,3:3]
    @test all(y[1:2].data.data.data .== x.data.data[:,1:3])
    @test y[1:2].data.bags.bags == [1:2,3:3,0:-1]
    @test all(y[2:3].data.data.data .== x.data.data[:,4:6])
    @test y[2:3].data.bags.bags == [0:-1,1:2,3:3]
end

@testset "testing ProductNode" begin
    x = ProductNode((ArrayNode(rand(3,2)),ArrayNode(rand(3,2)),ArrayNode(randn(3,2))))
    y = ProductNode((ArrayNode(rand(3,2)),ArrayNode(rand(3,2)),ArrayNode(randn(3,2))))
    @test all(catobs(x,y).data[1].data .== hcat(x.data[1].data,y.data[1].data))
    @test all(reduce(catobs, [x,y]).data[1].data .== hcat(x.data[1].data,y.data[1].data))
    @test all(catobs(x,y).data[2].data .== hcat(x.data[2].data,y.data[2].data))
    @test all(reduce(catobs, [x,y]).data[2].data .== hcat(x.data[2].data,y.data[2].data))
    @test all(catobs(x,y).data[3].data .== hcat(x.data[3].data,y.data[3].data))
    @test all(reduce(catobs, [x,y]).data[3].data .== hcat(x.data[3].data,y.data[3].data))
    @test all(cat(x,y, dims = ndims(x)).data[3].data .== hcat(x.data[3].data,y.data[3].data))

    @test all(k[1].data[1].data.data .==  wb[1].data.data)
    @test all(k[2].data[1].data.data .==  wb[2].data.data)
    @test all(l[2:3].data[1].data.data .==  wc[2:3].data.data)
end

@testset "testing sparsify" begin
    @test sparsify(zeros(10, 10), 0.05) isa SparseMatrixCSC
    @test sparsify(randn(10, 10), 0.05) isa Matrix
    @test sparsify(randn(10), 0.05) isa Vector
end

@testset "testing sparsify and mapdata" begin
    x = ProductNode((ProductNode((ArrayNode(randn(5,5)), ArrayNode(zeros(5,5)))), ArrayNode(zeros(5,5))))
    xs = mapdata(i -> sparsify(i, 0.05), x)
    @test xs.data[2].data isa SparseMatrixCSC
    @test xs.data[1].data[2].data isa SparseMatrixCSC
    @test all(xs.data[1].data[1].data .== x.data[1].data[1].data)
end

@testset "testing missing mapdata" begin
    x = ProductNode((ProductNode((ArrayNode(randn(5,5)), ArrayNode(zeros(5,5)))), ArrayNode(zeros(5,5))), BagNode(missing, AlignedBags([0:-1]), nothing))
    xs = mapdata(i -> sparsify(i, 0.05), x)
    @test xs.data[2].data isa SparseMatrixCSC
    @test xs.data[1].data[2].data isa SparseMatrixCSC
    @test all(xs.data[1].data[1].data .== x.data[1].data[1].data)
end

@testset "testing terseprint on" begin
    Mill.terseprint(true)
    buf = IOBuffer()
    show(buf, typeof([h,i]))
    str_repr = String(take!(buf))
    @test str_repr == "Array{ProductNode{…},1}"
end

@testset "testing terseprint off" begin
    Mill.terseprint(false)
    buf = IOBuffer()
    show(buf, typeof([h,i]))
    str_repr = String(take!(buf))
    @test str_repr == "Array{ProductNode{T,Nothing} where T,1}"
end

@testset "testing equals and hash" begin
    a2 = deepcopy(a)
    i2 = deepcopy(i)
    k2 = ProductNode((a = wb, b = b))
    metadata1 = fill("metadata", 4)
    metadata2 = "Oh, Hi Mark"
    r = rand(3,4)
    x1 = BagNode(ArrayNode(r),[1:4], metadata1)
    x2 = BagNode(ArrayNode(r),[1:4], metadata1)
    x3 = BagNode(ArrayNode(r),[1:4], metadata2)

    @test a != b
    @test a != i
    @test a == a2
    @test i == i2
    @test i != k
    @test k == k2
    @test x1 == x2
    @test x1 != x3
    @test hash(a) !== hash(b)
    @test hash(a) === hash(a2)
    @test hash(a) !== hash(i)
    @test hash(i) === hash(i2)
    @test hash(i) !== hash(k)
    @test hash(k) === hash(k2)
    @test hash(x1) === hash(x2)
    @test hash(x1) !== hash(x3)
end

function buf_show(data; kwargs...)
    buf = IOBuffer()
	Base.show(buf, data; kwargs...)
	String(take!(buf))
end

@testset "testing more obscure terseprint things" begin
    function experiment(ds::LazyNode{T}) where {T<:Symbol}
    	@show ds
    	@show T
    end
    a = methods(experiment)
    a_method = a.ms[1]
    b = getfield(a_method.sig, 2)
    c = getfield(b, 3)
    d = c[2]
    e = d.body

    t = UnionAll(TypeVar(:t), LazyNode)
    u = typeof(LazyNode(:oh_hi, ["Mark"]))

    orig_terse = Mill._terseprint[]

    Mill.terseprint(true)
    @test occursin("(ds::LazyNode{…}) where T<:Symbol", buf_show(methods(experiment)))
    @test buf_show(t) == "LazyNode{…}"
    @test buf_show(u) == "LazyNode{…}"
    @test buf_show(d) == "LazyNode{…}"
    @test buf_show(e) == "LazyNode{…}"

    # extremely weird behavior, see https://github.com/pevnak/Mill.jl/issues/45
	Mill.terseprint(false)
    @test_throws ErrorException startswith("(ds::LazyNode{T,D} where D) where T<:Symbol", buf_show(methods(experiment)))
    @test_broken buf_show(methods(experiment))
    @test buf_show(t) == "LazyNode"
    @test buf_show(u) == "LazyNode{:oh_hi,Array{String,1}}"
    @test_throws ErrorException buf_show(d)
    @test_broken buf_show(d) == "LazyNode{T<:Symbol,D}"
	@test_throws ErrorException buf_show(e)
    @test_broken buf_show(e) == "LazyNode{T<:Symbol,D}"

    Mill.terseprint(orig_terse)
end


text:
function encode(input)

end

function decode(input)

end


text:
# KPreso with DVC

include("misc.jl")
include("primitives.jl")
include("kptype.jl")
include("kpbase.jl")

# variable manipulation

const SEPSYM='_'

issym(x) = isa(x, Symbol) || isa(x, Number)
iscid(x) = issym(x) && string(x)[1] == 'C'
isrid(x) = issym(x) && string(x)[1] == 'R'
hasid(x) = occursin(SEPSYM, string(x))

function newvar(xid, var)
 if isuppercase(string(var)[1])
   return var
 elseif hasid(var)
   return Symbol(var, xid)
 else
   return Symbol(var, SEPSYM, xid)
 end
end

function newvarlist(xid, varlist)
 nvars = map(var->newvar(xid, var), varlist)
end

# Resolvent

function rename_term(xid, vars, term::Term)
  if issym(term)
    if isevar(term, vars)
      return term
    elseif isvar(term, vars)
      return newvar(xid, term)
    else
      return term
    end
  else
    params = term.args
    nterm = deepcopy(term)
    nterm.args = map(param->rename_term(xid, vars, param), params)
    return nterm
  end
end

function rename_atom(xid, vars, atom::KPAtom)
 for v in keys(atom.args)
  atom.args[v] = rename_term(xid, vars, atom.args[v])
 end
 atom
end

function rename_literal(xid, vars, lit::KPLiteral) #::KPLiteral)
  lit.atom = rename_atom(xid, vars, lit.atom)
  lit
end

function rename_clause(xid, vars, body)
 for ix in 1:length(body)
  body[ix] = rename_literal(xid, vars, body[ix])
 end
 nvars=newvarlist(xid, vars)
# CForm2(xid, newvarlist(xid, vars), map(lit->rename_term(xid, vars, lit), body))
 cf=CForm2(xid, nvars, map(lit->rename_literal(xid, vars, lit), body))
 cf
end

###
function entrylit(rid, nlid, lid, core)
  olit = literalof(lid, core).body
  ovars = varsof(cidof(lid,core),core)
  nlit=LForm2(nlid, rename_literal(rid, ovars, olit))
#  nlit=LForm2(nlid, rename_term(rid, ovars, olit))
  core.ldb[nlid] = nlit
end

function rename_lid(rid, lid, core)
  Symbol(origof(lid), SEPSYM, rid)
end

function rename_lids(rid, lids, core)
 nlids=[]
 for lid in lids
  nlid=rename_lid(rid, lid, core)
  entrylit(rid, nlid, lid, core)
  push!(nlids, nlid)
 end
 return nlids
end

## for renamed vars explosion
#fitting vars in literal
function fitting_vars_term(vars, term)
  if issym(term)
    if isvar(term, vars)
      return [term]
    else
      return []
    end
  else
    return fitting_vars_args(vars, term.args[2:end])
  end
end

function fitting_vars_args(vars, args)
  nvars=[]
  for term in args
    append!(nvars, fitting_vars_term(vars, term))
  end
  union(nvars, nvars) # uniquefy
end

function fitting_vars_lit(vars, lit)
  fitting_vars_term(vars, lit.args[2])
end

function fitting_vars_atom(vars, atom)
  evars=[]
  for k in atom.args
   append!(evars, fitting_vars_term(vars, k[2]))
  end
  evars
end

function fitting_vars(vars, lits, core)
 evars = []
 for lit in lits
   if lit isa KPLiteral
     avars = fitting_vars_atom(vars, lit.atom)
     append!(evars, avars)
   else
    append!(evars, fitting_vars_lit(vars, lit))
   end
 end
 nvars=union(evars, evars)
 return nvars
end

## evaluation

function leval(lit)
  try
    val=eval(lit.args[2])
    if val==true
      if lit.args[1] == :+; return true end
      if lit.args[1] == :-; return false end
    elseif val==false
      if lit.args[1] == :-; return true end
      if lit.args[1] == :+; return false end
    end
    return lit
  catch
    return lit
  end
end

function evaluate_literals(lids, lits)
 rlits=[]
 rlids=[]
 for lix in 1:length(lids)
   lid=lids[lix]
   lit=lits[lix]
   val = leval(lit)
   if val == true; return true,true end
   if val == false; continue end
   push!(rlits,lit)
   push!(rlids,lid)
 end
 return rlids,rlits
end

## resolution
ovarsof(l1,l2,core)=vcat(varsof(cidof(l1, core), core), varsof(cidof(l2, core), core))

function kp_resolution(l1,l2,core)
 vars1 = varsof(cidof(l1, core), core)
 vars2 = varsof(cidof(l2, core), core)
 lit1  = literalof(l1, core).body
 lit2  = literalof(l2, core).body
 ovars = vcat(vars1,vars2)
 (sign1, psym1) = psymof(l1, core)
 (sign2, psym2) = psymof(l2, core)

 if sign1 == sign2; return :FAIL end
 if psym1 != psym2; return :FAIL end
 if length(lit1.atom.args) != length(lit2.atom.args); return :FAIL end
# if length(lit1.args[2].args) != length(lit2.args[2].args); return :FAIL end

 try
   core.trycnt[1] += 1
   sigmai = kpunify(ovars, lit1.atom, lit2.atom)
#   sigmai = kpunify(ovars, lit1.args[2], lit2.args[2])
   rem1 = lidsof(cidof(l1, core),core)
   rem1 = setdiff(rem1, [l1])
   rem2 = lidsof(cidof(l2, core),core)
   rem2 = setdiff(rem2, [l2])

   rem = vcat(rem1, rem2)
   vars = ovars

# rename rlid
   rid =  newrid(core)
   nrem = rename_lids(rid, rem, core)
   nbody = literalsof(rem, core)
   nbody1 = kpapply(ovars, nbody, sigmai)
   if evalon
     rb = evaluate_literals(nrem, nbody1)
     if rb[1] == true
       println("Valid")
       return :FAIL
     end
     nrem, nbody1 = rb
   end
  ovars = newvarlist(rid, ovars)
   vars = fitting_vars(ovars, nbody1, core)
   body = rename_clause(rid, vars, nbody1)
 rename_subst = [vars, body.vars]

 ## settlement

 core.succnt[1] += 1

# cdb[rid] to vars
  core.cdb[rid] = VForm2(rid, body.vars)

# clmap[rid] to rlid* = nrem
  core.clmap[rid] = nrem

# ldb[rlid] to full

  for i in 1:length(nrem)
    core.ldb[nrem[i]] = LForm2(nrem[i], body.body[i])
  end

# lcmap[rlid] to rid
  for rlid in nrem
    core.lcmap[rlid] = rid
  end
  core.proof[rid] = STEP(rid, l1, l2, sigmai, rename_subst)
  rcf2=CForm2(rid, body.vars, body.body)
  return rcf2

  catch e
    println("FAIL = $e")
    return :FAIL
  end
end

## Template
function psymof(lid, core)
 lit = literalof(lid, core).body
 (lit.sign, lit.atom.Psym)
end

lsym(sign, psym) = Symbol(sign, psym)

function lsymof(lid, core)
 psym = psymof(lid, core)
 string(psym[1],psym[2])
end

inverseof(sign) = if sign == :- ; :+ else :- end

function litis(sign, psym, lid, core)
 (osign, opsym) = psymof(lid, core)
 return osign == inverseof(sign) && opsym == psym
end

function oppositof(sign, psym, core)
 rems = []

 for olid in keys(core.ldb)
  lfm2 = core.ldb[olid]
  if litis(sign, psym, olid, core)
    push!(rems, olid)
  end
 end
 rems

end

function templateof(sign, psym, core)
 eqs = Dict()
 oppos = oppositof(sign, psym, core)
 body = []
 for lid in oppos
  cid = cidof(lid, core)
  lids = lidsof(cid, core)
  rem = setdiff(lids, [lid])
  push!(body, [cid, varsof(cid, core), lid, rem])
 end
 core.level0[lsym(sign, psym)]= body
end

function alltemplateof(core)
 allpsym = core.allpsym
 alltemp = []
 for psym in allpsym
   push!(alltemp, templateof(:+, psym, core))
   push!(alltemp, templateof(:-, psym, core))
 end
 core.level0
end

## apply template to lid
function applytemp(lid, core)
 (sign, psym) = psymof(lid, core)
 templs = templateof(sign, psym, core)
 rids = []
 for templ in templs
   reso = kp_resolution(lid, templ[3], core)
   if isa(reso, CForm2)
     push!(rids, reso.cid)
   else
     continue
   end
 end
 map(rid->lidsof(rid, core),rids)
end

### step forward on template

function prod(tmpl1, tmpl2)
  newtmpl = []
  for t1 in tmpl1
  for t2 in tmpl2
    push!(newtmpl, vcat(t1, t2))
  end
  end
  newtmpl
end

function prodm(tmpls)
 mtempl=tmpls[1]
 for rtmpl in tmpls[2:end]
  mtempl=prod(mtempl, rtmpl)
 end
 return mtempl
end

function rotate(lids, i=0)
 i >= length(lids) && error("required: 0<=$i<length($lids)")
 if i == 0
   lids
 else
   vcat(lids[i+1:end], lids[1:i])
 end
end

"""
 goal = [L1,L2,...]
"""
function dostepagoal(goal, core)
 newlids = []
 for lid in goal
   nlids = applytemp(lid, core)
   append!(newlids, nlids)
 end
 newlids
end

"""
 goals = [[L1,L2,...],...]
"""
function dostep1goals(goals, core)
 nextg = []
 for g in goals
  ngs = dostepagoal(g, core)
  append!(nextg, ngs)
 end
println("$goals => $nextg")
 nextg
end

"""
 goal = [L1,L2,...]
"""
function dostepagoal1(goal, core)
 nlids = []
 if isempty(goal); return [] end
 for shift in 0:length(goal)-1
  wgoal = rotate(goal,shift)
  lid = wgoal[1]
  nlids = applytemp(lid, core)
  if !isempty(nlids); break end
 end
 if isempty(nlids); return :FAIL end
 return nlids
end

"""
 goals = [[L1,L2,...],...]
"""
function dostepgoals1(goals, core)
 nextg = []
 for g in goals
  ngs = dostepagoal1(g, core)
  if ngs == :FAIL; continue end
  append!(nextg, ngs)
 end
println("$goals => $nextg")
 nextg
end

"""
 find the resolvents of body empty
"""
function contradictionsof(core)
 conds=[]
 for cid in keys(core.cdb)
   if length(lidsof(cid, core)) == 0
     push!(conds, cid)
   end
 end
 conds
end

function pairwiseeq(d1,d2)
 all(map((x,y)->origof(x)==origof(y),d1,d2))
end

# proof should be an array of pair(Lid)
function findrepeat(proof)
 for i in 1:(length(proof)-1)
  for j in (i+1):length(proof)
   if pairwiseeq(proof[i], proof[j])
     println("$(proof[i]) is upon is repeated")
     return true
   end
  end
 end
 return false
end

isresolvent(cid) = isrid(cid)

function proofstep(cid, core)
  step = core.proof[cid]
  return [step.leftp, step.rightp]
end

function prooftree(cid, core)
  if !isresolvent(cid);return [] end

  step = core.proof[cid]
  lefttree  = prooftree(cidof(step.leftp,core), core)
  righttree = prooftree(cidof(step.rightp,core), core)

  parenttree = []
  !isempty(lefttree) && append!(parenttree, lefttree)
  !isempty(righttree) && append!(parenttree, righttree)
  append!(parenttree, [proofstep(cid, core)])
  return parenttree
end

function isrepeatproof(cid, core)
  aproof = prooftree(cid, core)
  findrepeat(aproof)
end

### prover
"""
simple prover find some contracictions, but not all
"""
function kpsimpleprover(cnf, steplimit, contralimit)
 println("in KP prover, literal evaluation is not supported. evalon is off")
 global evalon = false

 println("reading core file")
 cdx=readcore(cnf)
 tdx=alltemplateof(cdx)
 gb=[lidsof(:C1, cdx)]
 conds = []
 nstep = 0;

 evalon && evalproc(cdx.proc)

 println("preparation is finished\n")

 while true
  ga=dostepgoals1(gb, cdx)
  nstep += 1
  conds = contradictionsof(cdx)
 # if !isempty(conds); break end
  if length(conds) >= contralimit;break end
  if nstep >= steplimit;break end
  gb = ga
 end
 return conds,cdx
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
function scaledbessel(xij, qi, qj, ti, tj, mu, eta, kappa)

dim, N = size(xij);
rcut = eta[1]
alpha = eta[2]
nbf = length(mu)

dij = sqrt.(xij[1,:].*xij[1,:] .+ xij[2,:].*xij[2,:] .+ xij[3,:].*xij[3,:]);
dr = xij./reshape(dij,(1,N));
x =  (1.0 .- exp.(-alpha*dij/rcut))/(1-exp(-alpha));
dx = -((alpha/rcut)*exp.(-(alpha*dij/rcut)))/(exp(-alpha) - 1);

eij = zeros(N);
fij = zeros(dim, N);
for i = 1:nbf
    a = i*pi;
    eij = eij + mu[i]*(sqrt(2/rcut)/i)*(sin.(a*x)./dij);
    b = (sqrt(2/rcut)/i)*(sin.(a*x)./(dij.^2) .- (a*cos.(a*x).*dx)./dij);
    fij = fij + -mu[i]*dr.*reshape(b,(1,N))
end

return eij, fij

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module EMLUserColumn
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...UserModels.EMLUserTray
	import EMLtypes.length
	@reexport using ...stage_separators.EMLreboiler
	import EMLtypes.length
	@reexport using ...stage_separators.EMLcondenser
	import EMLtypes.length
	@reexport using ...mixers_splitters.EMLsplitter
	import EMLtypes.length
	@reexport using ...stage_separators.EMLtank
	import EMLtypes.length
	@reexport using ...pressure_changers.EMLpump
	import EMLtypes.length
	include("UserColumn/User_Section_ColumnBasic.jl")
	include("UserColumn/User_Section_Column.jl")
	include("UserColumn/User_Distillation_kettle_cond.jl")
	include("UserColumn/User_Distillation_thermosyphon_subcooling.jl")
	include("UserColumn/User_Distillation_thermosyphon_cond.jl")
end


text:
"""
A package for fitting data from auditory evidence accumulation task (Poisson clicks task)
to evidence accumulation model.
"""
__precompile__()

module PBupsModel

# 3rd party
import Base.convert

using MAT
using ForwardDiff
using Optim
using GeneralUtils

import ForwardDiff.DiffBase
# using DiffBase

export

# data_handle
    LoadData,
    TrialData,
    WriteData,

# model_likelihood
    logProbRight,
    LogLikelihood,

# all_trials
    ComputeLL,
    ComputeLL_bbox,
    ComputeGrad,
    ComputeHess,
    TrialsLikelihood,

# optimization
    InitParams,
    ModelFitting,
    FitSummary

include("data_handle.jl")
include("model_likelihood.jl")
include("all_trials.jl")
include("model_optimization.jl")

end # module


text:
###########################################################################
#
#  SerialMC runner: consumes repeatedly a sampler and returns a MCMCChain
#
#
###########################################################################

export SerialMC

println("Loading SerialMC(steps, burnin, thinning) runner")

immutable SerialMC <: MCMCRunner
  burnin::Int
  thinning::Int
  len::Int
  r::Range

  function SerialMC(steps::Range{Int})
    r = steps

    burnin = first(r)-1
    thinning = r.step
    len = last(r)

    @assert burnin >= 0 "Burnin rounds ($burnin) should be >= 0"
    @assert len > burnin "Total MCMC length ($len) should be > to burnin ($burnin)"
    @assert thinning >= 1 "Thinning ($thinning) should be >= 1"

    new(burnin, thinning, len, r)
  end
end

SerialMC(steps::Range1{Int}) = SerialMC(first(steps):1:last(steps))

SerialMC(; steps::Int=100, burnin::Int=0, thinning::Int=1) = SerialMC((burnin+1):thinning:steps)

function run_serialmc(t::MCMCTask)
  tic() # start timer

  # array allocations to store results
  samples        = fill(NaN, t.model.size, length(t.runner.r))
  gradients      = fill(NaN, t.model.size, length(t.runner.r))
  diags          = {"step" => collect(t.runner.r)}
  logtargets     = fill(NaN, length(t.runner.r))

  # sampling loop
  j = 1
  for i in 1:t.runner.len
    newprop = consume(t.task)
    if in(i, t.runner.r)
      samples[:, j] = newprop.ppars
      logtargets[j] = newprop.plogtarget

      if newprop.pgrads != nothing
        gradients[:, j] = newprop.pgrads
      end

      # save diagnostics
      for (k,v) in newprop.diagnostics
        # if diag name not seen before, create column
        if !haskey(diags, k)
          diags[k] = Array(typeof(v), length(diags["step"]))
        end

        diags[k][j] = v
      end

      j += 1
    end
  end

  # create Chain
  MCMCChain(t.runner.r, samples', logtargets, gradients', diags, t, toq())
end

function run_serialmc_exit(t::MCMCTask)
  chain = run_serialmc(t)
  stop!(chain)
  return chain
end

function resume_serialmc(t::MCMCTask; steps::Int=100)
  @assert typeof(t.runner) == SerialMC
    "resume_serialmc can not be called on an MCMCTask whose runner is of type $(fieldtype(t, :runner))"
  run(t.model, t.sampler, SerialMC(steps=steps, thinning=t.runner.thinning))
end


text:
@testset "azimuth.jl" begin
    p1 = rad2deg(Point(0.5, 1.0))
    p2 = rad2deg(Point(1.0, 0.5))
    p3 = rad2deg(Point(0.6, 1.0))  # above p1

    rs = Arc(p1,p2)
    @test azimuth(p1, p2) ≈ rad2deg(5.81412807071) atol = 0.001
    @test azimuth(rs) ≈ rad2deg(5.81412807071) atol = 0.001
    @test azimuth(p3, p1) ≈ 180.0 atol = 0.001

    @test final_azimuth(p1, p2) ≈ rad2deg(5.45864813531) atol = 0.001
    @test final_azimuth(rs) ≈ rad2deg(5.45864813531) atol = 0.001

    angular_distance_x = 34.49
    angular_distance_y = 28.11
    angular_distance_z = 29.0

    @test spherical_angle(angular_distance_z, angular_distance_y, angular_distance_x) ≈ 76.6 atol = 0.1
    @test spherical_angle(angular_distance_x, angular_distance_z, angular_distance_y) ≈ 54.1 atol = 0.1
    @test spherical_angle(angular_distance_y, angular_distance_x, angular_distance_z) ≈ 56.4 atol = 0.1

    @test spherical_excess(angular_distance_z, angular_distance_y, angular_distance_x) ≈ 7.1 atol = 0.1

    angular_distance_x2 = 35.0
    angular_distance_y2 = 35.0
    angular_distance_z2 = 0.1

    @test spherical_angle(angular_distance_z2, angular_distance_y2, angular_distance_x2) ≈ 89.9 atol = 0.1
    @test spherical_angle(angular_distance_x2, angular_distance_z2, angular_distance_y2) ≈ 89.9 atol = 0.1
    @test spherical_angle(angular_distance_y2, angular_distance_x2, angular_distance_z2) ≈ 0.174 atol = 0.001

end


text:
module Retester

end # module


text:
# Snooping functions
function _snoopi_bench_cmd(snoop_script)
    tmin = 0.0 # For benchmarking
    return quote
        global SnoopCompile_ENV = true
        using SnoopCompileCore

        data = @snoopi tmin=$tmin begin
            $snoop_script
        end

        global SnoopCompile_ENV = false

        using CompileBot: timesum
        @info( "\nTotal inference time (ms): \t" * string(timesum(data, :ms)))
    end
end

function _snoopv_bench_cmd(snoop_script, package_name)
    return quote
        using $package_name
        @info("Script execution time:")
        @timev begin
            $snoop_script
        end
    end
end

function _snoop_bench(config::BotConfig, snoop_script::Expr, test_modul::Module = Main; snoop_mode::Symbol)

    package_name = config.package_name
    package_path = config.package_path

    # automatic (based on Julia version)
    if snoop_mode == :auto
        if VERSION < v"1.2"
            snoop_mode = :run_time
        else
            snoop_mode = :snoopi
        end
    end

    if snoop_mode == :snoopi
        snooping_code = _snoopi_bench_cmd(snoop_script)
    elseif snoop_mode == :run_time
        snooping_code = _snoopv_bench_cmd(snoop_script, package_name)
    else
        error("snoop_mode $snoop_mode is unkown")
    end
    snooping_code = toplevel_string(snooping_code)

    ################################################################
    if in(get(ENV, "SnoopCompile_coverage_ENV", false), [true, "true"])
        code_coverage = "user"
    else
        code_coverage = "none"
    end

    julia_cmd = `julia --code-coverage=$code_coverage --project=@. -e $snooping_code`

    addpkg_ifnotfound(:SnoopCompileCore, test_modul)
    out = quote
        package_sym = Symbol($package_name)
        ################################################################
        using CompileBot

        @info("""------------------------
        Benchmark Started
        ------------------------
        """)
        ################################################################
        @info("""------------------------
        Precompile Deactivated Benchmark
        ------------------------
        """)
        CompileBot.precompile_deactivator($package_path);
        ### Log the compiles
        run($julia_cmd)
        ################################################################
        @info("""------------------------
        Precompile Activated Benchmark
        ------------------------
        """)
        CompileBot.precompile_activator($package_path);
        ### Log the compiles
        run($julia_cmd)
        @info("""------------------------
        Benchmark Finished
        ------------------------
        """)
    end
    return out
end

function _snoop_bench(config::BotConfig, test_modul::Module = Main; snoop_mode::Symbol)

    package_name = config.package_name
    package_rootpath = dirname(dirname(config.package_path))

    package = Symbol(package_name)
    runtestpath = "$package_rootpath/test/runtests.jl"

    snoop_script = quote
        using $(package);
        include($runtestpath);
    end
    out = _snoop_bench(config, snoop_script, test_modul; snoop_mode = snoop_mode)
    return out
end

################################################################
"""
    snoop_bench(config::BotConfig, path_to_example_script::String, test_modul::Module=Main; snoop_mode=:auto)

Benchmark the impact of precompile statements, by running a script with and without the precompiles.
`config` can be generated by [`BotConfig`](@ref).
`path_to_example_script` is preferred to be an absolute path.
The example script will be run in the module specified by `test_modul`.
`snoop_mode` can be `:auto`, `:snoopi` (to test with [`SnoopCompileCore.@snoopi`](@ref)),
or `:runtime` (to measure total script run time with `@timev`).
`:auto` chooses `:snoopi` on supported versions of Julia.

See the [online documentation](https://timholy.github.io/SnoopCompile.jl/stable/bot/)
for a more complete overview.

# Extended help

## Example

In this case, the benchmarking script is placed in the same directory as the
precompile script, so we can use `@__DIR__` to find it:

```julia
using CompileBot

snoop_bench(BotConfig("MatLang"), "\$(@__DIR__)/example_script.jl")
```

As an alternative to `@__DIR__` (for example, if you store your benchmarking script(s) outside
the package itself), you can find the package with [`pathof_noload`](@ref).
"""
function snoop_bench(config::BotConfig, path_to_example_script::String, test_modul::Module=Main; snoop_mode::Symbol=:auto)
    # search for the script! - needed because of confusing paths when referencing pattern_or_file in CI
    path_to_example_script = searchdirsboth([pwd(),dirname(dirname(config.package_path))], path_to_example_script)
    snoop_script = quote
        include($path_to_example_script)
    end
    out = _snoop_bench(config, snoop_script, test_modul; snoop_mode = snoop_mode)
    Core.eval( test_modul, out )
end

"""
    snoop_bench(config::BotConfig, test_modul::Module = Main; snoop_mode::Symbol = :auto)

Benchmark your precompile files using the package's `runtests.jl` file.

During snooping, `snoop_bench` sets the global variable `SnoopCompile_ENV` to `true`.
If needed, your `runtests.jl` can check for the existence and value of this variable to
customize test behavior specifically for snooping.
```
"""
function snoop_bench(config::BotConfig, test_modul::Module = Main; snoop_mode::Symbol = :auto)
    out = _snoop_bench(config, test_modul; snoop_mode = snoop_mode)
    Core.eval( test_modul, out )
end

"""
    snoop_bench(config::BotConfig, expression::Expr, test_modul::Module = Main; snoop_mode::Symbol = :auto)

Benchmark your precompile files by evaluating an expression, for example `:(using MyPackage)`.
Interpolation and macros are not supported.
"""
function snoop_bench(config::BotConfig, snoop_script::Expr, test_modul::Module  = Main; snoop_mode::Symbol = :auto)
    out = _snoop_bench(config, snoop_script, test_modul; snoop_mode = snoop_mode)
    Core.eval( test_modul, out )
end


text:
function _packdensematrix!(A::GBVecOrMat{T}, M::DenseVecOrMat; desc = nothing) where {T}
    desc = _handledescriptor(desc)
    Csize = length(A) * sizeof(T)
    values = Ref{Ptr{Cvoid}}(pointer(M))
    isuniform = false
    libgb.GxB_Matrix_pack_FullC(
        A.p,
        values,
        Csize,
        isuniform,
        desc
    )
    return A
end

function _packcscmatrix!(
    A::GBVecOrMat{T},
    colptr::Vector{libgb.GrB_Index},
    rowidx::Vector{libgb.GrB_Index},
    values::Vector{T};
    desc = nothing
    ) where {T}
    colptr .-= 1
    rowidx .-= 1
    colptrsize = length(colptr) * sizeof(libgb.GrB_Index)
    rowidxsize = length(rowidx) * sizeof(libgb.GrB_Index)
    valsize = length(values) * sizeof(T)
    colptr = Ref{Ptr{libgb.GrB_Index}}(pointer(colptr))
    rowidx = Ref{Ptr{libgb.GrB_Index}}(pointer(rowidx))
    values = Ref{Ptr{Cvoid}}(pointer(values))
    desc = _handledescriptor(desc)

    x = libgb.GxB_Matrix_pack_CSC(
        A,
        colptr,
        rowidx,
        values,
        colptrsize,
        rowidxsize,
        valsize,
        false,
        false,
        desc
    )
    return A
end

function _packcsrmatrix!(
    A::GBVecOrMat{T},
    rowptr::Vector{libgb.GrB_Index},
    colidx::Vector{libgb.GrB_Index},
    values::Vector{T};
    desc = nothing
    ) where {T}
    rowptr .-= 1
    colidx .-= 1
    rowptrsize = length(rowptr) * sizeof(libgb.GrB_Index)
    colidxsize = length(colidx) * sizeof(libgb.GrB_Index)
    valsize = length(values) * sizeof(T)
    rowptr = Ref{Ptr{libgb.GrB_Index}}(pointer(rowptr))
    colidx = Ref{Ptr{libgb.GrB_Index}}(pointer(colidx))
    values = Ref{Ptr{Cvoid}}(pointer(values))
    desc = _handledescriptor(desc)

    libgb.GxB_Matrix_pack_CSC(
        A,
        rowptr,
        colidx,
        values,
        rowptrsize,
        colidxsize,
        valsize,
        false,
        false,
        desc
    )
    return A
end


text:
# comment
println("Hi!")

max_stars_repo_name: monaqa/satysfi-code-printer

max_issues_repo_name: monaqa/satysfi-code-printer

max_forks_repo_name: monaqa/satysfi-code-printer


text:
using CSVFiles
using DataFrames
using DrWatson

ddir = projectdir("plots", "figure_gene_expr", "data")

# dmc_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-06-17_directmc.csv")))
# smc_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-06-21_smc.csv")))
comb_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-10-04.csv")))
dmc_data = comb_data[comb_data.Alg .== "Direct MC", :]
smc_data = comb_data[comb_data.Alg .== "SMC", :]
ti_data = comb_data[(comb_data.Alg .== "TI") .& (comb_data.Duration .== comb_data.DiscreteTimes), :]

# ti_data  = DataFrame(load(joinpath(ddir, "gene-expression_2021-06-17_thermodynamic-integration.csv")))
# ti_data2  = DataFrame(load(joinpath(ddir, "gene-expression_2021-06-18_thermodynamic-integration.csv")))
zechner_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-06-18_zechner.csv")))
gaussian_data = DataFrame(load(joinpath(ddir, "gene-expression_2021-07-19_gaussian.csv")))

smc_estimate = groupby(smc_data, :M)[(;M=128)]
# ti_estimates = ti_data[ti_data.Duration .== ti_data.DiscreteTimes, :]
# ti_estimates2 = ti_data2[ti_data.Duration .== ti_data.DiscreteTimes, :]

using CairoMakie

# presentation_theme = Theme(font = "Noto Sans Regular", linewidth=3)

# fig1 = with_theme(presentation_theme) do
#     f = Figure(resolution=(600, 800))
#     ax = Axis(f[1, 1], xlabel="trajectory duration", ylabel="mutual information (nats)")
#     dmc_d = groupby(dmc_data, :M)[end]
#     dmc_lines = lines!(ax, dmc_d.DiscreteTimes, dmc_d.Mean, label="direct")
#     smc_lines = lines!(ax, smc_estimate.DiscreteTimes, smc_estimate.Mean, color="black")
#     zec_lines = lines!(zechner_data.Duration, zechner_data.PMI, linewidth=2, linestyle=:dash, color=:red)

#     legend = Legend(
#         f[2, 1],
#         [dmc_lines, smc_lines, zec_lines],
#         ["direct", "segment-by-segment", "Duso, Zechner (2019)"],
#         orientation = :horizontal,
#         tellheight = true
#         # [group_color, ti_color, smc_color, approx_entry],
#         # [m_label, ["256 MCMC steps", "4096 MCMC steps"],
#         # ["M = 128"], ["Duso, Zechner (2019)", "Gaussian"]],
#         # ["brute force", "thermodynamic\nintegration", "particle filter", "approximations"],
#         # rowgap=0,
#         # titlegap=2,
#         # groupgap=5,
#         # patchsize=(10,10),
#         # padding=(3,3,3,3),
#         # linewidth=1.0,
#         # gridshalign=:left,
#         # titlehalign=:left,
#         # titlefont="Noto Sans Medium",
#         # framevisible=false
#     )
#     xlims!(ax, 0, 8)
#     ylims!(ax, 0, 20)

#     f
# end
# save(projectdir("plots", "figure_gene_expr", "gene-expr-figure-p.pdf"), fig1, pt_per_unit = 1)

fontsize_theme = Theme(
    font = "Noto Sans Regular",
    fontsize = 8,
    colgap=3,
    linewidth=2,
    markersize=6,
    spinewidth=0.5,
    Axis = (
        spinewidth=0.5,
        xgridwidth=0.5,
        ygridwidth=0.5,
        xtickwidth=0.5,
        ytickwidth=0.5,
        xticksize=3,
        yticksize=3,
        xlabelpadding=0,
        ylabelpadding=5,
        xticklabelpad=0,
        yticklabelpad=3,
    ),
    Legend = (
        framewidth=0.5,
        rowgap=0,
        titlegap=2,
        groupgap=5,
        patchsize=(11,10),
        padding=(3,3,3,3),
        linewidth=1.25,
        gridshalign=:left,
        titlehalign=:left,
        titlefont="Noto Sans SemiBold Italic",
        framevisible=false
    )
)

fig = with_theme(fontsize_theme) do
    f = Figure(resolution=(246, 1.4*146), figure_padding = 2)
    Axis(f[1, 1], xlabel="trajectory duration", ylabel="mutual information (nats)")

    groups = collect(groupby(dmc_data, :M; sort=true))

    m_label = ["M = $(g.M[1])" for g in groups]
    cmap = cgrad(:linear_kgy_5_95_c69_n256, rev=true)
    colors = [cmap[i / length(groups) * 0.7] for i = 1:length(groups)]
    ti_colors = [cgrad(:lajolla)[i] for i = range(0.2, 0.6, length=3)]
    pathmi_color = :orangered

    for (i, group) in enumerate(groups)
        label = "$(group.M[1])"
        lines!(group.DiscreteTimes, group.Mean, linewidth=0.75, label=label, color=colors[i])
    end

    group_color = [LineElement(color=color) for color in colors]
    smc_color = [LineElement(color=:black)]

    lines!(smc_estimate.DiscreteTimes, smc_estimate.Mean, label="SMC", linewidth=1.7, color=:black)
    lines!(zechner_data.Duration, zechner_data.PMI, linewidth=1.5, linestyle=:dash, color=pathmi_color)
    lines!(gaussian_data.DiscreteTimes, gaussian_data.Value, linewidth=1.5, linestyle=:dot, color=:MidnightBlue)

    ti_label = []
    ti_color = []
    for (i, group) in enumerate(groupby(ti_data, :M, sort=true))
        scatter!(group.DiscreteTimes, group.Mean, marker=:cross, color=ti_colors[i])
        push!(ti_label, "M = $(group.M[1])")
        push!(ti_color, MarkerElement(marker=:cross, color=ti_colors[i]))
    end

    approx_entry = [LineElement(linestyle=:dash, color=pathmi_color), LineElement(linestyle=:dot, color=:MidnightBlue)]

    legend = Legend(
        f[1, 2],
        [group_color, ti_color, smc_color, approx_entry],
        [m_label, ti_label,
        ["M = 128"], ["Duso, 𝘦𝘵. 𝘢𝘭. (2019)", "Gaussian"]],
        ["DPWS", "TI-PWS", "RR-PWS", "approximations"],
    )
    xlims!(0.0, 10.2)
    ylims!(0.0, 42)

    inset_lims = (3.3, 8.5)
    lines!([(0,0), (inset_lims[1], 0), inset_lims, (0, inset_lims[2]), (0,0)], color=:black, linewidth=0.75, linestyle=:dot)

    inset_ax = Axis(f[1, 1],
        width=Relative(0.55),
        height=Relative(0.33),
        halign=0.25,
        valign=0.95,
        xtickalign=1.0,
        ytickalign=1.0,
        backgroundcolor=:white,
        xticksize=2,
        xticklabelsize=7,
        xgridvisible=true,
        yticksize=2,
        yticklabelsize=7,
    )

    for (i, group) in enumerate(groups)
        label = "$(group.M[1])"
        lines!(inset_ax, group.DiscreteTimes, group.Mean, linewidth=0.5, label=label, color=colors[i])
    end
    lines!(inset_ax, smc_estimate.DiscreteTimes, smc_estimate.Mean, linewidth=1.0, color=:black)
    lines!(inset_ax, zechner_data.Duration, zechner_data.PMI, linewidth=0.75, linestyle=:dash, color=pathmi_color)
    lines!(inset_ax, gaussian_data.DiscreteTimes, gaussian_data.Value, linewidth=1.5, linestyle=:dot, color=:MidnightBlue)

    for (i, group) in enumerate(groupby(ti_data, :M, sort=true))
        scatter!(inset_ax, group.DiscreteTimes, group.Mean, marker=:cross, color=ti_colors[i])
    end

    xlims!(inset_ax, 0, inset_lims[1])
    ylims!(inset_ax, 0, inset_lims[2])
    translate!(inset_ax.scene, 0, 0, 10)
    translate!(inset_ax.elements[:background], 0, 0, 9)
    translate!(inset_ax.elements[:xgridlines], 0, 0, 9)
    translate!(inset_ax.elements[:ygridlines], 0, 0, 9)

    f
end

save(projectdir("plots", "figure_gene_expr", "gene-expr-figure.pdf"), fig, pt_per_unit = 1)


text:
#RSDeltaSigmaPort: utilities for transfer functions
#-------------------------------------------------------------------------------

"""`y = evalRPoly(roots,x,k=1)`
Compute the value of a polynomial which is given in terms of its roots.
"""
function evalRPoly(roots,x::AbstractArray,k=1)
	y = k * ones(size(x))
	finiteroots = findall(isfinite.(roots))
	roots = roots[finiteroots]
	for i in 1:length(roots)
		y = y.*(x .- roots[i])
	end
	return y
end
evalRPoly(roots,x,k=1) = evalRPoly(roots, [x], k)[1]

"""`h = evalTF(tf,z)`

Evaluates the rational function described by the struct tf
at the point(s) given in the z vector.

TF must be either a zpk object or a struct containing
 - form: 'zp' or 'coeff'
 - zeros,poles,k: if form=='zp'
 - num,den: if form=='coeff'

In Matlab 5, the ss/freqresp() function does nearly the same thing.
"""
function evalTF(tf::ZPKData,z)
	h = tf.k * evalRPoly(tf.z,z) ./ evalRPoly(tf.p,z)
	return h
end

function evalTF(tf,z)
	throw(ErrorException("Transfer function of this type not currently supported."))
#elseif any(strcmp(fieldnames(tf),'form'))
#    if strcmp(tf.form,'zp')
#        h = tf.k * evalRPoly(tf.zeros,z) ./ evalRPoly(tf.poles,z);
#    elseif strcmp(tf.form,'coeff')
#        h = polyval(tf.num,z) ./ polyval(tf.den,z);
#    else
#        fprintf(1,'%s: Unknown form: %s\n', mfilename, tf.form);
#    end
#else	% Assume zp form
#    h = tf.k * evalRPoly(tf.zeros,z) ./ evalRPoly(tf.poles,z);
#end
end

"""`H=evalTFP(Hs, Hz, f)`

Compute the value of a transfer function product Hs*Hz at a frequency f,
where Hs is a cts-time TF and Hz is a discrete-time TF.
Both Hs and Hz are SISO zpk objects.
This function attempts to cancel poles in Hs with zeros in Hz.
"""
function evalTFP(Hs,Hz,f::AbstractArray)
	(szeros, spoles, sk) = _zpkdata(Hs, 1)
	(zzeros, zpoles, zk) = _zpkdata(Hz, 1)

	slim = min(1e-3, max(1e-5, eps(1.0)^(1/(1+length(spoles)))))
	zlim = min(1e-3, max(1e-5, eps(1.0)^(1/(1+length(zzeros)))))

	H = zeros(size(f))
	w = 2*pi*f;	s = j*w;	z=exp.(s)
	for i in 1:length(f)
		wi = w[i];	si = s[i];	zi = z[i]
		cancel = false
		if !isempty(spoles)
			cancel = abs.(si .- spoles) .< slim
		end
		if !any(cancel)
			#wi is far from a pole, so just use the product Hs*Hz
			#H[i] = evalTF(Hs,si) * evalTF(Hz,zi);
			#use evalRPoly directly with z/p/k data, much faster. (David Alldred, Feb 4 2013)
			H[i] = sk * evalRPoly(szeros,si) / evalRPoly(spoles,si) *
				zk * evalRPoly(zzeros,zi) / evalRPoly(zpoles,zi)
		else
			#cancel pole(s) of Hs with corresponding zero(s) of Hz
			cancelz = abs.(zi .- zzeros) .< zlim
			sumcancel = sum(1*cancel); sumcancelz = sum(1*cancelz)
			if sumcancelz > sumcancel
				H[i] = 0.0
			elseif sumcancelz < sumcancel
				H[i] = Inf
			else
				H[i] = evalRPoly(szeros,si,sk) *
					zi^sumcancel * evalRPoly(zzeros[.!cancelz],zi,zk) /
					(evalRPoly(spoles[.!cancel],si,1)*evalRPoly(zpoles,zi,1))
			end
		end
	end
	return H
end
evalTFP(Hs,Hz,f) = evalTFP(Hs,Hz,[f])[1] #Single value

"""`evalMixedTF(tf,f, df=1e-5)`

Compute the mixed transfer function tf at a frequency f.
tf is a struct array with fields Hs and Hz wich represent
a continuous-time/discrete-time tfs which must be multiplied together
and then added up.
"""
function evalMixedTF(tf::CTDTPrefilters, f; df=1e-5)
	H = zeros(size(f))
	for i = 1:length(tf)
		H = H + evalTFP(tf[i].Hs, tf[i].Hz, f)
	end

	err = findall(isnan.(H) .| isinf.(H))
	if !isempty(err)
		#Need to fill in the holes. !! Use a "nearby" frequency point
		for i in err
			H[i] = evalMixedTF(tf, f[i] .+ df, df=df*10)
		end
	end
	return H
end

#==Calculators
===============================================================================#

"""`g = rmsGain(H,f1,f2,N=100)`

Compute the root mean-square gain of the discrete-time
tf H in the frequency band (f1,f2)
"""
function rmsGain(H,f1,f2,N::Int=100)
	w = collect(range(2*pi*f1, stop=2*pi*f2, length=N))
	g = norm( evalTF(H, exp.(j*w)) ) / sqrt(N)
	return g
end

"""`(pGain, Nimp) = powerGain(num,den,Nimp0=100)`

Calculate the power gain of a TF given in coefficient form.

 - `Nimp` is the recommended number of impulse response samples for use
   in future calls and `Nimp0` is the suggested number to use.
"""
function powerGain(num, den, Nimp0::Int=100)
	Nimp = Nimp0
	unstable = false

	sys = _tf(num,den,1)
	(imp,) = _impulse(sys, Nimp)
	if sum(abs.(imp[Nimp-10:Nimp])) < 1e-8 && Nimp > 50 #Use fewer samples.
		Nimp = round(Int, Nimp/1.3)
	else
		while sum(abs.(imp[Nimp-10:Nimp])) > 1e-6
			Nimp = Nimp*2
			(imp,) = _impulse(sys, Nimp)
			if sum(abs.(imp[Nimp-10:Nimp])) >= 50 || Nimp >= 1e4
				#H is close to being unstable
				unstable = true
				break
			end
		end
	end
	pGain = unstable ? Inf : sum(imp .^ 2)
	return (pGain, Nimp)
end

"""`zpk2 = cancelPZ(zpk1, tol=1e-6)`

Cancel zeros/poles in a zpk system
"""
function cancelPZ(zpk1, tol::Float64=1e-6)
	(z, p) = _zpkdata(zpk1, 1)
	zpk2 = deepcopy(zpk1) #Captures k, and whether tf is in s or z

	#Need to go in reverse order because z gets smaller with each cancellation
	for i = length(z):-1:1
		d = z[i] .- p
		cancel = findall(abs.(d) .< tol)
		if !isempty(cancel)
			deleteat!(p, cancel[1])
			deleteat!(z, i)
		end
	end
	zpk2.z = z
	zpk2.p = p
	return zpk2
end

#==infnorm
===============================================================================#
"""`nabsH(w, H)`

Computes the negative of the absolute value of H
at the specified frequency on the unit circle.

This function is used by infnorm().
"""
function nabsH(w, H)
	z = exp.(j*w)
	return -abs.(evalTF(H, z))
end

"""`(Hinf,fmax) = infnorm(H)`

Find the infinity norm of a z-domain transfer function.

Get a rough idea of the location of the maximum.
"""
function infnorm(H)
	N = 129
	w = range(0, stop=2π, length=N); dw = 2π/(N-1)
	Hval = abs.(evalTF(H, exp.(j*w)))
	wi = argmax(abs.(Hval))[1]
	Hinf = Hval[wi]

	local wmax
	if true #NEEDS_OPTIMIZATION
		opt = Optim.Options(x_tol=1e-8, f_tol=1e-6)
		f2min(w) = nabsH(w, H) #Function to minimize (optimum point)
		lower = w[wi]-dw; upper = w[wi]+dw
		wmax = optimize(f2min, lower, upper, GoldenSection()).minimum
	else
		msg = "Hinf: Warning. Optimization toolbox functions not found.\n" *
			" The result returned may not be very accurate."
		@warn(msg)
		wmax = w[wi]
	end

	Hinf = -nabsH(wmax, H)
	fmax = wmax/(2*pi)
	return (Hinf, wmax)
end

#==calculateTF
===============================================================================#
#Use a loophole to set complex poles and zeros in zpk objects
#MALaforge: I think this tries to implement minreal()
function setPolesAndZeros(z,p,k)
	tol = 3e-5 #tolerance for pole-zero cancellation
	z = deepcopy(z) #Don't change original
	Z = _zpk(ComplexF64[0.0],ComplexF64[],1.0,1.0)
	ztf = _zpk(ComplexF64[],ComplexF64[],k,1.0)
	for i in 1:length(p)
		match = abs.(p[i] .- z) .< tol
		if any(match)
			f = findall(match)
			deleteat!(z, f[1])
		else
			#ztf = ztf/(Z-p[i]) #Original code
			push!(ztf.p, p[i]) #VERIFYME: Think this is intent of original code
		end
	end
	for i in 1:length(z)
		#ztf = ztf*(Z-z[i]) #Original code
		push!(ztf.z, z[i]) #VERIFYME: Think this is intent of original code
	end
	return ztf
end

"""`(NTF,STF) = calculateTF(ABCD, k=1, wantSTF=true)`

Calculate the NTF and STF of a delta-sigma modulator whose loop filter
is described by the ABCD matrix, assuming a quantizer gain of k.
The NTF and STF are zpk objects.
"""
function calculateTF(ABCD, k::Vector=[1], wantSTF::Bool=true)
	nq = length(k)
	n = size(ABCD,1) - nq
	nu = size(ABCD,2) - size(ABCD,1)

	A,B,C,D = partitionABCD(ABCD, nu+nq)
	B1 = B[:,1:nu]
	B2 = B[:,nu+1:end]
	D1 = D[:,1:nu]
	if any(any( D[:,nu+1:end] .!= 0 ))
		throw("D2 must be zero")
	end
	K = diagm(k)
	#Find the noise transfer function by forming the closed-loop
	#system (sys_cl) in state-space form.
	Acl = A + B2*K*C
	Bcl = hcat(B1 + B2*K*D1, B2)
	Ccl = K*C
	Dcl = [K*D1 eye(nq)]
	#map(display, [Acl, Bcl, Ccl, Dcl])

	local NTF, STF
	tol = min(1e-3, max(1e-6, eps(1.0)^(1/(size(ABCD,1)))))
	if false # true || all(imag.(ABCD) .== 0) #real modulator
		#REQUESTHELP
		sys_cl = _ss(Acl,Bcl,Ccl,Dcl,1)
		sys_tf = ControlSystems.tf(sys_cl)
		mtfs = _minreal(sys_tf, tol)
		display(mtfs)
		@show (mtfs.nu, mtfs.ny), (nu, nq)
		#MALaforge: Can't make sense of the following; not sure what structure should be:
		STF = mtfs[:,1:nu]
		NTF = mtfs[:,nu+1:nu+nq]
	else #quadrature modulator and less advanced language features
		p = eig(Acl)
		NTFz = eig(A)
		NTF = setPolesAndZeros(NTFz, p, 1)
		if wantSTF
	        @warn("Sorry. calculateTF cannot compute the STF for a complex modulator with this version of the code.")
   	 end
	    STF=nothing
	end
	return (NTF, STF)
end

#Last line


text:
#=
Best possible score of a PWM

# Input
    `pwm::Matrix{Real}`: a 4 x m matrix
# Output
    the best possible score this matrix can get
=#
function best_score(pwm::Matrix{T}) where T <: Real
    sum(maximum(pwm[:,i]) for i = 1:size(pwm,2));
end

function best_score(pwm_col::Vector{T}) where T <: Real
    maximum(pwm_col);
end

#=
Worst possible score of a PWM

# Input
    `pwm::Matrix{Real}`: a 4 x m matrix
# Output
    the worst possible score this matrix can get
=#
function worst_score(pwm::Matrix{T}) where T <: Real
    sum(minimum(pwm[:,i]) for i = 1:size(pwm,2));
end

function worst_score(pwm_col::Vector{T}) where T <: Real
    minimum(pwm_col);
end
#=
Return a column-permuted PWM that minimize the score range so that δ₁ ≥ δ₂ ≥ … ≥ δₘ
where δᵢ = best_score(pwm[:,i])-worst_score(pwm[:,i]).

# Input
    `pwm::Matrix{Real}`: a 4 x m matrix
# Output
    a column-permuted pwm
=#
function min_score_range(pwm::Matrix{T}) where T <: Real
    pwm[:,sortperm([best_score(pwm[:,i])-worst_score(pwm[:,i]) for i = 1:size(pwm,2)],rev=true)];
end

#=
"Round the PWM"

# Input
    `pwm::Matrix{Real}`: a 4 x m matrix
    `granularity`: a small positive real number e.g. 0.01 or 0.001, etc.
# Output
    a rounded pwm of the input pwm
=#
function round_pwm(pwm::Matrix{T}, granularity::Real) where T <: Real
    floor.(pwm ./ granularity) * granularity;
end

#=
The maximum error induced by the rounded pwm M_ϵ
(see definition 3 in https://almob.biomedcentral.com/articles/10.1186/1748-7188-2-15; this is the quantity E)

# Input
    `pwm::Matrix{Real}`: a 4 x m matrix
    `granularity`:
# Output
    A positve real number that's the maximum error induced by the rounded pwm M_ϵ
=#
function calc_E(pwm, pwm_rounded)
    sum(maximum(pwm[:,i]-pwm_rounded[:,i]) for i = 1:size(pwm,2));
end

#=
Note: Use a nested dictionary to represent the distribution Q
    Since Q is used to reference the probability of (M[1…i],score),
    the keys in the first layer is i, and the value in the first layer
    are dictionaries with scores as keys and probability as values

    call create_Q(m) to initialize such a distribution Q
        where m is the "width" of the PWM
=#
create_Q(m::Integer) = Dict{Int16,SortedDict{Double64,Double64}}(i==0 ? i=>SortedDict(0=>1) : i=>SortedDict() for i=0:m);

#=
Input:
    pwm: a 4 x m matrix
    α, β: score interval [α, β]
    bg: 4 x 1 vector that specifies the multinomial genomic background; default to flat background.
Output:
    Q: a probability mass table
        e.g. Q[m] shows all the weights of P[pwm_score = η] for α ≤ η ≤ β
=#
function score_distribution(pwm_::Matrix{T}, α::Real, β::Real, bg=[.25,.25,.25,.25]) where T <: Real
    m = size(pwm_,2);
    Q = create_Q(m);
    @inbounds for i = 1:m
        bs = i+1 > m ? 0 : best_score(pwm_[:,i+1:m]);
        ws = i+1 > m ? 0 : worst_score(pwm_[:,i+1:m]);
        for score in keys(Q[i-1])
            for j = 1:4
                t = score + pwm_[j,i];
                if α - bs ≤ t ≤ β - ws
                    if haskey(Q[i], t)
                        Q[i][t] += Q[i-1][score]*bg[j];
                    else
                        Q[i][t] = Q[i-1][score]*bg[j];
                    end
                end
            end
        end
    end
    return Q
end

#=
Return the probability of the background model that has score ≥ α with respect to the input pwm

Input:
    pwm: a 4 x m matrix
    α: a score threshold
    bg: 4 x 1 vector that specifies the multinomial genomic background; default to flat background.
Output:
    pval: the probability that the background can have score ≥ α with respect to pwm
=#
function fast_pvalue(pwm::Matrix{T}, α::Real, bg=[.25,.25,.25,.25]) where T <: Real
    m = size(pwm,2);
    Q = create_Q(m);
    pval = 0f0;
    @inbounds for i = 1:m
        bs = i+1 > m ? 0 : best_score(pwm[:,i+1:m]);
        ws = i+1 > m ? 0 : worst_score(pwm[:,i+1:m]);
        for (score,_) in Q[i-1]
            for j = 1:4
                t = score + pwm[j,i];
                if α - ws ≤ t
                    pval = pval + Q[i-1][score]*bg[j];
                elseif α - bs ≤ t
                    if haskey(Q[i], t)
                        Q[i][t] += Q[i-1][score]*bg[j];
                    else
                        Q[i][t] = Q[i-1][score]*bg[j];
                    end
                end
            end
        end
    end
    return pval
end

# return the sum of all the weights
Q_sum(Q_m::SortedDict{Double64,Double64}) = sum(v for v in values(Q_m));


text:
using GaussianDistributions, Random
import Random.rand
import Distributions.logpdf
import GaussianDistributions: whiten, unwhiten

"""
    StartingPtPrior

Types inheriting from the abstract type `StartingPtPrior` indicate the prior
that is put on the starting point of the observed path following the dynamics
of some stochastic differential equation.
"""
abstract type StartingPtPrior{T} end

"""
    KnownStartingPt{T} <: StartingPtPrior

Indicates that the starting point is known and stores its value in `y`
"""
struct KnownStartingPt{T} <: StartingPtPrior{T}
    y::T
    KnownStartingPt(y::T) where T = new{T}(y)
end

"""
    GsnStartingPt{T,S} <: StartingPtPrior

Indicates that the starting point is equipped with a Gaussian prior with
mean `μ` and covariance matrix `Σ`. It also stores the most recently sampled
white noise `z` used to compute the starting point and a precision matrix
`Λ`:=`Σ`⁻¹. `μ₀` and `Σ₀` are the mean and covariance of the white noise

    GsnStartingPt(y::T, μ::T, Σ::S)

Base constructor explicitly initialising starting point `y` (on top of mean `μ`
and covariance matrix `Σ`) which is internally transformed to its corresponding
white noise `z`

    GsnStartingPt(μ::T, Σ::S)

Base constructor v2. It initialises the mean `μ` and covariance `Σ` parameters
and samples white noise `z`. `Λ` is set according to `Λ`:=`Σ`⁻¹

    GsnStartingPt(z::T, G::GsnStartingPt{T,S})

Copy constructor. It copies the mean `μ`, covariance `Σ`, precision `Λ` as well
as mean of the white noise `μ₀` and covariance of the white noise `Σ₀` from
the old prior `G` and sets a new white noise to `z`.
"""
struct GsnStartingPt{T,S} <: StartingPtPrior{T} where {S}
    z::T
    μ::T
    Σ::S
    Λ::S
    μ₀::T
    Σ₀::UniformScaling

    function GsnStartingPt(y::T, μ::T, Σ::S) where {T,S}
        z = whiten(Σ, y-μ)
        new{T,S}(z, μ, Σ, inv(Σ), 0*μ, I)
    end

    function GsnStartingPt(μ::T, Σ::S) where {T,S}
        μ₀ = 0*μ
        Σ₀ = I
        z = rand(Gaussian(μ₀, Σ₀))
        new{T,S}(z, μ, Σ, inv(Σ), μ₀, Σ₀)
    end

    function GsnStartingPt(z::T, G::GsnStartingPt{T,S}) where {T,S}
        new{T,S}(z, G.μ, G.Σ, G.Λ, G.μ₀, G.Σ₀)
    end
end

"""
    rand(G::GsnStartingPt, ρ)

Sample new white noise using Crank-Nicolson scheme with memory parameter `ρ` and
a previous value of the white noise stored inside object `G`
"""
function rand(G::GsnStartingPt, ρ=0.0)
    zᵒ = rand(Gaussian(G.μ₀, G.Σ₀))
    z = √(1-ρ)*zᵒ + √(ρ)*G.z # preconditioned Crank-Nicolson
    GsnStartingPt(z, G)
end

"""
    rand(G::KnownStartingPt, ::Any)

If starting point is known then nothing can be sampled
"""
rand(G::KnownStartingPt, ::Any=nothing) = G

"""
    startPt(G::GsnStartingPt, P)

Compute a new starting point from the white noise for a given posterior
distribution obtained from combining prior `G` and the likelihood encoded by the
object `P`.
"""
function startPt(G::GsnStartingPt, P::GuidPropBridge)
    μₚₒₛₜ = (P.H[1] + G.Λ) \ (P.Hν[1] + G.Λ * G.μ)
    Σₚₒₛₜ = inv(P.H[1] + G.Λ)
    Σₚₒₛₜ = 0.5 * (Σₚₒₛₜ + Σₚₒₛₜ') # remove numerical inaccuracies
    unwhiten(Σₚₒₛₜ, G.z) + μₚₒₛₜ
end

"""
    startPt(G::GsnStartingPt, P)

Compute a new starting point from the white noise for a given prior
distribution `G`
"""
startPt(G::GsnStartingPt) = unwhiten(G.Σ, G.z) + G.μ

"""
    startPt(G::KnownStartingPt, P)

Return starting point
"""
startPt(G::KnownStartingPt, P::GuidPropBridge) = G.y

"""
    startPt(G::KnownStartingPt, P)

Return starting point
"""
startPt(G::KnownStartingPt) = G.y

"""
    invStartPt(y, G::GsnStartingPt, P::GuidPropBridge)

Compute the driving noise that is needed to obtain starting point `y` under
prior `G` and likelihood `P`. Return a new starting point object
"""
function invStartPt(y, G::GsnStartingPt, P::GuidPropBridge)
    μₚₒₛₜ = (P.H[1] + G.Λ) \ (P.Hν[1] + G.Λ * G.μ)
    Σₚₒₛₜ = inv(P.H[1] + G.Λ)
    Σₚₒₛₜ = 0.5 * (Σₚₒₛₜ + Σₚₒₛₜ')
    GsnStartingPt(whiten(Σₚₒₛₜ, y-μₚₒₛₜ), G)
end

"""
    invStartPt(y, G::KnownStartingPt, P::GuidPropBridge)

Starting point known, no need for dealing with white noise
"""
invStartPt(y, G::KnownStartingPt, P::GuidPropBridge) = G

"""
    logpdf(::KnownStartingPt, y)

Nothing to do so long as `y` is equal to the known starting point inside `G`
"""
function logpdf(G::KnownStartingPt, y)
    (G.y == y) ? 0.0 : error("Starting point is known, but a different value ",
                             "was passed to logpdf.")
end

"""
    logpdf(::GsnStartingPt, y)

log-probability density function evaluated at `y` of a prior distribution `G`
"""
logpdf(G::GsnStartingPt, y) = logpdf(Gaussian(G.μ, G.Σ), y)


text:
using
  PassiveTracerFlows,
  Test,
  Statistics,
  Random

import # use 'import' rather than 'using' for submodules to keep namespace clean
  PassiveTracerFlows.TracerAdvectionDiffusion

# the devices on which tests will run
devices = CUDA.has_cuda() ? (CPU(), GPU()) : (CPU(),)

const rtol_traceradvectiondiffusion = 1e-12 # tolerance for rtol_traceradvdiff tests

# Run tests
testtime = @elapsed begin

for dev in devices

  println("testing on "*string(typeof(dev)))

  @testset "TracerAdvectionDiffusion" begin
    include("test_traceradvectiondiffusion.jl")

    stepper = "RK4"
    dt, nsteps  = 1e-2, 40
    @test test_constvel(stepper, dt, nsteps, dev)
    dt, tfinal  = 0.002, 0.1
    @test test_timedependentvel(stepper, dt, tfinal, dev)
    dt, tfinal  = 0.005, 0.1
    @test test_diffusion(stepper, dt, tfinal, dev; steadyflow=true)
    dt, tfinal  = 0.005, 0.1
    @test test_diffusion(stepper, dt, tfinal, dev; steadyflow=false)
    dt, tfinal  = 0.005, 0.1
    @test test_hyperdiffusion(stepper, dt, tfinal, dev)

    @test TracerAdvectionDiffusion.noflow(π) == 0
  end

end

end #time

println("Total test time: ", testtime)

max_stars_repo_name: FourierFlows/PassiveTracerFlows

max_issues_repo_name: FourierFlows/PassiveTracerFlows

max_forks_repo_name: FourierFlows/PassiveTracerFlows


text:
using ApproxFunOrthogonalPolynomials, StaticArrays, Test

@testset "Clenshaw" begin
    test_function = Fun(sin, 50)
    @test test_function(0.9) ≈ clenshaw(0.9, SVector{50}(test_function.coefficients))
end