text:
using Test
using Pluto
import Pluto.ExploreExpression: SymbolsState, compute_symbolreferences

function testee(expr, ref, def, funccalls, funcdefs, verbose=true)
    funcdefs_dict = map(funcdefs) do (k, v)
        k => if v isa SymbolsState
            v
        else
            SymbolsState(Set(v[1]), Set(v[2]), Set(v[3]))
        end
    end |> Dict
    expected = SymbolsState(Set(ref), Set(def), Set(funccalls), funcdefs_dict)
    result = compute_symbolreferences(expr)

    # Anonymous function are given a random name, which looks like anon67387237861123
    # To make testing easier, we rename all such functions to anon
    new_name(sym) = startswith(string(sym), "anon") ? :anon : sym

    result.assignments = Set(new_name.(result.assignments))
    result.funcdefs = let
        newfuncdefs = Dict()
        for (k,v) in result.funcdefs
            newfuncdefs[new_name(k)] = v
        end
        newfuncdefs
    end

    if verbose && expected != result
        println()
        println("FAILED TEST")
        println(expr)
        println()
        dump(expr, maxdepth = 20)
        println()
        @show expected
        @show result
        println()
    end
    return expected == result
end

# nowarn tests are for functionality that is not yet implemented
# (but we don't want any expressions to error/warn)
# Once the functionality has been implemented, they should be changed to normal tests

@testset "Explore Expressions" begin
    @testset "Basics" begin
        @test testee(:(a), [:a], [], [], [])
        @test testee(:(1 + 1), [:+], [], [:+], [])
        @test testee(:(x = 3), [], [:x], [], [])
        @test testee(:(x = x), [:x], [:x], [], [])
        @test testee(:(x = 1 + y), [:+, :y], [:x], [:+], [])
        @test testee(:(x = +(a...)), [:+, :a], [:x], [:+], [])

        @test_nowarn testee(:(x::Int64 = 3), [], [:x, :Int64], [], [], false)
    end
    @testset "Lists and structs" begin
        @test testee(:(1:3), [:(:)], [], [:(:)], [])
        @test testee(:(a[1:3,4]), [:a, :(:)], [], [:(:)], [])
        @test testee(:(a[1:3,4] = b[5]), [:b], [], [], [])
        @test testee(:(a.property), [:a], [], [], [])
        @test testee(:(a.property = 1), [], [], [], [])
        @test testee(:(struct a; c; d; end), [], [:a], [], [
            :a => ([], [], [])
            ])

        @test_nowarn testee(:(struct a <: b; c; d::Int64; end), [:b, :Int64], [:a], [], [], false)
        @test testee(:(module a; f(x) = x; z = r end), [], [:a], [], [])

    end
    @testset "Modifiers" begin
        @test testee(:(a = a + 1), [:a, :(+)], [:a], [:+], [])
        @test testee(:(a += 1), [:a, :(+)], [:a], [:+], [])
        @test testee(:(a[1] += 1), [:a, :(+)], [], [:+], [])
        @test testee(:(x = let a = 1; a += b end), [:(+), :b], [:x], [:+], [])
    end
    @testset "`for` & `while`" begin
        @test testee(:(for k in 1:n; k + s; end), [:n, :s, :+, :(:)], [], [:+, :(:)], [])
        @test testee(:(for k in 1:2, r in 3:4; global z = k + r; end), [:+, :(:)], [:z], [:+, :(:)], [])
        @test testee(:(while k < 2; r = w; global z = k + r; end), [:k, :(<), :w, :+], [:z], [:+, :(<)], [])
    end
    @testset "Comprehensions" begin
        @test testee(:([sqrt(s) for s in 1:n]), [:sqrt, :n, :(:)], [], [:sqrt, :(:)], [])
        @test testee(:([s + j + r + m for s in 1:3 for j in 4:5 for (r, l) in [(1, 2)]]), [:+, :m, :(:)], [], [:+, :(:)], [])

        @test_nowarn testee(:([a for a in a]), [:a], [], [], [], false)
        @test_nowarn testee(:(a = [a for a in a]), [:a], [:a], [], [], false)
    end
    @testset "Multiple expressions" begin
        @test testee(:(x = let r = 1; r + r end), [:+], [:x], [:+], [])
        @test testee(:(begin let r = 1; r + r end; r = 2 end), [:+], [:r], [:+], [])
        @test testee(:(a, b = 1, 2), [], [:a, :b], [], [])
        @test testee(:((a, b) = 1, 2), [], [:a, :b], [], [])
        @test testee(:((a[1], b.r) = (1, 2)), [], [], [], [])
        @test testee(:((k = 2; 123)), [], [:k], [], [])
        @test testee(:((a = 1; b = a + 1)), [:+], [:a, :b], [:+], [])
        @test testee(:((a = b = 1)), [], [:a, :b], [], [])
        @test testee(:(let k = 2; 123 end), [], [], [], [])
        @test testee(:(let k() = 2 end), [], [], [], [])

        @test_nowarn testee(:(a::Int64, b::String = 1, "2"), [:Int64, :String], [:a, :b], [], [], false)
    end
    @testset "Functions" begin
        @test testee(:(function g() r = 2; r end), [], [:g], [], [
            :g => ([], [], [])
        ])
        @test testee(:(function f(x, y = 1; r, s = 3 + 3) r + s + x * y * z end), [], [:f], [], [
            :f => ([:z, :+, :*], [], [:+, :*])
        ])
        @test testee(:(function f(x) x * y * z end), [], [:f], [], [
            :f => ([:y, :z, :*], [], [:*])
        ])
        @test testee(:(function f(x) x = x / 3; x end), [], [:f], [], [
            :f => ([:/], [], [:/])
        ])
        @test testee(:(function f(x) a end; function f(x, y) b end), [], [:f], [], [
            :f => ([:a, :b], [], [])
        ])
        @test testee(:(f(x, y = a + 1) = x * y * z), [], [:f], [], [
            :f => ([:*, :z], [], [:*])
        ])
        @test testee(:(minimum(x) do (a, b); a + b end), [:x, :minimum], [:anon], [:minimum], [
            :anon => ([:(+)], [], [:+])
        ])
        @test testee(:(f = x->x * y), [], [:f, :anon], [], [
            :anon => ([:y, :*], [], [:*])
        ])
        @test testee(:(f = (x, y)->x * y), [], [:f, :anon], [], [
            :anon => ([:*], [], [:*])
        ])
        @test testee(:(f = (x, y = a + 1)->x * y), [], [:f, :anon], [], [
            :anon => ([:*], [], [:*])
        ])
        @test testee(:((((a, b), c), (d, e))->a * b * c * d * e * f), [], [:anon], [], [
            :anon => ([:*, :f], [], [:*])
        ])

        @test testee(:(func(b)), [:func, :b], [], [:func], [])
        @test testee(:(funcs[i](b)), [:funcs, :i, :b], [], [], [])

        @test_nowarn testee(:(function f(y::Int64 = a)::String string(y) end), [], [], [], [], false)
        @test_nowarn testee(:(function f(x::T; k = 1) where T return x + 1 end), [], [], [], [], false)
        @test_nowarn testee(:(function f(::MIME"text/html") 1 end), [], [], [], [], false)
    end
    @testset "Scope modifiers" begin
        @test testee(:(let global a, b = 1, 2 end), [], [:a, :b], [], [])
        @test testee(:(let global k = 3 end), [], [:k], [], [])
        @test testee(:(let global k += 3 end), [:+, :k], [:k], [:+], [])
        @test testee(:(let global k; k = 4 end), [], [:k], [], [])
        @test testee(:(let global k; b = 5 end), [], [], [], [])

        @test testee(:(begin local a, b = 1, 2 end), [], [], [], [])
        @test testee(:(begin local k = 3 end), [], [], [], [])
        @test testee(:(begin local k += 3 end), [:+], [], [:+], [])
        @test testee(:(begin local k; k = 4 end), [], [], [], [])
        @test testee(:(begin local k; b = 5 end), [], [:b], [], [])

        @test testee(:(function f(x) global k = x end), [], [:f], [], [
            :f => ([], [:k], [])
        ])
        @test testee(:((begin x = 1 end, y)), [:y], [:x], [], [])
        @test testee(:(x = let global a += 1 end), [:(+), :a], [:x, :a], [:+], [])
    end
    @testset "`import` & `using`" begin
        @test testee(:(using Plots), [], [:Plots], [], [])
        @test testee(:(using JSON, UUIDs), [], [:JSON, :UUIDs], [], [])
        @test testee(:(import Pluto), [], [:Pluto], [], [])
        @test testee(:(import Pluto: wow, wowie), [], [:wow, :wowie], [], [])
        @test testee(:(import .Pluto: wow), [], [:wow], [], [])
        @test testee(:(import ..Pluto: wow), [], [:wow], [], [])
    end
    @testset "Macros" begin
        @test testee(:(@time a = 2), [Symbol("@time")], [:a], [], [])
        @test testee(:(html"a $(b = c)"), [Symbol("@html_str")], [], [], [])
        @test testee(:(md"a $(b = c)"), [Symbol("@md_str"), :c], [:b], [], [])
        @test testee(:(md"a \$(b = c)"), [Symbol("@md_str")], [], [], [])
    end
    @testset "String interpolation" begin
        @test testee(:("a $b"), [:b], [], [], [])
        @test testee(:("a $(b = c)"), [:c], [:b], [], [])
    end
end


text:
##################
# Differentiation
##################

function unsafe_eval_element_derivative(b::ChebyshevT{T}, idx::PolynomialDegree, x, order) where {T}
	@assert order >= 0
	if order == 0
		unsafe_eval_element(b, idx, x)
	elseif order == 1
	    d = degree(idx)
	    if d == 0
	        zero(T)
	    else
	        d * unsafe_eval_element(ChebyshevU{T}(length(b)), idx-1, x)
	    end
	else
		error("Higher order derivatives of Chebyshev polynomials not implemented.")
	end
end

# Chebyshev differentiation is so common that we make it its own type
struct ChebyshevDifferentiation{T} <: DictionaryOperator{T}
	src		::	Dictionary
	dest	::	Dictionary
    order   ::  Int
end

ChebyshevDifferentiation(src::Dictionary, dest::Dictionary, order::Int = 1) =
	ChebyshevDifferentiation{operatoreltype(src,dest)}(src, dest, order)

ChebyshevDifferentiation(src::Dictionary, order::Int = 1) =
    ChebyshevDifferentiation(src, src, order)

order(op::ChebyshevDifferentiation) = op.order

string(op::ChebyshevDifferentiation) = "Chebyshev differentiation matrix of order $(order(op)) and size $(length(src(op)))"

similar_operator(op::ChebyshevDifferentiation, src, dest) =
    ChebyshevDifferentiation(src, dest, order(op))

unsafe_wrap_operator(src, dest, op::ChebyshevDifferentiation) = similar_operator(op, src, dest)

function adjoint(op::ChebyshevDifferentiation)
    @warn "Inefficient adjoint of `ChebyshevDifferentiation`"
    ArrayOperator(adjoint(Matrix(op)), dest(op), src(op))
end

conj(op::ChebyshevDifferentiation) = op

# TODO: this allocates lots of memory...
function apply!(op::ChebyshevDifferentiation, coef_dest, coef_src)
    #	@assert period(dest)==period(src)
    n = length(coef_src)
    T = eltype(coef_src)
    tempc = coef_src[:]
    tempr = coef_src[:]
    for o = 1:order(op)
        tempr = zeros(T,n)
        # 'even' summation
        s = 0
        for i=(n-1):-2:2
            s = s+2*i*tempc[i+1]
            tempr[i] = s
        end
        # 'uneven' summation
        s = 0
        for i=(n-2):-2:2
            s = s+2*i*tempc[i+1]
            tempr[i] = s
        end
        # first row
        s = 0
        for i=2:2:n
            s = s+(i-1)*tempc[i]
        end
        tempr[1]=s
        tempc = tempr
    end
    coef_dest[1:n-order(op)] = tempr[1:n-order(op)]
    coef_dest
end

struct ChebyshevAntidifferentiation{T} <: DictionaryOperator{T}
	src		::	Dictionary
	dest	::	Dictionary
    order   ::  Int
end

ChebyshevAntidifferentiation(src::Dictionary, dest::Dictionary, order::Int = 1) =
	ChebyshevAntidifferentiation{operatoreltype(src,dest)}(src, dest, order)

ChebyshevAntidifferentiation(src::Dictionary, order::Int = 1) =
    ChebyshevAntidifferentiation(src, src, order)

order(op::ChebyshevAntidifferentiation) = op.order

string(op::ChebyshevAntidifferentiation) = "Chebyshev antidifferentiation matrix of order $(order(op)) and size $(length(src(op)))"

similar_operator(op::ChebyshevAntidifferentiation, src, dest) =
    ChebyshevAntidifferentiation(src, dest, order(op))

unsafe_wrap_operator(src, dest, op::ChebyshevAntidifferentiation) = similar_operator(op, src, dest)

# TODO: this allocates lots of memory...
function apply!(op::ChebyshevAntidifferentiation, coef_dest, coef_src)
    #	@assert period(dest)==period(src)
    T = eltype(coef_src)
    tempc = zeros(T,length(coef_dest))
    tempc[1:length(coef_src)] = coef_src[1:length(coef_src)]
    tempr = zeros(T,length(coef_dest))
    tempr[1:length(coef_src)] = coef_src[1:length(coef_src)]
    for o = 1:order(op)
        n = length(coef_src)+o
        tempr = zeros(T,n)
        tempr[2]+=tempc[1]
        tempr[3]=tempc[2]/4
        tempr[1]+=tempc[2]/4
        for i = 3:n-1
            tempr[i-1]-=tempc[i]/(2*(i-2))
            tempr[i+1]+=tempc[i]/(2*(i))
            tempr[1]+=real(-1im^i)*tempc[i]*(2*i-2)/(2*i*(i-2))
        end
        tempc = tempr
    end
    coef_dest[:]=tempr[:]
    coef_dest
end

function differentiation(::Type{T}, src::ChebyshevT, dest::ChebyshevT, order::Int; options...) where {T}
	@assert order >= 0
	if order == 0
		IdentityOperator{T}(src, dest)
	else
		ChebyshevDifferentiation{T}(src, dest, order)
	end
end

function antidifferentiation(::Type{T}, src::ChebyshevT, dest::ChebyshevT, order::Int; options...) where {T}
	@assert order >= 0
	if order == 0
		IdentityOperator{T}(src, dest)
	else
    	ChebyshevAntidifferentiation{T}(src, dest, order)
	end
end

ext: jl


text:
using DataFrames, JSON
using StringDistances, Statistics, DataStructures, Random, Pipe
import ParquetFiles
import Base: isless

include("utils.jl")
include("matchtypes.jl")
include("constructors.jl")
include("steps/init.jl")
include("steps/extend.jl")
include("steps/match.jl")
include("steps/integrate.jl")
include("steps/common.jl")

include("../test/testutils.jl")

## GEN

esp_dic = Dict("fing"=> generate_dfpair(3,1,0))
relp = RelPair([[1 1; 1 2;2 2], [2 2; 2 1]], ("fing", "fing"))
crs = CorefSystem(esp_dic, [relp])

shortside = 10
longside = 10
max_depth = 3
preplen = 3

crr = CorefResolver(crs, shortside, longside, max_depth, preplen);

steps = [
    InitParams("fing"),
    ExtendParams(1, 1),
    MatchParams(2),
    IntegrateParams(),
]

step!(crr, steps[1])

crr.leftspace_main
crr.space_pairs[1]
crr.space_pairs[2]

for erspace in crr |> main |> spaces
    (erspace.width == 1) |> println
    (erspace.height == 3) |> println
    (erspace.space[1:3,1] == [1, 2, 3]) |> println
end

step!(crr, steps[2])

(crr |> main).erspace1.width == 2
(crr |> main).erspace1.space[1:3, 1:2] == [1 1; 1 2; 2 2]

(crr |> main).erspace1.space[1:7,:]

step!(crr, steps[3])

crr.leftspace_main
crr.space_pairs[1]
crr.space_pairs[2]

(crr |> main).erspace1.space
(crr |> main).erspace2.space

step!(crr, steps[4])

crr.crspace.latent_match_vars

### FILE:
crsysarr = [CorefSystem("./test/test_data/test_sys_$i") for i in 0:2]

crs = crsysarr[1]

shortside = 600
longside = 200
max_depth = 3
preplen = 5

crr = CorefResolver(crs, shortside, longside, max_depth, preplen);

steps = [
    InitParams("e"),
    MatchParams(8),
    IntegrateParams(),
]

for params in steps
    step!(crr, params)
end

step!(crr, steps[end])

(crr |> main).erspace1.width#space[1:8,:]

res = crr |> result

tres0 = load_result("./test/test_data/test_sys_0")

res["e"] == tres0["e"]


text:
using DashBootstrapComponents, DashHtmlComponents

toast = html_div([
    dbc_button(
        "Open toast",
        id = "positioned-toast-toggle",
        color = "primary",
        n_clicks = 0,
    ),
    dbc_toast(
        "This toast is placed in the top right",
        id = "positioned-toast",
        header = "Positioned toast",
        is_open = false,
        dismissable = true,
        icon = "danger",
        # top: 66 positions the toast below the navbar
        style = Dict("position" => "fixed", "top" => 66, "right" => 10, "width" => 350),
    ),
]);

callback!(
    app,
    Output("positioned-toast", "is_open"),
    Input("positioned-toast-toggle", "n_clicks"),
) do n
    return n > 0
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using MITgcmTools, MeshArrays
using Test

@testset "MITgcmTools.jl" begin

    #format conversions
    (γ,Γ)=GridOfOnes("CubeSphere",30,30)
    @test isa(convert2gcmfaces(Γ["XC"]),Array)
    @test isa(convert2array(Γ["XC"]),Array)

    #physical oceanography
    (ρP,ρI,ρR) = SeaWaterDensity(3.,35.5,3000.)
    @test isapprox(ρI,1041.83267, rtol=1e-6)

    D=collect(0.0:1.0:500.0); tmp=(1.0.-tanh.(5*(-1 .+ 2/D[end]*D)));
    T=2.0 .+ 8.0*tmp; S=34.0 .+ 0.5*tmp;
    (ρP,ρI,ρR) = SeaWaterDensity(T,S,D);
    mld=MixedLayerDepth(T,S,D,"BM");
    @test isapprox(mld,134.0)

end


text:
# this file is only for julia older than 1.3

using Libdl

if "FLINT_PATH" in keys(ENV)

   if VERSION >= v"1.3.0-rc4"
      error("please use Pkg Overrides.toml for custom flint installations in julia >= 1.3")
   end

   flint_path = ENV["FLINT_PATH"]
   if !isfile(flint_path)
      error("FLINT_PATH: file does not exists")
   end

   open(joinpath(@__DIR__,"deps.jl"), "w") do f
      println(f, """
using Libdl
libflint = "$flint_path"
if Libdl.dlopen_e(libflint) in (C_NULL, nothing)
    error("$(libflint) cannot be opened, Please check FLINT_PATH environment variable and re-run Pkg.build("LoadFlint"), then restart Julia.")
end
# there must be some libgmp loaded for libflint
libgmp = filter(x->occursin(r"libgmp[.-]", x), dllist())[1]
libmpfr = filter(x->occursin(r"libmpfr[.-]", x), dllist())[1]

""")
   end

elseif VERSION < v"1.3.0-rc4"

   using Pkg, BinaryProvider

   # This does not work on julia >= 1.3, but there we use the *jll package anyway.
   ver = Pkg.API.__installed(PKGMODE_MANIFEST)["FLINT_jll"]

   # Parse some basic command-line arguments
   const verbose = "--verbose" in ARGS

   # GMP and MPFR might not be needed on unix as julia should have those loaded already
   # but on windows flint will not load without them so we leave them here for simplicity
   dependencies = [
     # This has to be in sync with the jll packages (using generate_build.jl and build_tarballs.jl from Yggdrasil)
     "build_GMP.v6.1.2.jl",
     "build_MPFR.v4.0.2.jl",
   ]

   if ver == v"200.700.0+0"
     push!(dependencies, "build_FLINT.v200.700.0.jl")
   elseif ver == v"200.700.100+0"
     push!(dependencies, "build_FLINT.v200.700.100.0.jl")
   elseif ver == v"200.800.0+1"
     push!(dependencies, "build_FLINT.v200.800.0.jl")
   elseif ver == v"200.800.100+0"
     push!(dependencies, "build_FLINT.v200.800.100.jl")
   elseif ver == v"200.800.300+0"
     push!(dependencies, "build_FLINT.v200.800.300.jl")
   elseif ver == v"200.800.400+0"
     push!(dependencies, "build_FLINT.v200.800.400.jl")
   else
     throw(error("Flint version $ver not supported for julia version <= 1.3"))
   end

   const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

   products = []

   # Execute the build scripts for the dependencies in an isolated module to avoid overwriting
   # any variables/constants here
   for file in dependencies
       build_file = joinpath(@__DIR__, file)
       m = @eval module $(gensym()); include($build_file); end
       append!(products, m.products)
   end

   write_deps_file(joinpath(@__DIR__, "deps.jl"), Array{Product,1}(products), verbose=verbose)

end


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

abstract type SolverState end
abstract type SolverMetric end
abstract type Solver end

struct Optimizer
    method::Any
    initial_state::Any
    update_state!::Function
    check_state::Function
    check_convergence::Function
end

OptimizationTrace{T<:SolverMetric} = Vector{T}

mutable struct DeconvolutionOptions
    max_iters::Int
    show_trace::Bool
    check_every::Int
    time_limit::Float64
    x0::AbstractArray
end

function DeconvolutionOptions(;
    max_iters::Integer = 100,
    show_trace::Bool = true,
    check_every::Integer = 10,
    time_limit = NaN,
    x0::AbstractArray = [],
)
    check_every = check_every > 0 ? check_every : 1
    DeconvolutionOptions(Int(max_iters), show_trace, check_every, time_limit, x0)
end

mutable struct OptimizationResult
    x::AbstractArray
    trace::OptimizationTrace
    converged::Bool
end


text:
export ignore

import Base: show

"""
    ignore(count::Int)

Creates a `ignore` operator, which returns an Observable
that skips the first `count` items emitted by the source Observable.

# Arguments
- `count::Int`: the number of times, items emitted by source Observable should be skipped.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from(1:5)

subscribe!(source |> ignore(2), logger())
;

# output

[LogActor] Data: 3
[LogActor] Data: 4
[LogActor] Data: 5
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
ignore(count::Int) = IgnoreOperator(count)

struct IgnoreOperator <: InferableOperator
    count :: Int
end

function on_call!(::Type{L}, ::Type{L}, operator::IgnoreOperator, source) where L
    return proxy(L, source, IgnoreProxy(operator.count))
end

operator_right(operator::IgnoreOperator, ::Type{L}) where L = L

struct IgnoreProxy <: ActorProxy
    count :: Int
end

actor_proxy!(::Type{L}, proxy::IgnoreProxy, actor::A)  where { L, A } = IgnoreActor{L, A}(proxy.count, actor)

mutable struct IgnoreActor{L, A} <: Actor{L}
    count         :: Int
    actor         :: A
    skipped_count :: Int

    IgnoreActor{L, A}(count::Int, actor::A) where { L, A } = new(count, actor, 0)
end

function on_next!(actor::IgnoreActor{L}, data::L) where L
    if actor.skipped_count < actor.count
        actor.skipped_count += 1
    else
        next!(actor.actor, data)
    end
end

function on_error!(actor::IgnoreActor, err)
    error!(actor.actor, err)
end

function on_complete!(actor::IgnoreActor)
    complete!(actor.actor)
end

Base.show(io::IO, ::IgnoreOperator)          = print(io, "IgnoreOperator()")
Base.show(io::IO, ::IgnoreProxy)             = print(io, "IgnoreProxy()")
Base.show(io::IO, ::IgnoreActor{L})  where L = print(io, "IgnoreActor($L)")


text:
# Bayesian logistic regression with independent Gaussian priors
using StatsFuns: logistic

type SteinLogisticRegressionGaussianPrior <: SteinLogisticRegressionPosterior
    X::Array{Float64, 2}
    y::Array{Float64, 1} # we assume the y = -1, 1
    priormu::Array{Float64, 1}
    priorvar::Array{Float64, 1}  # we assume components are independent
    c1::Float64
    c2::Float64
    c3::Float64
end

SteinLogisticRegressionGaussianPrior(X::Array{Float64, 2}, y::Array{Float64, 1}) =
SteinLogisticRegressionGaussianPrior(X, y,
                                     repmat([0.0], size(X, 2)),
                                     repmat([1.0], size(X, 2)))

# Constructor with canonical setting of (c1,c2,c3) = (1,1,1)
SteinLogisticRegressionGaussianPrior(X::Array{Float64, 2},
                                     y::Array{Float64, 1},
                                     priormu::Array{Float64, 1},
                                     priorvar::Array{Float64, 1}) =
SteinLogisticRegressionGaussianPrior(X, y, priormu, priorvar,
                                     1.0, 1.0, 1.0)

function logprior(d::SteinLogisticRegressionGaussianPrior,
                  beta::Array{Float64, 1})
    # Independent Gaussian prior for each component
    logpriorvalue =
        -0.5 * sum((beta .- d.priormu') .^ 2 ./ d.priorvar') -
        -0.5 * sum(log(2 * pi .* d.priorvar))
    logpriorvalue
end

function gradlogprior(d::SteinLogisticRegressionGaussianPrior,
                      beta::Array{Float64, 1})
    -(beta .- d.priormu) ./ d.priorvar
end

function numdimensions(d::SteinLogisticRegressionGaussianPrior)
    length(d.priormu)
end

function priordiaghessian(d::SteinLogisticRegressionGaussianPrior,
                          beta::Array{Float64, 1})
    -1.0 ./ d.priorvar
end


text:
seed!(42)
L = 3
D = 10
N = 20
nSamples = 10
b = 5
x = rand(N, D)
y = rand(N)
@testset "MCIntegrationVI" begin
    C = 10.0
    i = AGP.tuple_inference(MCIntegrationVI(clipping = C), 1, D, N, N, [], [])
    l = GaussianLikelihood()
    # @test AGP.grad_quad(l, 200.0, 0.0, 1.0, i) == (C, -C)
end


text:
"""
Solve the least squares problem using a SVD of the design matrix.
"""
struct MMSVD <: AbstractMMAlg end

# Initialize data structures.
function __mm_init__(::MMSVD, problem::MVDAProblem, ::Nothing)
    @unpack coeff = problem
    X = get_design_matrix(problem)
    n, p, c = probdims(problem)
    T = floattype(problem)
    nparams = ifelse(problem.kernel isa Nothing, p, n)

    # thin SVD of X
    U, s, V = __svd_wrapper__(X)
    r = length(s) # number of nonzero singular values

    # worker arrays
    Z = Matrix{T}(undef, n, c-1)
    buffer = Matrix{T}(undef, r, c-1)

    # diagonal matrices
    Ψ = Diagonal(Vector{T}(undef, r))

    return (;
        projection=StructuredL0Projection(nparams),
        U=U, s=s, V=V,
        Z=Z, Ψ=Ψ,
        buffer=buffer,
    )
end

# Check for data structure allocations; otherwise initialize.
function __mm_init__(::MMSVD, problem::MVDAProblem, extras)
    if :projection in keys(extras) && :buffer in keys(extras) # TODO
        return extras
    else
        __mm_init__(MMSVD(), problem, nothing)
    end
end

# Update data structures due to change in model size, k.
__mm_update_sparsity__(::MMSVD, problem::MVDAProblem, ϵ, ρ, k, extras) = nothing

# Update data structures due to changing ρ.
__mm_update_rho__(::MMSVD, problem::MVDAProblem, ϵ, ρ, k, extras) = update_diagonal(problem, ρ, extras)

# Update data structures due to changing λ.
__mm_update_lambda__(::MMSVD, problem::MVDAProblem, ϵ, λ, extras) = update_diagonal(problem, λ, extras)

function update_diagonal(problem::MVDAProblem, λ, extras)
    @unpack s, Ψ = extras
    n, _, _ = probdims(problem)
    a², b² = 1/n, λ

    # Update the diagonal matrices Ψ = (a² Σ²) / (a² Σ² + b² I).
    for i in eachindex(Ψ.diag)
        sᵢ² = s[i]^2
        Ψ.diag[i] = a² * sᵢ² / (a² * sᵢ² + b²)
    end

    return nothing
end

# Apply one update.
function __mm_iterate__(::MMSVD, problem::MVDAProblem, ϵ, ρ, k, extras)
    @unpack intercept, coeff, proj = problem
    @unpack buffer, projection = extras
    @unpack Z, Ψ, U, s, V = extras
    Σ = Diagonal(s)
    T = floattype(problem)

    # need to compute Z via residuals...
    apply_projection(projection, problem, k)
    __evaluate_residuals__(problem, ϵ, extras, true, false, true)

    # Update parameters:B = P + V * Ψ * (Σ⁻¹UᵀZ - VᵀP)
    B = coeff.all
    P = proj.all
    mul!(buffer, U', Z)
    ldiv!(Σ, buffer)
    mul!(buffer, V', P, -one(T), one(T))
    lmul!(Ψ, buffer)
    mul!(B, V, buffer)
    axpy!(one(T), P, B)

    return nothing
end

# Apply one update in reguarlized problem.
function __reg_iterate__(::MMSVD, problem::MVDAProblem, ϵ, λ, extras)
    @unpack intercept, coeff = problem
    @unpack buffer = extras
    @unpack Z, Ψ, U, s, V = extras
    Σ = Diagonal(s)

    # need to compute Z via residuals...
    __evaluate_residuals__(problem, ϵ, extras, true, false, true)

    # Update parameters: B = V * Ψ * Σ⁻¹ * Uᵀ * Z
    B = coeff.all
    mul!(buffer, U', Z)
    ldiv!(Σ, buffer)
    lmul!(Ψ, buffer)
    mul!(B, V, buffer)

    return nothing
end


text:
function readDBpars( mpSim::MPsim, ws::WS, filePath::String )

    isConfigFromDB = !( ws["B12"] isa Missing ) && ( ws["B11"] == "YES" )
    configDBname = Base.source_path()
    configDBname = configDBname isa Nothing ? "" : dirname( configDBname )
    tmpDBname = ws["B4"]

    if tmpDBname isa Missing
        tmpDBname = ""
    else
        tmpDBname = string( joinpath( filePath,
            split( tmpDBname, r"[/\\]" )... ),
            endswith( tmpDBname, ".sqlite" ) ? "" : ".sqlite" )
        # ! The split( )... part ensures that the pathname is correct for both Windows and Linux systems.
    end  # if tmpDBname isa Missing

    if isConfigFromDB
        tmpConfigName = joinpath( split( sheet["B12"], r"[/\\]" )... )

        if tmpConfigName isa Missing
            @warn "No database entered to get configuration from. Configuring simulation from Excel sheet."
            isConfigFromDB = false
        else
            configDBname = string( joinpath( configDBname, tmpConfigName ),
                endswith( configDBname, ".sqlite" ) ? "" : ".sqlite" )

            if ispath( configDBname )
                isRerun = ( configDBname == tmpDBname ) &&
                    ( ws["B14"] == "YES" )
                configureSimFromDatabase( mpSim, configDBname, isRerun )

                # Don't make any further changes if the database is the same.
                if configDBname == tmpDBname
                    return "sameDB"
                end  # if configDBname == tmpDBname
            else
                @warn string( "Database '", configDBname,
                    "' does not exist. Configuring simulation from ",
                    "Excel sheet." )
                isConfigFromDB = false
            end  # if ispath( configDBname )
        end  # if tmpConfigName isa Missing
    end  # if isConfigFromDB

    # This block creates the database file if necessary and opens a link to it.
    mkpath( filePath )

    if mpSim.showInfo
        println( "Database file '", tmpDBname, "' ",
            isfile( tmpDBname ) ? "exists" : "does not exist", "." )
    end  # if mpSim.showInfo

    setSimulationDatabase!( mpSim, tmpDBname )

    # ! This line ensures that foreign key logic works.
    DBInterface.execute( mpSim.simDB, "PRAGMA foreign_keys = ON" )

    simName = ws["B5"] isa Missing ? "simulation" : string( ws["B5"] )
    setSimulationName!( mpSim, simName )

    # Check if databases are present and issue a warning if so.
    dbTableList = get( SQLite.tables( mpSim.simDB ), :name, Vector{String}() )

    if ( mpSim.persDBname ∈ dbTableList ) ||
        ( mpSim.histDBname ∈ dbTableList ) ||
        ( mpSim.transDBname ∈ dbTableList )
        @warn string( "Results for a simulation called '", mpSim.simName,
            "' already in database. These will be overwritten." )
    end  # if ( mpSim.personnelDBname ∈ dbTableList ) || ...

    return isConfigFromDB ? "otherDB" : "Excel"

end  # readDBpars( mpSim, ws, filePath )

function readGeneralPars( mpSim::MPsim, ws::WS, filePath::String )

    # Read the name of the catalogue file.
    catalogueName = ws["B3"] isa Missing ? "" : string( ws["B3"] )
    catalogueName = string( catalogueName,
        endswith( catalogueName, ".xlsx" ) ? "" : ".xlsx" )
    catalogueName = normpath( joinpath( filePath, "..", catalogueName ) )

    if !ispath( catalogueName )
        @error string( "Catalogue file '", catalogueName, "' does not exist." )
    end  # if !ispath( catalogueName )

    mpSim.catFileName = catalogueName

    # Set general simulation parameters.
    if ws["B6"] isa Integer
        setSimulationPersonnelTarget!( mpSim, ws["B6"] )
    end  # if ws["B6"] isa Integer

    if ws["B8"] isa Real
        setSimulationLength!( mpSim, ws["B8"] * 12.0 )
    end  # if ws["B8"] isa Real

    # setDatabaseCommitTime( mpSim, sheet["B8"] * 12.0 / sheet["B9"] )

end  # readGeneralPars( mpSim, ws, filePath )

function readGeneralPars( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    generalPars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'General'" ) ) )

    if isempty( generalPars )
        return
    end  # if isempty( generalPars )

    parInds = map( ["Sim name", "ID key", "Personnel target", "Sim length",
        "Current time", "DB commits"] ) do parName
        return findfirst( parName .== generalPars[:, :parName] )
    end  # map( ... ) do parName

    # Read simulation name.
    if parInds[1] isa Int
        setSimulationName!( mpSim, generalPars[parInds[1], :parValue] )
    end  # if parInds[1] isa Int

    # Read ID key.
    if parInds[2] isa Int
        setSimulationKey!( mpSim, generalPars[parInds[2], :parValue] )
    end  # if parInds[2] isa Int

    # Read global personnel target.
    if parInds[3] isa Int
        parVal = tryparse( Int, generalPars[parInds[3], :parValue] )

        if parVal isa Int
            setSimulationPersonnelTarget!( mpSim, parVal )
        end  # if parval isa Int
    end  # if parInds[3] isa Int

    # Read simulation length.
    if parInds[4] isa Int
        parVal = tryparse( Float64, generalPars[parInds[4], :parValue] )

        if parVal isa Float64
            setSimulationLength!( mpSim, parVal )
        end  # if parval isa Float64
    end  # if parInds[4] isa Int

    # Read current simulation time.
    if parInds[5] isa Int
        parVal = tryparse( Float64, generalPars[parInds[5], :parValue] )

        if ( parVal isa Float64 ) && ( parVal > 0.0 )
            run( mpSim.sim, parVal )
        end  # if parval isa Float64
    end  # if parInds[5] isa Int

    # Read number of database commits.
    if parInds[6] isa Int
        parVal = tryparse( Int, generalPars[parInds[6], :parValue] )

        if parVal isa Int
            # setSimulationPersonnelTarget!( mpSim, parVal )
        end  # if parval isa Int
    end  # if parInds[6] isa Int

end  # readGeneralPars( mpSim, configDB, configName )

function runSim( mpSim::MPsim, filename::String, showInfo::Bool )

    XLSX.openxlsx( filename ) do xf
        xs = xf["General"]

        # Run sim only if requested.
        if uppercase( xs["B14"] ) != "YES"
            return
        end  # if uppercase( xs["B14"] ) != "YES"

        seed = xs["B15"]
        seed = seed isa Integer ? max( -1, seed ) : -1
        sysEnt = uppercase( xs["B16"] ) == "YES"
        saveConfig = uppercase( xs["B17"] ) == "YES"
        run( mpSim, showInfo, saveConfig=saveConfig, seed=seed, sysEnt=sysEnt )
    end  # XLSX.openxlsx( filename ) do xf

    return true

end  # runSim( mpSim, filename )


text:
module MathToolkit

export recvec
export pslq

function recvec(v; compute=BigInt, output=Int64)
  n = length(v)
  l = Array(Rational{compute}, n)
  l[:] = v[:]
  z = div(n, 2)
  q1 = Array(Rational{compute}, n)
  q2 = Array(Rational{compute}, n)
  q1[1] = 0
  q2[1] = 1
  sq1 = 1
  sq2 = 1
  q  = Array(Rational{compute}, n)
  m  = Array(Rational{compute}, n)
  b = 1
  while sq2 <= z
    # look for first non zero term in the series taylor(l)
    while l[b] == 0
      if (b += 1) > n
        c = 1
        g = num(q2[1])
        for k = 1:sq2
          c = lcm(c, den(q2[k]))
          g = gcd(g, num(q2[k]))
        end
        if q2[1] < 0
          c = -c
        end
        q = Array(output, sq2)
        q[:] = (c//g)*q2[1:sq2]
        return q
      end
    end
    # reciprocal of the series: x^(b-1) / taylor(l)
    m[1] = 1 / l[b]
    for k = (b+1):n
      c = zero(compute)
      for j = b:(k-1)
        c -= l[j+1]*m[k-j]
      end
      m[ k-b+1 ] = c / l[b]
    end
    l, m = m, l
    n -= b-1
    # compute (only) denominator of Pade approximant
    q[1:sq2] = l[1] * q2[1:sq2]
    q[sq2+1:b+sq1-1] = 0
    q[b:b+sq1-1] += q1[1:sq1]
    sq = max(sq2, b+sq1-1)
    while q[sq] == 0
      sq -= 1
    end
    q1, q2, q = q2, q, q1
    sq1, sq2 = sq2, sq
    # substract term of degree 0 before computing next reciprocal
    b = 2
  end
  return Array(output, 0)
end

function pslq(x, prec; maxiter=256)
  const n = length(x)
  const teps = prec * 16.0
  const gam = 1.2
  s = Array(BigFloat, n)
  s[n] = abs(x[n])
  t = x[n]^2
  for i = (n-1):-1:1
    t += x[i]^2
    s[i] = sqrt(t)
  end
  t = s[1]
  y = x/t
  s = s/t

  h = [ zeros(BigFloat, n-1) for _ = 1:n ]
  for i = 1:n, j = 1:min(i,n-1)
    h[i][j] = (i==j)?s[j+1]/s[j]:-y[i]*y[j]/(s[j]*s[j+1])
  end

  a = [ zeros(Integer, n) for _ = 1:n ]
  for i = 1:n
    a[i][i] = 1
  end
  b = deepcopy(a)

  for i = 2:n, j = (i-1):-1:1
    t = round(h[i][j]/h[j][j])
    y[j] += t*y[i]
    for k =1:j
      h[i][k] -= t*h[j][k]
    end
    t = convert(Integer,t)
    for k = 1:n
      a[i][k] -= t*a[j][k]
      b[j][k] += t*b[i][k]
    end
  end

  m = 0
  for itr =1:maxiter
    mval = -1
    g = gam
    for i = 1:(n-1)
      t = abs(g*h[i][i])
      if t > mval
        mval = t
        m = i
      end
      g *= gam
    end
    y[m], y[m+1] = y[m+1], y[m]
    h[m], h[m+1] = h[m+1], h[m]
    a[m], a[m+1] = a[m+1], a[m]
    b[m], b[m+1] = b[m+1], b[m]
    if m < n-1
      t0 = sqrt(h[m][m]^2 + h[m][m+1]^2)
      t1 = h[m][m]/t0
      t2 = h[m][m+1]/t0
      for i = m:n
        t3 = h[i][m]
        t4 = h[i][m+1]
        h[i][m] = t1*t3 + t2*t4
        h[i][m+1] = t1*t4 - t2*t3
      end
    end
    for i = (m+1):n, j=min(i-1,m+1):-1:1
      t = round(h[i][j] / h[j][j])
      y[j] += t*y[i]
      for k=1:j
        h[i][k] -= t*h[j][k]
      end
      t = convert(Integer,t)
      for k=1:n
        a[i][k] -= t*a[j][k]
        b[j][k] += t*b[i][k]
      end
    end

    mval = 0
    for j=1:(n-1)
      t = abs(h[j][j])
      if t > mval
        mval = t
      end
    end
    for i=1:n
      t = abs(y[i])
      if t < teps
        return b[m]
      end
    end
  end
  return zeros(Integer, n)
end

end # module


text:
"""
    Interface to SNOPT. Modified from official repository
    (https://github.com/snopt/SNOPT7.jl) to have increased workspace size.
    SNOPT license is required.
"""

const libsnopt7 = "libsnopt7"

mutable struct snoptWorkspace
    status::Int

    # Workspace
    leniw::Int
    lenrw::Int
    iw::Vector{Int32}
    rw::Vector{Float64}

    leniu::Int
    lenru::Int
    iu::Vector{Int32}
    ru::Vector{Float64}

    x::Vector{Float64}
    lambda::Vector{Float64}
    obj_val::Float64

    num_inf::Int
    sum_inf::Float64

    iterations::Int
    major_itns::Int
    run_time::Float64

    function snoptWorkspace(leniw::Int,lenrw::Int)
        prob = new(0,leniw, lenrw,
                   zeros(Int,leniw), zeros(Float64,lenrw), 0, 0)
        finalizer(freeWorkspace!,prob)
        prob
    end
end

function freeWorkspace!(prob::snoptWorkspace)
    ccall((:f_snend, libsnopt7),
          Cvoid, (Ptr{Cint}, Cint, Ptr{Float64}, Cint),
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
end

# Exit codes
SNOPT_status = Dict(
    1=>:Solve_Succeeded,
    2=>:Feasible_Point_Found,
    3=>:Solved_To_Acceptable_Level,
    4=>:Solved_To_Acceptable_Level,
    5=>:Solved_To_Acceptable_Level,
    6=>:Solved_To_Acceptable_Level,
    11=>:Infeasible_Problem_Detected,
    12=>:Infeasible_Problem_Detected,
    13=>:Infeasible_Problem_Detected,
    14=>:Infeasible_Problem_Detected,
    15=>:Infeasible_Problem_Detected,
    16=>:Infeasible_Problem_Detected,
    21=>:Unbounded_Problem_Detected,
    22=>:Unbounded_Problem_Detected,
    31=>:Maximum_Iterations_Exceeded,
    32=>:Maximum_Iterations_Exceeded,
    34=>:Maximum_CpuTime_Exceeded,
    41=>:Numerical_Difficulties,
    42=>:Numerical_Difficulties,
    43=>:Numerical_Difficulties,
    44=>:Numerical_Difficulties,
    45=>:Numerical_Difficulties,
    51=>:User_Supplied_Function_Error,
    52=>:User_Supplied_Function_Error,
    56=>:User_Supplied_Function_Error,
    61=>:User_Supplied_Function_Undefined,
    62=>:User_Supplied_Function_Undefined,
    63=>:User_Supplied_Function_Undefined,
    71=>:User_Requested_Stop,
    72=>:User_Requested_Stop,
    73=>:User_Requested_Stop,
    74=>:User_Requested_Stop,
    81=>:Insufficient_Memory,
    82=>:Insufficient_Memory,
    83=>:Insufficient_Memory,
    91=>:Invalid_Problem_Definition,
    92=>:Invalid_Problem_Definition,
    999=>:Internal_Error)

# Callbacks
function obj_wrapper!(mode_::Ptr{Cint}, nnobj::Cint, x_::Ptr{Float64},
                      f_::Ptr{Float64}, g_::Ptr{Float64},
                      status::Cint)
    x    = unsafe_wrap(Array, x_, Int(nnobj))
    mode = unsafe_load(mode_)

    if mode == 0 || mode == 2
        obj = convert(Float64, eval_f(x)) :: Float64
        unsafe_store!(f_,obj)
    end

    if mode == 1 || mode == 2
        g = unsafe_wrap(Array,g_,Int(nnobj))
        eval_grad_f(x,g)
    end
    return
end

function con_wrapper!(mode_::Ptr{Cint}, nncon::Cint, nnjac::Cint, negcon::Cint,
                      x_::Ptr{Float64}, c_::Ptr{Float64}, J_::Ptr{Float64},
                      status::Cint)

    x    = unsafe_wrap(Array, x_, Int(nnjac))
    mode = unsafe_load(mode_)

    if mode == 0 || mode == 2
        c = unsafe_wrap(Array, c_, Int(nncon))
        eval_g(x, c)
    end

    if mode == 1 || mode == 2
        J = unsafe_wrap(Array,J_,Int(negcon))
        eval_jac_g(x, J)
    end
    return
end

# SNOPT7 routines
function initialize(printfile::String, summfile::String)
    prob = snoptWorkspace(30500,3000)

    ccall((:f_sninitx, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Ptr{UInt8}, Cint,
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          printfile, length(printfile), summfile, length(summfile),
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    return prob
end

function initialize(printfile::String, summfile::String,
                    leniw::Int, lenrw::Int)
    prob = snoptWorkspace(leniw, lenrw)

    ccall((:f_sninitx, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Cint,
           Ptr{UInt8}, Cint, Cint,
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          printfile, plen, prob.iprint, summfile, slen, prob.isumm,
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    return prob
end

# make larger workspace
function initialize(printfile::String, summfile::String,ws::Tuple)
    prob = snoptWorkspace(ws[1],ws[2])

    ccall((:f_sninitx, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Ptr{UInt8}, Cint,
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          printfile, length(printfile), summfile, length(summfile),
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    return prob
end

function readOptions(prob::snoptWorkspace, specsfile::String)
    status = [0]
    ccall((:f_snspecf, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Ptr{Cint},
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          specsfile, length(specsfile), status,
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    prob.status = status[1]
    return Int(prob.status)
end

function snopt!(prob::snoptWorkspace, start::String, name::String,
                m::Int, n::Int, nnCon::Int, nnObj::Int, nnJac::Int,
                fObj::Float64, iObj::Int,
                confun::Function, objfun::Function,
                #eval_f::Function, eval_grad_f::Function,
                #eval_g::Function, eval_jac_g::Function,
                J::SparseMatrixCSC, bl::Vector{Float64}, bu::Vector{Float64},
                hs::Vector{Int}, x::Vector{Float64})

    @assert n+m == length(x) == length(bl) == length(bu)
    @assert n+m == length(hs)

    prob.iu = [0]
    prob.ru = [0.]

    prob.x      = copy(x)
    prob.lambda = zeros(Float64,n+m)
    pi          = zeros(Float64,m)

    obj_callback = @cfunction($objfun, Cvoid,
                             (Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                              Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}))
    con_callback = @cfunction($confun, Cvoid,
                             (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                              Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}))

    valJ = J.nzval
    indJ = convert(Array{Cint}, J.rowval)
    locJ = convert(Array{Cint}, J.colptr)
    neJ  = length(valJ)

    status  = [0]
    nS      = [0]
    nInf    = [0]
    sInf    = [0.0]
    obj_val = [0.0]
    miniw   = [0]
    minrw   = [0]

    ccall((:f_snoptb, libsnopt7), Cvoid,
          (Ptr{UInt8}, Ptr{UInt8},
           Cint, Cint, Cint, Cint, Cint, Cint,
           Cint, Cdouble,
           Ptr{Cvoid}, Ptr{Cvoid},
           Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint},
           Ptr{Float64}, Ptr{Float64}, Ptr{Cint},
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cint}, Ptr{Cint},
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint,
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          start, name, m, n, neJ, nnCon, nnObj, nnJac,
          iObj, fObj,
          con_callback, obj_callback,
          valJ, indJ, locJ,
          bl, bu, hs, prob.x, pi, prob.lambda,
          status, nS, nInf, sInf, obj_val,
          miniw, minrw,
          prob.iu, prob.leniu, prob.ru, prob.lenru,
          prob.iw, prob.leniw, prob.rw, prob.lenrw)

    prob.status  = status[1]
    prob.obj_val = obj_val[1]

    prob.num_inf = nInf[1]
    prob.sum_inf = sInf[1]

    prob.iterations = prob.iw[421]
    prob.major_itns = prob.iw[422]

    prob.run_time   = prob.rw[462]

    return Int(prob.status)
end

function setOption!(prob::snoptWorkspace, optstring::String)
    # Set SNOPT7 option via string
    if !isascii(optstring)
        error("SNOPT7: Non-ASCII parameters not supported")
    end

    errors = [0]
    ccall((:f_snset, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Ptr{Cint},
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          optstring, length(optstring), errors,
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    return errors[1]
end

function setOption!(prob::snoptWorkspace, keyword::String, value::Int)
    # Set SNOPT7 integer option
    if !isascii(keyword)
        error("SNOPT7: Non-ASCII parameters not supported")
    end

    errors = [0]
    ccall((:f_snseti, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Cint, Ptr{Cint},
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          optstring, length(optstring), value, errors,
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    return errors[1]
end

function setOption!(prob::snoptWorkspace, keyword::String, value::Float64)
    # Set SNOPT7 real option
    if !isascii(keyword)
        error("SNOPT7: Non-ASCII parameters not supported")
    end

    errors = [0]
    ccall((:f_snseti, libsnopt7), Cvoid,
          (Ptr{UInt8}, Cint, Cdouble, Ptr{Cint},
           Ptr{Cint}, Cint, Ptr{Cdouble}, Cint),
          optstring, length(optstring), value, errors,
          prob.iw, prob.leniw, prob.rw, prob.lenrw)
    return errors[1]
end

mutable struct VariableData
    lower::Float64
    upper::Float64
    x::Float64
    g::Vector{Float64}
end
VariableData() = VariableData(-Inf, Inf, 0.0, [])

mutable struct Optimizer <: MOI.AbstractOptimizer
    sense::MOI.OptimizationSense
    objective::Union{MOI.SingleVariable,MOI.ScalarAffineFunction{Float64},MOI.ScalarQuadraticFunction{Float64},Nothing}
    variables::Vector{VariableData}
    linear_constraints::Vector{VariableData}
    nonlin::Union{MOI.NLPBlockData,Nothing}

    n::Int
    nncon::Int
    nlcon::Int

    workspace::Union{snoptWorkspace,Nothing}
    options::Dict{String, Any}
end
Optimizer(;options...) = Optimizer(MOI.MIN_SENSE, nothing, [], [], nothing, 0, 0, 0, nothing, options)

function MOI.empty!(model::Optimizer)
    model.sense = MOI.MIN_SENSE
    model.objective = nothing
    empty!(model.variables)
    empty!(model.linear_constraints)
    model.nonlin = nothing
    model.n = 0
    model.nncon = 0
    model.nlcon = 0
end

function MOI.is_empty(model::Optimizer)
    return model.objective == nothing &&
           isempty(model.variables) &&
           isempty(model.linear_constraints) &&
           model.nonlin == nothing &&
           model.n == 0 && model.nncon == 0 && model.nlcon == 0
end

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    return MOI.Utilities.default_copy_to(model, src, copy_names)
end

# MOI objects that SNOPT supports
MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,::Type{MOI.VariableIndex}) = true  #?
MOI.supports(::Optimizer, ::MOI.RawParameter) = true

MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Interval{Float64}}) = true

MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.Interval{Float64}}) = true

MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.variables)

function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    model.options[p.name] = value
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.sense = sense
    return
end

# Linear/quadratic objective
function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction,
                 func::Union{MOI.SingleVariable, MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction})
    model.objective = func
    return
end

# Variables
function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, vi::MOI.VariableIndex, value::Float64)
    model.variables[vi.value].x = value
    return
end

function MOI.add_variable(model::Optimizer)
    push!(model.variables, VariableData())
    return MOI.VariableIndex(length(model.variables))
end

function MOI.add_variables(model::Optimizer, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, lt::MOI.LessThan{Float64})
    # Upper bound on x
    vi = v.variable
    model.variables[vi.value].upper = lt.upper
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, gt::MOI.GreaterThan{Float64})
    # Lower bound on x
    vi = v.variable
    model.variables[vi.value].lower = gt.lower
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, eq::MOI.EqualTo{Float64})
    # Fixed variable
    vi = v.variable
    model.variables[vi.value].lower = eq.value
    model.variables[vi.value].upper = eq.value
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, bds::MOI.Interval{Float64})
    # Lower/upper bounds on x
    vi = v.variable
    model.variables[vi.value].lower = bds.lower
    model.variables[vi.value].upper = bds.upper
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}(vi.value)
end

# Linear constraints
macro define_add_constraint(set_type)
    quote
        function MOI.add_constraint(model::Optimizer, v::MOI.ScalarAffineFunction{Float64},
                                    set::$set_type)
            # Bounds on constraint
            x = VariableData()
            if typeof(set) == MOI.LessThan{Float64}
                x.upper = set.upper - v.constant
            elseif typeof(set) == MOI.GreaterThan{Float64}
                x.lower = set.lower - v.constant
            elseif typeof(set) == MOI.EqualTo{Float64}
                x.lower = set.value - v.constant
                x.upper = set.value - v.constant
            elseif typeof(set) == MOI.Interval{Float64}
                x.lower = set.lower - v.constant
                x.upper = set.upper - v.constant
            end

            x.g = [ t.coefficient for t in v.terms ]

            push!(model.linear_constraints, x)
            return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},$set_type}(length(model.linear_constraints))
        end
    end
end
@define_add_constraint(MOI.LessThan{Float64})
@define_add_constraint(MOI.GreaterThan{Float64})
@define_add_constraint(MOI.EqualTo{Float64})
@define_add_constraint(MOI.Interval{Float64})

# Nonlinear constraints
function MOI.set(model::Optimizer, ::MOI.NLPBlock, nonlin::MOI.NLPBlockData)
    model.nonlin = nonlin
    return
end

# Solve
function MOI.optimize!(model::Optimizer;
        ws::Tuple=(),outfile::String="snopt.out")

    n     = length(model.variables)
    nncon = model.nonlin !== nothing ? length(model.nonlin.constraint_bounds) : 0
    nlcon = length(model.linear_constraints)

    m     = nncon + length(model.linear_constraints)

    model.n = n
    model.nncon = nncon
    model.nlcon = nlcon

    # Bounds/initial point
    #    Variables
    bl = [v.lower for v in model.variables]
    bu = [v.upper for v in model.variables]
    x  = [v.x     for v in model.variables]
    #    Nonlinear constraints
    if nncon > 0
        for bound in model.nonlin.constraint_bounds
            push!(bl, bound.lower)
            push!(bu, bound.upper)
            push!(x,  0.0)
        end
    end

    #    Linear constraints
    for v in model.linear_constraints
        push!(bl, v.lower)
        push!(bu, v.upper)
        push!( x, v.x)
    end

    # Nonlinear data
    if nncon > 0
        features = MOI.features_available(model.nonlin.evaluator)
        has_objgrad = (:Grad in features)
        has_jacobian = (:Jac in features)

        MOI.initialize(model.nonlin.evaluator, features)

        # Jacobian
        sparseJ = MOI.jacobian_structure(model.nonlin.evaluator)  # returns tuple of coordinates
        nnjac = n

        j_rows = [sparseJ[i][1] for i in 1:length(sparseJ)]
        j_cols = [sparseJ[i][2] for i in 1:length(sparseJ)]
        j_vals = zeros(length(sparseJ))
    else
        nnjac  = 0
        j_rows = []
        j_cols = []
        j_vals = []
    end

    #   Add linear constraint Jacobian
    i = nncon
    for v in model.linear_constraints
        i = i+1
        for j in 1:n
            push!(j_rows, i)
            push!(j_cols, j)
            push!(j_vals, v.g[j])
        end
    end
    J = sparse(j_rows, j_cols, j_vals)

    # Minimize or maximize?

    if model.nonlin !== nothing
        # Nonlinear constraints/objective
        nnobj = n

        # Callbacks
        function eval_obj(mode_::Ptr{Cint}, n_::Ptr{Cint}, x_::Ptr{Float64},
                          f_::Ptr{Float64}, g_::Ptr{Float64}, status::Ptr{Cint})
            mode = unsafe_load(mode_)
            n    = unsafe_load(n_)

            x    = unsafe_wrap(Array, x_, Int(n))

            if mode == 0 || mode == 2
                if model.nonlin.has_objective
                    obj = MOI.eval_objective(model.nonlin.evaluator, x)
                elseif model.objective !== nothing
                    obj = eval_function(model.objective, x)
                else
                    obj = 0.0
                end
                unsafe_store!(f_,obj)
            end

            if mode == 1 || mode == 2
                g = unsafe_wrap(Array,g_, Int(n))
                if model.nonlin.has_objective
                    MOI.eval_objective_gradient(model.nonlin.evaluator, g, x)
                end
            end
            return
        end

        function eval_con(mode_::Ptr{Cint},
                          nncon_::Ptr{Cint}, nnjac_::Ptr{Cint}, negcon_::Ptr{Cint},
                          x_::Ptr{Float64}, c_::Ptr{Float64}, J_::Ptr{Float64},
                          status::Ptr{Cint})
            mode   = unsafe_load(mode_)
            nncon  = unsafe_load(nncon_)
            nnjac  = unsafe_load(nnjac_)
            negcon = unsafe_load(negcon_)

            x = unsafe_wrap(Array, x_, Int(nnjac))
            if mode == 0 || mode == 2
                c = unsafe_wrap(Array, c_, Int(nncon))
                MOI.eval_constraint(model.nonlin.evaluator, c, x)
            end

            if mode == 1 || mode == 2
                J = unsafe_wrap(Array,J_, Int(negcon))

                # MOI returns Jacobian in order specified by j_rows, j_cols
                # SNOPT needs it in sparse-by-col format
                JJ = zeros(negcon)
                MOI.eval_constraint_jacobian(model.nonlin.evaluator, JJ, x)
                vals = sparse(j_rows,j_cols,JJ).nzval
                for k in 1:negcon
                    J[k] = vals[k]
                end
            end
            return
        end

        # SNOPT
        fobj = 0.0
        iobj = 0
        hs   = zeros(Int,n+m)

        # make larger workspace
        if ws == ()
            # set workspace size (10x recommended by manual)
            ws = (1000*(n+m),2000*(n+m))
        end
        model.workspace = SNOPT7.initialize(outfile,"screen",ws)

        # model.workspace = SNOPT7.initialize("snopt.out","screen")

        for (name,value) in model.options
            # Replace underscore with space
            keyword = replace(String(name), "_" => " ")
            SNOPT7.setOption!(model.workspace, "$keyword $value")
        end
        SNOPT7.setOption!(model.workspace, "Summary frequency 1")
        SNOPT7.setOption!(model.workspace, "Print frequency 1")

        start = "Cold"
        name  = "prob"
        status = SNOPT7.snopt!(model.workspace, start, name,
                              m, n, nncon, nnobj, nnjac,
                              fobj, iobj, eval_con, eval_obj,
                              J, bl, bu, hs, x)
        SNOPT7.freeWorkspace!(model.workspace)

    else
        # Linear/quadratic objective

    end

end

# Return info
function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.workspace === nothing
        return MOI.NO_SOLUTION
    end

    status = SNOPT_status[model.workspace.status]
    if status in (:Solve_Succeeded,
                  :Feasible_Point_Found)
        return MOI.LOCALLY_SOLVED
    elseif status ==:Infeasible_Problem_Detected
        return MOI.LOCALLY_INFEASIBLE
    elseif status ==:Unbounded_Problem_Detected
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif status == :Solved_To_Acceptable_Level
        return MOI.ALMOST_LOCALLY_SOLVED
    elseif status == :User_Requested_Stop
        return MOI.INTERRUPTED
    elseif status == :Maximum_Iterations_Exceeded
        return MOI.ITERATION_LIMIT
    elseif status == :Maximum_CpuTime_Exceeded
        return MOI.TIME_LIMIT
    elseif status == (:Invalid_Problem_Definition,
                      :User_Supplied_Function_Error,
                      :User_Supplied_Function_Undefined)
        return MOI.INVALID_MODEL
    elseif status == :Insufficient_Memory
        return MOI.MEMORY_LIMIT
    elseif status == :Numerical_Difficulties
        return MOI.NUMERICAL_ERROR
    else
        error("Unrecognized SNOPT7 status $status")
    end
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.workspace !== nothing) ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if model.workspace === nothing
        return MOI.NO_SOLUTION
    end
    status = SNOPT_status[model.workspace.status]
    if status == :Solve_Succeeded
        return MOI.FEASIBLE_POINT
    elseif status == :Feasible_Point_Found
        return MOI.FEASIBLE_POINT
    elseif status == :Solved_To_Acceptable_Level
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == :Infeasible_Problem_Detected
        return MOI.INFEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, ::MOI.DualStatus)
    if model.workspace === nothing
        return MOI.NO_SOLUTION
    end
    status = SNOPT_status[model.workspace.status]
    if status == :Solve_Succeeded
        return MOI.FEASIBLE_POINT
    elseif status == :Feasible_Point_Found
        return MOI.FEASIBLE_POINT
    elseif status == :Solved_To_Acceptable_Level
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == :Infeasible_Problem_Detected
        return MOI.UNKNOWN_RESULT_STATUS # unbounded?
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    if model.workspace === nothing
        @error("ObjectiveValue not available.")
    end
    return model.workspace.obj_val
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.workspace === nothing
        @error("VariablePrimal not available.")
    end
    return model.workspace.x[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,S}) where {S}
    if model.workspace === nothing
        @error("ConstraintPrimal not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    return model.workspace.x[vi.value]
end

macro define_constraint_primal(set_type)
    quote
        function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                         ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, $set_type})
            if model.workspace === nothing
                @error("ConstraintPrimal not available.")
            end
            if !(1 <= ci.value <= model.nlcon)
                @error("Invalid constraint index ", ci.value)
            end
            return model.workspace.x[ci.value + model.n+model.nncon]
        end
    end
end
@define_constraint_primal(MOI.LessThan{Float64})
@define_constraint_primal(MOI.GreaterThan{Float64})
@define_constraint_primal(MOI.EqualTo{Float64})
@define_constraint_primal(MOI.Interval{Float64})

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,S}) where {S}
    if model.workspace === nothing
        @error("ConstraintDual not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    return model.workspace.lambda[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}) where {S}
    if model.workspace === nothing
        @error("ConstraintDual not available.")
    end
    @assert 1 <= ci.value <= model.nlcon
    return model.workspace.lambda[ci.value + model.n + model.nncon]
end

function MOI.get(model::Optimizer, ::MOI.NLPBlockDual)
    if model.workspace === nothing
        @error("NLPBlockDual not available.")
    end
    return model.workspace.lambda[model.n+1:model.n+model.nncon]
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
## Angle type definition and elementary functions
type Angle
  #type for angles. of form o + p * pi where p is rational
  o::Rational{Int64}
  p::Rational{Int64}
end

Base.show(io::IO, a::Angle) = print(io, a.o, "+", a.p, "π")

function Base.:+(a::Angle, b::Angle)
  return Angle(a.o + b.o, a.p + b.p)
end

function Base.:-(a::Angle)
  return Angle(-a.o, -a.p)
end

function Base.:-(a::Angle, b::Angle)
  return a + -b
end

function Base.:*(a::Angle, b::Angle)
  if a.p == 0//1 || b.p == 0//1 #only one nonzero pi coefficient => closed under multiplication
    return Angle(a.o * b.o, (a.o * b.p + b.o * a.p)) #pi^2 squared coefficient (a.p*b.p) is always zero
  end
  error("Cannot multiply. Must have at least one zero p term for type stability")
end

function Base.:*(a::Angle, b::Rational)
  return Angle(a.o * b, a.p * b)
end

function Base.:*(a::Rational, b::Angle)
  return b*a
end

function Base.:*(a::Array{Angle,1}, b::Angle)
  return b.*a
end

function Base.:*(a::Angle, b::Array{Angle,1})
  return b*a
end

function Base.:/(a::Angle, b::Angle)
  if b.p == 0//1 #invertable
    return Angle(a.o / b.o, a.p / b.o)
  end
  if a.o == 0//1 && b.o == 0//1 #pi terms will cancel
    return Angle(a.p / b.p, 0//1)
  end
  warn("Answer not expressible as angle data typel.")
  return (a.o + a.p * pi)/(b.o + b.p * pi)
end

function Base.:/(a::Array{Angle,1}, b::Angle)
  return a ./ b
end

function Base.:/(a::Int, b::Angle)
  return Angle(a//1, 0//1)/b
end

function Base.inv(a::Angle)
  return 1/a
end

function Base.zero(::Type{Angle})
  return Angle(zero(Rational), zero(Rational))
end

function Base.zero(a::Angle)
  return zero(Angle)
end

function Base.one(::Type{Angle})
  return Angle(one(Rational), zero(Rational))
end

import Base.== #the other syntax was throwing errors
==(a::Angle, b::Angle) = a.o == b.o && a.p == b.p
==(a::Angle, n:: Int) = a.o == n && a.p == 0

function Base.transpose(a::Angle)
  return a
end

function Base.norm(a::Angle)
  return norm(a.o + a.p * pi)
end

function Base.convert(::Type{Angle}, a::Int)
  return Angle(a//1, 0//1)
end

## Side type definition and elementary functions
type Side
  #type for sides. of form q + r *sqrt(m)
  q::Rational{Int64}
  r::Rational{Int64}
  m::Int

  function Side(q, r, m)
    if !isnull(m) && issquare(m) != 0 #m is square
      new(q + r * issquare(m), 0//1, m) #simplify to rational
    else
      new(q, r, m)
    end
  end
end

Base.show(io::IO, a::Side) = print(io, a.q, "+", a.r, "√", a.m)

function issquare(m)
  #returns if m is a perfect square. m will always be small and > 2
  #returns 0 if not square and the square root if square
  i = 2 #m is always greater than 2
  while i^2 <= m
    if i^2 == m
      return i
    end
    i += 1
  end
  return 0
end

function mfix(a::Side, b::Side)
  #checks if the value for m is null
  if isnull(a.m)
    if isnull(b.m)
      error("No non-null values for m.")
    else
      a.m = b.m
    end
  else
    if isnull(b.m)
      b.m = a.m
    end
  end
    return a, b
end

function Base.:+(a::Side, b::Side)
  #could add error check if m is the same
  a, b = mfix(a, b)
  return Side(a.q + b.q, a.r + b.r, a.m)
end

function Base.:-(a::Side)
  return Side(-a.q, -a.r, a.m)
end

function Base.:-(a::Side, b::Side)
  return a + -b
end

function Base.:*(a::Side, b::Rational)
  return Side(a.q * b, a.r * b, a.m)
end

function Base.:*(a::Rational, b::Side)
  return b*a
end

function Base.:*(a::Side, b::Side)
  a, b = mfix(a, b)
  return Side(a.q * b.q + a.r * b.r * a.m, a.r * b.q + a.q * b.r, a.m)
end

function Base.:*(a::Array{Side,2}, b::Side)
  return a .* b
end

function Base.:/(a::Side, b::Side)
  a, b = mfix(a, b)
  return Side((a.q*b.q-a.r*b.r*a.m)/(b.q^2-a.m*b.r^2),
   -(a.q*b.r-a.r*b.q)/(b.q^2-a.m*b.r^2), a.m) #never divides by zero because m cannot be square
end

function Base.:/(a::Int, b::Side)
  return Side(a//1, 0, b.m) / b
end

function Base.:/(a::Array{Side,1}, b::Side)
  return a ./ b
end

function Base.isless(a::Side, b::Side)
  a, b = mfix(a, b)
  return (a.q + a.r * sqrt(a.m)) < (b.q + b.r * sqrt(b.m))
end

function Base.inv(a::Side)
  return 1/a
end

import Base.== #the other syntax was throwing errors
==(a::Side, b::Side) = a.q == b.q &&  a.r == b.r #does not check for m equality
==(a::Side, b::Int) = a.q == b && a.r == 0//1

function Base.convert(::Type{Side}, a::Int, m)
  return Side(a//1, 0//1, m)
end

function Base.convert(::Float64, a::Side)
  return (a.p + a.r * a.m)::Float64
end

function Base.convert(::Type{Array{Side,2}}, a::Array{Int64,2}, m::Int64)
  output = zeros(Side, size(a), mvalue=m)
  for i in 1:length(a)
    output[i] = convert(Side, a[i], m)
  end
  return output
end

function Base.zero(::Type{Side})
  return Side(zero(Rational), zero(Rational), Nullable{Int}())
end

function Base.zero(::Type{Side}, m::Int)
  return Side(zero(Rational), zero(Rational), m)
end

function Base.zero(a::Side)
  return zero(Side)
end

function Base.zeros(::Type{Side}, size; mvalue::Int=Nullable{Int64})
  return repmat([zero(Side, mvalue)], size...)
end

function Base.zeros(a::Array{Side})
  return zeros(Side, size(a), mvalue=findm(a))
  error("All m-values in matrix null.")
end

function Base.one(::Type{Side})
  return Side(one(Rational), zero(Rational), Nullable{Int}())
end

function Base.one(::Type{Side}, m::Int)
  return Side(one(Rational), zero(Rational), m)
end

function Base.transpose(a::Side)
  return a
end

function Base.norm(a::Side)
  return a.q + a.r * sqrt(a.m)
end

function findm(a::Array{Side})
 #finds the first non-null value for m in the matrix
  for i in a
    if !isnull(i.m)
      return i.m
    end
  end
end

function sideeye(::Type{Side}, n, m)
  #returns nxn identity matrix of size with m
  output = zeros(Side, (n,n), mvalue=m)
  for i in 1:n
    output[i,i] = one(Side, m)
  end
  return output
end

## rref
function rref(mat)
  #returns modified reduced row eschelon form of mat
  #probably could be improved further
  j, k = size(mat)
  for i in 1:min(j,k) #loop through rows
    if mat[i,:] != zeros(mat[i,:]) #non zero row.
      pivotind = findfirst(mat[i,:])
      mat[i,:] /= mat[i,pivotind] #normalize leading element
      minus = repmat(mat[i,:]', j, 1) #base to zero the other rows
      minus .*= mat[:,pivotind] #scale to zero elements below pivot
      minus[i,:] = zeros(minus[i,:]) #zero pivot row to avoid cancelling it
      mat -= minus
    end
  end
  return rrefsort(mat)

end

function rrefsort(mat)
  #sorts rows to their aligned pivot location
  output = zeros(mat) #storage for output
  j, k = size(mat)
  for i in 1:min(j, k) #loop through each row
   if mat[i,:] != zeros(mat[i,:]) #non zero row.
    output[findfirst(mat[i,:]),:] = mat[i,:]
   end
  end
  return output
end

max_stars_repo_name: ericbumbalough/Reptiles

max_issues_repo_name: ericbumbalough/Reptiles

max_forks_repo_name: ericbumbalough/Reptiles


text:
module TestJaynes

using Test

include("../src/Jaynes.jl")
using .Jaynes
using Distributions

@time @testset "Execution contexts." begin
    include("contexts.jl")
end

@time @testset "Vectorized interfaces." begin
    include("vectorized.jl")
end

@time @testset "Importance sampling." begin
    include("importance_sampling.jl")
end

@time @testset "Particle filtering." begin
    include("particle_filter.jl")
end

@time @testset "Metropolis-Hastings." begin
    include("metropolis_hastings.jl")
end

@time @testset "Black-box extensions." begin
    include("blackbox.jl")
end

@time @testset "Selection query language." begin
    include("selection_query_language.jl")
end

@time @testset "Gradient learning." begin
    include("backpropagation.jl")
end

end #module


text:
#--------------------------------------
## LOAD PACKAGES
using RegNeuralDE, OrdinaryDiffEq, Flux, DiffEqFlux, Tracker
using YAML, Dates, BSON, Random, Statistics, Printf
using CUDA
using Flux: @functor
import Base.show

CUDA.allowscalar(false)
#--------------------------------------

#--------------------------------------
## CONFIGURATION
config_file = joinpath(pwd(), "experiments", "configs", "latent_ode.yml")
config = YAML.load_file(config_file)

Random.seed!(config["seed"])

hparams = config["hyperparameters"]
BATCH_SIZE = hparams["batch_size"]
EPOCHS = hparams["epochs"]
REGULARIZE = hparams["regularize"]
REG_TYPE = hparams["type"]
STEER = hparams["steer"]
identifier =
    REGULARIZE ? "$(string(now()))_$(REGULARIZE)_$(REG_TYPE)" : "$(string(now()))_vanilla"
EXPERIMENT_LOGDIR = joinpath(pwd(), "results", "latent_ode", identifier)
MODEL_WEIGHTS = joinpath(EXPERIMENT_LOGDIR, "weights.bson")
FILENAME = joinpath(EXPERIMENT_LOGDIR, "results.yml")

# Create a directory to store the results
isdir(EXPERIMENT_LOGDIR) || mkpath(EXPERIMENT_LOGDIR)
cp(config_file, joinpath(EXPERIMENT_LOGDIR, "config.yml"))
#--------------------------------------

#--------------------------------------
## NEURAL NETWORK
# Latent GRU Model
struct LatentGRU{U,R,N,D}
    update_gate::U
    reset_gate::R
    new_state::N
    latent_dim::D
end

function LatentGRU(in_dim::Int, h_dim::Int, latent_dim::Int)
    update_gate = Chain(
        Dense(latent_dim * 2 + in_dim * 2 + 1, h_dim, CUDA.tanh),
        Dense(h_dim, latent_dim, σ),
    )
    reset_gate = Chain(
        Dense(latent_dim * 2 + in_dim * 2 + 1, h_dim, CUDA.tanh),
        Dense(h_dim, latent_dim, σ),
    )
    new_state = Chain(
        Dense(latent_dim * 2 + in_dim * 2 + 1, h_dim, CUDA.tanh),
        Dense(h_dim, latent_dim * 2),
    )
    return LatentGRU(update_gate, reset_gate, new_state, latent_dim)
end

@functor LatentGRU

function single_run(p::LatentGRU, y_mean, y_std, x)
    # y -> latent_dim x B
    # x -> in_dim x B
    # x is the concatenation of a data and mask and Δtime
    y_concat = vcat(y_mean, y_std, x)

    update_gate = p.update_gate(y_concat)
    reset_gate = p.reset_gate(y_concat)

    concat = vcat(y_mean .* reset_gate, y_std .* reset_gate, x)

    new_state = p.new_state(concat)
    new_state_mean = new_state[1:p.latent_dim, :]
    # The easy-neural-ode paper uses abs, which is indeed a strange choice.
    # Instead we treat this value as logσ²
    new_state_std = new_state[p.latent_dim+1:end, :]

    new_y_mean = @. (1 - update_gate) * new_state_mean + update_gate * y_mean
    new_y_std = @. (1 - update_gate) * new_state_std + update_gate * y_std

    mask = sum(x[(size(x, 1)÷2+1):end, :], dims = 1) .> 0

    new_y_mean = @. mask * new_y_mean + (1 - mask) * y_mean
    new_y_std = @. mask * new_y_std + (1 - mask) * y_std

    return new_y_mean, new_y_std
end

function (p::LatentGRU)(x::AbstractArray{T,3}) where {T}
    z = TrackedArray(CUDA.zeros(T, p.latent_dim, size(x, 3)))
    y_mean, y_std = z, z
    for t = size(x, 2):-1:1
        y_mean, y_std = single_run(p, y_mean, y_std, x[:, t, :])
    end
    return vcat(y_mean, y_std)
end
#--------------------------------------

#--------------------------------------
## MODEL + DATASET + TRAINING UTILS
# Get the dataset
train_dataloader, test_dataloader =
    load_physionet(BATCH_SIZE, "data/physionet.bson", 0.8, x -> cpu(x))

opt = Flux.Optimise.Optimiser(InvDecay(1e-5), AdaMax(0.01))

# Setup the models
gru_rnn = LatentGRU(37, 40, 50) |> track |> gpu
rec_to_gen = Chain(Dense(100, 50, CUDA.tanh), Dense(50, 2 * 20)) |> track |> gpu
gen_dynamics =
    Chain(
        x -> CUDA.tanh.(x),
        Dense(20, 50, CUDA.tanh),
        Dense(50, 20, CUDA.tanh),
        Dense(20, 50, CUDA.tanh),
        Dense(50, 20, CUDA.tanh),
        Dense(20, 50, CUDA.tanh),
        Dense(50, 20, CUDA.tanh),
        Dense(20, 50, CUDA.tanh),
        Dense(50, 20, CUDA.tanh),
    ) |>
    track |>
    gpu

# AutoTsit5(Tsit5()) is simply Tsit5() since we don't want to switch to a
# stiff solver. This "hack" allows us to construct a CompositeAlgorithm and
# allows us to get the stiffness estimate from the solver itself.
solver =
    REGULARIZE ?
    (
        REG_TYPE == "stiff_est" || REG_TYPE == "error_stiff_est" ? AutoTsit5(Tsit5()) :
        Tsit5()
    ) : Tsit5()
saveat = train_dataloader.data[5][1, :, 1] |> f32
node = TrackedNeuralODE(
    gen_dynamics,
    [0.0f0, 1.0f0],
    false,
    REGULARIZE,
    solver,
    saveat = saveat,
    reltol = 1.4f-8,
    abstol = 1.4f-8,
)
gen_to_data = Dense(20, 37) |> track |> gpu

model = LatentTimeSeriesModel(gru_rnn, rec_to_gen, node, gen_to_data)
ps = Flux.trainable(model)

agg = mean
if REG_TYPE == "error_est"
    # Anneal the regularization so that it doesn't overpower the
    # the main objective
    λᵣ₀ = 1.0f3
    λᵣ₁ = 1.0f2
    save_func(u, t, integrator) = integrator.EEst * integrator.dt
    global agg = mean
elseif REG_TYPE == "stiff_est"
    # No annealing is generally needed for stiff_est
    λᵣ₀ = 1.0f1
    λᵣ₁ = 1.0f1
    const stability_size =
        Tracker.TrackedReal(1 / Float32(OrdinaryDiffEq.alg_stability_size(Tsit5())))
    function save_func(u, t, integrator)
        stiff_est = abs(integrator.eigen_est)
        return stability_size * ((iszero(stiff_est) || isnan(stiff_est)) ? 0 : stiff_est)
    end
    global agg = maximum
elseif REG_TYPE == "error_stiff_est"
    λᵣ₀ = 1.0f1
    λᵣ₁ = 1.0f1
    const mul_val = Tracker.TrackedReal(1.0f0)
    const stability_size =
        Tracker.TrackedReal(1 / Float32(OrdinaryDiffEq.alg_stability_size(Tsit5())))
    function save_func(u, t, integrator)
        err_est = integrator.EEst * integrator.dt
        eest = Tracker.data(err_est)
        stiff_est = integrator.eigen_est
        sest = Tracker.data(stiff_est)
        return (
            ((iszero(eest) || isnan(eest)) ? 0 : err_est) +
            0.1f0 * stability_size * ((iszero(sest) || isnan(sest)) ? 0 : stiff_est)
        ) * mul_val
    end
    global agg = mean
end
kᵣ = log(λᵣ₀ / λᵣ₁) / EPOCHS
# Exponential Decay
λᵣ_func(t) = λᵣ₀ * exp(-kᵣ * t)

λₖ_func(t) = max(0, 1 - 0.99f0^(t - 10))

## Steer Regularization
sample_tbounds(t, δt) = clamp!(
    vcat(t[1:1, :], t[2:end, :] .+ (2 .* rand(eltype(δt), size(δt)) .- 1) .* δt ./ 2),
    0.0f0,
    1.0f0,
)

function sample_tbounds(t)
    δt = t[2:end, :] .- t[1:end-1, :] .+ eps(Float32)
    return sample_tbounds(t, δt), δt
end

_, δsaveat = sample_tbounds(saveat)

## Loss Functions
function log_likelihood(∇pred, mask)
    function _sum_stable_infer(x)::typeof(x)
        return sum(x, dims = (1, 2))
    end
    σ_ = 0.01f0
    sample_likelihood =
        -CUDA.pow.(∇pred, 2) ./ (2 * σ_^2) .- log(σ_) .- log(Float32(2π)) ./ 2
    return reshape(_sum_stable_infer(sample_likelihood) ./ _sum_stable_infer(mask), :)
end

# Holds only for a Standard Gaussian Prior
kl_divergence(μ, logσ²) =
    reshape(mean(exp.(logσ²) .+ CUDA.pow.(μ, 2) .- 1 .- logσ², dims = 1) ./ 2, :)

function loss_function(
    data,
    mask,
    _t,
    model,
    p1,
    p2,
    p3,
    p4;
    λᵣ = 1.0f2,
    λₖ = 1.0f0,
    notrack = false,
    saveat = nothing,
)
    x_ = vcat(data, mask, _t)
    result, μ₀, logσ², nfe, sv =
        model(x_, p1, p2, p3, p4; func = save_func, saveat = saveat)

    data_ = data .* mask
    pred_ = result .* mask
    ∇pred = pred_ .- data_

    _log_likelihood = log_likelihood(∇pred, mask)
    _kl_div = λₖ .* kl_divergence(μ₀, logσ²)
    reg = REGULARIZE ? λᵣ * agg(sv.saveval) : zero(eltype(pred_))
    total_loss = -mean(_log_likelihood .- _kl_div) + reg

    if !notrack
        ll_un = -mean(_log_likelihood |> untrack)
        kl_un = mean(_kl_div |> untrack)
        rg_un = reg |> untrack
        tl_un = total_loss |> untrack
        logger(
            false,
            Dict(
                "Total Loss" => tl_un,
                "Negative Log Likelihood" => ll_un,
                "KL Divergence" => kl_un,
                "Regularization" => rg_un,
            ),
        )
    end

    return total_loss
end

function total_loss_on_dataset(model, dataloader)
    loss = 0.0f0
    count = 0
    for (i, (d_, m_, _, _, _, _)) in enumerate(dataloader)
        d = d_ |> track |> gpu
        m = m_ |> track |> gpu

        x_ = vcat(d, m, _t)
        result, _, _, _, _ = model(x_)

        data_ = d .* m
        pred_ = result .* m
        ∇pred = pred_ .- data_

        count += size(d, 3)
        loss += sum(sum(∇pred .^ 2, dims = (1, 2)) ./ sum(m, dims = (1, 2))) |> untrack

        d = m = nothing
        GC.gc(true)
    end
    return loss ./ count
end
#--------------------------------------

#--------------------------------------
## LOGGING UTILITIES
nfe_counts = Vector{Float64}(undef, EPOCHS + 1)
train_loss = Vector{Float64}(undef, EPOCHS + 1)
test_loss = Vector{Float64}(undef, EPOCHS + 1)
train_runtimes = Vector{Float64}(undef, EPOCHS + 1)
inference_runtimes = Vector{Float64}(undef, EPOCHS + 1)
train_runtimes[1] = 0

logger = table_logger(
    [
        "Epoch Number",
        "NFE Count",
        "Train Loss",
        "Test Loss",
        "Train Runtime",
        "Inference Runtime",
    ],
    ["Total Loss", "Negative Log Likelihood", "KL Divergence", "Regularization"],
)
#--------------------------------------

#--------------------------------------
## TESTING THE MODEL
# Dummy Input for first run
d, m, _, _, t_d, _ = iterate(train_dataloader)[1]
d = d |> gpu
m = m |> gpu
function get_t_saveat(t = nothing)
    if !STEER
        t = t |> f32
        tt = saveat
    else
        tt = sample_tbounds(saveat, δsaveat)
        t = repeat(reshape(tt, 1, :, 1), 1, 1, BATCH_SIZE)
    end
    _t = hcat(t[:, 2:end, :] .- t[:, 1:end-1, :], zeros(1, 1, size(t, 3))) |> gpu
    return t, tt, _t
end
t, tt, _t = get_t_saveat(t_d)
x_ = vcat(d, m, _t |> track)
dummy_data = x_
stime = time()
result, μ₀, logσ², _nfe, sv = model(x_; saveat = tt)
inference_runtimes[1] = time() - stime
train_runtimes[1] = 0.0
nfe_counts[1] = _nfe
train_loss[1] = total_loss_on_dataset(model, train_dataloader)
test_loss[1] = total_loss_on_dataset(model, test_dataloader)

logger(
    false,
    Dict(),
    0.0,
    nfe_counts[1],
    train_loss[1],
    test_loss[1],
    train_runtimes[1],
    inference_runtimes[1],
)

t, tt, _t = get_t_saveat(t_d)
loss_function(d, m, _t |> track, model, ps...; notrack = true, saveat = tt)

Tracker.gradient(
    (p1, p2, p3, p4) -> loss_function(
        d,
        m,
        _t |> track,
        model,
        p1,
        p2,
        p3,
        p4;
        notrack = true,
        saveat = tt,
    ),
    ps...,
)
#--------------------------------------

#--------------------------------------
## TRAINING
for epoch = 1:EPOCHS
    λᵣ = λᵣ_func(epoch - 1)
    λₖ = λₖ_func(epoch - 1)

    timing = 0

    for (i, (d_, m_, _, _, t_d, _)) in enumerate(train_dataloader)
        local d = d_ |> gpu |> track
        local m = m_ |> gpu |> track
        # If STEER is true then we get stochastic saveat positions
        local t, tt, _t = get_t_saveat(t_d)

        start_time = time()
        gs = Tracker.gradient(
            (p1, p2, p3, p4) -> loss_function(
                d,
                m,
                _t |> track,
                model,
                p1,
                p2,
                p3,
                p4;
                λᵣ = λᵣ,
                λₖ = λₖ,
                notrack = false,
                saveat = tt,
            ),
            ps...,
        )
        update_parameters!(ps, gs, opt)
        timing += time() - start_time

        d = m = nothing
        GC.gc(true)
    end

    # Record the time per epoch
    train_runtimes[epoch+1] = timing

    # Record the NFE count
    start_time = time()
    _, _, _, nfe, _ = model(dummy_data)
    inference_runtimes[epoch+1] = time() - start_time
    nfe_counts[epoch+1] = nfe

    # Test and Train Accuracy
    train_loss[epoch+1] = total_loss_on_dataset(model, train_dataloader)
    test_loss[epoch+1] = total_loss_on_dataset(model, test_dataloader)

    logger(
        false,
        Dict(),
        epoch,
        nfe_counts[epoch+1],
        train_loss[epoch+1],
        test_loss[epoch+1],
        train_runtimes[epoch+1],
        inference_runtimes[epoch+1],
    )
end
logger(true, Dict())
#--------------------------------------

#--------------------------------------
## STORE THE RESULTS
results = Dict(
    :nfe_counts => nfe_counts,
    :train_loss => train_loss,
    :test_loss => test_loss,
    :train_runtimes => train_runtimes,
    :inference_runtimes => inference_runtimes,
)

weights = Flux.params(model) .|> cpu .|> untrack
BSON.@save MODEL_WEIGHTS weights

YAML.write_file(FILENAME, results)
#--------------------------------------


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# simplelazy2.jl
#
# Builds on the simplelazy example by demonstrating how the logic of
# determing a new constraint can be seperated from the logic of the callback
# handler itself.
#############################################################################

using JuMP
using Gurobi
using Base.Test

function cornerChecker(x_val, y_val)
    # This function does not depend on the model, and could
    # be written anywhere. Instead, it returns a tuple of
    # values (newcut, x_coeff, y_coeff, rhs) where newcut is a
    # boolean if a cut is needed, x_coeff is the coefficient
    # on the x variable, y_coeff is the coefficient on
    # the y variable, and rhs is the right hand side
    TOL = 1e-6
    if y_val - x_val > 1 + TOL
        return (true, -1.0, 1.0, 1.0)  # Top left
    elseif y_val + x_val > 3 + TOL
        return (true,  1.0, 1.0, 3.0)  # Top right
    else
        return (false, 0.0, 0.0, 0.0)  # No cut
    end
end

# A unit test for the cornerChecker function
function test_cornerChecker()
    # Test the four corners - only two should produce cuts
    newcut, x_coeff, y_coeff, rhs = cornerChecker(0, 0)
    @test !newcut

    newcut, x_coeff, y_coeff, rhs = cornerChecker(2, 0)
    @test !newcut

    newcut, x_coeff, y_coeff, rhs = cornerChecker(0, 2)
    @test newcut
    @test x_coeff == -1.0
    @test y_coeff ==  1.0
    @test rhs == 1.0

    newcut, x_coeff, y_coeff, rhs = cornerChecker(2, 2)
    @test newcut
    @test x_coeff ==  1.0
    @test y_coeff ==  1.0
    @test rhs == 3.0
end

function solveProblem()
    m = Model(solver=GurobiSolver())

    @variable(m, 0 <= x <= 2, Int)
    @variable(m, 0 <= y <= 2, Int)
    @objective(m, Max, y)

    # Note that the callback is now a stub that passes off
    # the work to the "algorithm"
    function corners(cb)
        x_val = getvalue(x)
        y_val = getvalue(y)
        println("In callback function, x=$x_val, y=$y_val")

        newcut, x_coeff, y_coeff, rhs = cornerChecker(x_val, y_val)

        if newcut
            @lazyconstraint(cb, x_coeff*x + y_coeff*y <= rhs)
        end
    end  # End of callback function

    addlazycallback(m, corners)
    solve(m)
    println("Final solution: [ $(getvalue(x)), $(getvalue(y)) ]")
end

# Run tests
test_cornerChecker()

# Solve it
solveProblem()


text:
"""
    pcdc!(rbm, data)

Trains the RBM on data using Persistent Contrastive divergence, with centered gradients.
See:

J. Melchior, A. Fischer, and L. Wiskott. JMLR 17.1 (2016): 3387-3447.
<http://jmlr.org/papers/v17/14-237.html>
"""
function pcdc!(rbm::RBM, data::AbstractArray;
    batchsize = 1,
    epochs = 1,
    optim = Flux.ADAM(),
    history::MVHistory = MVHistory(),
    wts = nothing,
    steps::Int = 1,
    center_h::Bool = true, # center hidden?
    damping::Real = true # running average for λh
)
    @assert size(data) == (size(rbm.visible)..., size(data)[end])
    @assert isnothing(wts) || _nobs(data) == _nobs(wts)

    # data statistics
    stats = RBMs.suffstats(rbm.visible, data; wts)
    avg_h = RBMs.batchmean(rbm.hidden, RBMs.mean_h_from_v(rbm, data); wts)

    # initialize fantasy chains by sampling visible layer
    vm = RBMs.transfer_sample(rbm.visible, falses(size(rbm.visible)..., batchsize))

    for epoch in 1:epochs
        batches = minibatches(data, wts; batchsize = batchsize)
        Δt = @elapsed for (vd, wd) in batches
            # update fantasy chains
            vm = RBMs.sample_v_from_v(rbm, vm; steps = steps)
            # compute centered gradients
            ∂ = ∂contrastive_divergence_centered(
                rbm, vd, vm; wd, stats, center_h, damping, avg_h
            )
            # update parameters using gradient
            RBMs.update!(optim, rbm, ∂)
            # store gradient norms
            push!(history, :∂, RBMs.gradnorms(∂))
        end
        push!(history, :epoch, epoch)
        push!(history, :Δt, Δt)
        @debug "epoch $epoch/$epochs ($(round(Δt, digits=2))s)"
    end

    return history
end

function ∂contrastive_divergence_centered(
    rbm::RBM, vd::AbstractArray, vm::AbstractArray;
    wd = nothing, wm = nothing,
    stats = RBMs.suffstats(rbm.visible, vd; wts = wd),
    center_h::Bool = true, damping::Real = 0.01, avg_h::AbstractArray
)
    ∂d = RBMs.∂free_energy(rbm, vd; wts = wd, stats)
    ∂m = RBMs.∂free_energy(rbm, vm; wts = wm)
    ∂ = RBMs.subtract_gradients(∂d, ∂m)

    # extract moment estimates from the gradients
    λv = grad2mean(visible(rbm), ∂d.visible)   # <v>_d, uses full data from sufficient_statistics
    λh = grad2mean(hidden(hidden), ∂d.hidden)  # <h>_d, uses minibatch

    # since <h>_d uses minibatches, we keep a running average
    @assert size(avg_h) == size(λh)
    λh .= avg_h .= (1 - damping) * avg_h .+ damping .* λh

    ∂c = center_gradients(rbm, ∂, λv, center_h * λh)
    return ∂c
end

function center_gradients(rbm::RBM, ∂::NamedTuple, λv::AbstractArray, λh::AbstractArray)
    @assert size(rbm.visible) == size(λv)
    @assert size(rbm.hidden) == size(λh)
    @assert size(∂.w) == size(rbm.w)

    ∂wmat = reshape(∂.w, length(rbm.visible), length(rbm.hidden))
    ∂cwmat = ∂wmat - vec(λv) * vec(∂.hidden.θ)' - vec(∂.visible.θ) * vec(λh)'
    ∂cw = reshape(∂cwmat, size(rbm.w))

    shift_v = reshape(∂cwmat  * vec(λh), size(rbm.visible))
    shift_h = reshape(∂cwmat' * vec(λv), size(rbm.hidden))
    ∂cv = center_gradients(rbm.visible, ∂.visible, shift_v)
    ∂ch = center_gradients(rbm.hidden,  ∂.hidden,  shift_h)

    return (visible = ∂cv, hidden = ∂ch, w = ∂cw)
end

function center_gradients(
    layer::Union{Binary,Spin,Potts,Gaussian,ReLU,pReLU,xReLU},
    ∂::NamedTuple, λ::AbstractArray
)
    @assert size(layer) == size(∂.θ) == size(λ)
    return (∂..., θ = ∂.θ - λ)
end

function center_gradients(layer::dReLU, ∂::NamedTuple, λ::AbstractArray)
    @assert size(layer) == size(∂.θp) == size(∂.θn) == size(λ)
    return (θp = ∂.θp - λ, θn = ∂.θn - λ, γp = ∂.γp, γn = ∂.γn)
end

# get moments from layer gradients, e.g. <v> = -derivative w.r.t. θ
grad2mean(::Union{Binary,Spin,Potts,Gaussian,ReLU,pReLU,xReLU}, ∂::NamedTuple) = -∂.θ
grad2mean(::dReLU, ∂::NamedTuple) = -(∂.θp + ∂.θn)