text:
using HDF5, JLD, Base.Test, Compat
import Compat.String

# Define variables of different types
x = 3.7
A = reshape(collect(1:15), 3, 5)
str = "Hello"
stringsA = Compat.ASCIIString["It", "was", "a", "dark", "and", "stormy", "night"]
stringsU = Compat.UTF8String["It", "was", "a", "dark", "and", "stormy", "night"]
empty_string = ""
empty_string_array = Compat.ASCIIString[]
empty_array_of_strings = Compat.ASCIIString[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
@compat c = Float32(3)+Float32(7)im
cint = 1+im  # issue 108
C = reinterpret(Complex128, B, (4,))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
try
    global MyStruct
    type MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
d = Dict([(syms[1], "aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
T = UInt8
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]] # issue #123
typevar = Array{Int}[[1]]
typevar_lb = Vector{TypeVar(:U, Integer)}[[1]]
typevar_ub = Vector{TypeVar(:U, Int, Any)}[[1]]
typevar_lb_ub = Vector{TypeVar(:U, Int, Real)}[[1]]
undef = cell(1)
undefs = cell(2, 2)
ms_undef = MyStruct(0)
# Immutable type:
rng = 1:5
# Type with a pointer field (#84)
objwithpointer = r"julia"
# Custom BitsType (#99)
bitstype 64 MyBT
bt = reinterpret(MyBT, 55)
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = sub([1:5;], 1:5)
# Array of empty tuples (to test tuple type params)
@compat arr_empty_tuple = Tuple{}[]

iseq(x,y) = isequal(x,y)
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(c1::Array{Base.Sys.CPUinfo}, c2::Array{Base.Sys.CPUinfo}) = length(c1) == length(c2) && all([iseq(c1[i], c2[i]) for i = 1:length(c1)])
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in Base.Sys.CPUinfo.names
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                warn("Error reading ", $(string(sym)))
                rethrow(e)
            end
            if !iseq(tmp, $sym)
                error("For ", $(string(sym)), ", read value does not agree with written value")
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
        end
    end
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @assert a == b
function checkexpr(a::Expr, b::Expr)
    @assert a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if (isa(a.args[i], Expr) && a.args[i].head == :line) || isa(a.args[i], LineNumberNode)
            i += 1
            continue
        end
        if (isa(b.args[j], Expr) && b.args[j].head == :line) || isa(b.args[j], LineNumberNode)
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(a.args) && j >= length(b.args)
end

for fn in ("v0.2.26.jld", "v0.2.28.jld")
    for mmap = (true, false)
        fidr = jldopen(joinpath(splitdir(@__FILE__)[1], fn), "r", mmaparrays=mmap)
        @check fidr x
        @check fidr A
        @check fidr str
        @check fidr stringsA
        @check fidr stringsU
        @check fidr empty_string
        @check fidr empty_string_array
        @check fidr empty_array_of_strings
        @check fidr tf
        @check fidr TF
        @check fidr AB
        @check fidr t
        @check fidr c
        @check fidr cint
        @check fidr C
        @check fidr emptyA
        @check fidr emptyB
        @check fidr ms
        @check fidr msempty
        @check fidr sym
        @check fidr syms
        @check fidr d
        exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
        checkexpr(ex, exr)
        @check fidr T
        @check fidr char
        @check fidr unicode_char
        @check fidr α
        @check fidr β
        @check fidr vv
        @check fidr rng
        @check fidr typevar
        @check fidr typevar_lb
        @check fidr typevar_ub
        @check fidr typevar_lb_ub

        # Special cases for reading undefs
        undef = read(fidr, "undef")
        if !isa(undef, Array{Any, 1}) || length(undef) != 1 || isdefined(undef, 1)
            error("For undef, read value does not agree with written value")
        end
        undefs = read(fidr, "undefs")
        if !isa(undefs, Array{Any, 2}) || length(undefs) != 4 || any(map(i->isdefined(undefs, i), 1:4))
            error("For undefs, read value does not agree with written value")
        end
        ms_undef = read(fidr, "ms_undef")
        if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
            error("For ms_undef, read value does not agree with written value")
        end

        @assert !in("objwithpointer", names(fidr))
        @check fidr bt
        @check fidr sa_asc
        @check fidr sa_utf8
        @check fidr arr_empty_tuple

        x1 = read(fidr, "group1/x")
        @assert x1 == Any[1]
        x2 = read(fidr, "group2/x")
        @assert x2 == Any[2]

        # load but don't check
        read(fidr, "cpus")

        # subarray (changed representation in julia 0.4)
        subarray_safe = JLD.JLD00.readsafely(fidr, "subarray")
        @test subarray_safe["indexes"] == (1:5,)
        @test subarray_safe["parent"] == [1:5;]

        close(fidr)
    end
end


text:
@inline function parseint(io::IO, len::Int)
    if len == 1
        return Int(Parsers.readbyte(io) - Parsers.ZERO)
    end
    v = zero(Int)
    b = Parsers.readbyte(io)
    negative = false
    if b === Parsers.MINUS # check for leading '-' or '+'
        negative = true
    elseif b !== Parsers.PLUS
        v = Int(b - Parsers.ZERO)
    end
    for _ = 1:(len - 1)
        v *= Int(10)
        v += Int(Parsers.readbyte(io) - Parsers.ZERO)
    end
    return ifelse(negative, -v, v)
end

function checkint2(str::String)
    len = sizeof(str)
    BUF.data = str
    BUF.ptr = 1
    BUF.size = len
    return checkint(BUF)
end

@inline function checkint(io::IO)
    eof(io) && return false
    b = Parsers.peekbyte(io)
    negative = false
    if b === Parsers.MINUS # check for leading '-' or '+'
        negative = true
        Parsers.readbyte(io)
        eof(io) && return false
        b = Parsers.peekbyte(io)
    elseif b === Parsers.PLUS
        Parsers.readbyte(io)
        eof(io) && return false
        b = Parsers.peekbyte(io)
    end
    parseddigits = false
    while Parsers.NEG_ONE < b < Parsers.TEN
        parseddigits = true
        b = Parsers.readbyte(io)
        eof(io) && break
        b = Parsers.peekbyte(io)
    end
    return parseddigits
end

ext: jl


text:
import Pkg

Pkg.activate(".")

import Wordle
import Test: @test, @test_throws, @testset

test_files = (
    "historical_words.jl",
    "filtering/filter_candidates.jl",
    "filtering/utils.jl",
    "filtering/update_constraints.jl",
    "filtering/satisfies_constraints.jl",
    "evaluation.jl",
    "core_types/candidate.jl",
    "core_types/constraints.jl",
    "core_types/result.jl",
    "core_types/trial.jl",
    "scoring/scoring.jl",
    "sourcing/candidate_pool.jl",
)

@testset begin
    for test_file in test_files
        include(test_file)
    end
end


text:
using ReachabilityAnalysis, SparseArrays, JLD2

LazySets.set_ztol(Float64, 1e-14)

const x25 = [zeros(24); 1.0; zeros(23)]
const x25e = vcat(x25, 0.0)

examples_dir = normpath(@__DIR__, "..", "..", "..", "examples")
building_path = joinpath(examples_dir, "Building", "building.jld2")

function building_BLDF01()
    @load building_path A B
    n = size(A, 1)
    U = Interval(0.8, 1.0)
    S = @system(x' = Ax + Bu, u ∈ U, x ∈ Universe(n))

    # initial states
    center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
    radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
    X0 = Hyperrectangle(center_X0, radius_X0)

    prob_BLDF01 = InitialValueProblem(S, X0)
end

using ReachabilityAnalysis: add_dimension

function building_BLDC01()
    @load building_path A B
    n = size(A, 1)
    U = Interval(0.8, 1.0)

    # initial states
    center_X0 = [fill(0.000225, 10); fill(0.0, 38)]
    radius_X0 = [fill(0.000025, 10); fill(0.0, 14); 0.0001; fill(0.0, 23)]
    X0 = Hyperrectangle(center_X0, radius_X0)

    Ae = add_dimension(A)
    Ae[1:n, end] = B
    prob_BLDC01 = @ivp(x' = Ae * x, x(0) ∈ X0 × U)
end

using Plots

prob_BLDF01 = building_BLDF01()

sol_BLDF01_dense = solve(prob_BLDF01, T=20.0, alg=LGG09(δ=0.004, vars=(25), n=48));

plot(sol_BLDF01_dense, vars=(0, 25), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0, xlab="t", ylab="x25")

ρ(x25, sol_BLDF01_dense)

ρ(x25, sol_BLDF01_dense) <= 5.1e-3 # BLDF01 - BDS01

ρ(x25, sol_BLDF01_dense) <= 4e-3 # BLDF01 - BDU01

ρ(x25, sol_BLDF01_dense(20.0))

ρ(x25, sol_BLDF01_dense(20.0)) <= -0.78e-3 # BLDF01 - BDU02

sol_BLDF01_discrete = solve(prob_BLDF01, T=20.0, alg=LGG09(δ=0.01, vars=(25), n=48, approx_model=NoBloating()));

plot(sol_BLDF01_discrete, vars=(0, 25), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0, xlab="t", ylab="x25")

ρ(x25, sol_BLDF01_discrete)

ρ(x25, sol_BLDF01_discrete) <= 5.1e-3 # BLDF01 - BDS01

ρ(x25, sol_BLDF01_discrete)

ρ(x25, sol_BLDF01_discrete) <= 4e-3 # BLDF01 - BDU01

ρ(x25, sol_BLDF01_discrete(20.0))

ρ(x25, sol_BLDF01_discrete(20.0)) <= -0.78e-3 # BLDF01 - BDU02

prob_BLDC01 = building_BLDC01()

sol_BLDC01_dense = solve(prob_BLDC01, T=20.0, alg=LGG09(δ=0.005, vars=(25), n=49))

plot(sol_BLDC01_dense, vars=(0, 25), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0, xlab="t", ylab="x25")

ρ(x25e, sol_BLDC01_dense)

ρ(x25e, sol_BLDC01_dense) <= 5.1e-3 # BLDC01 - BDS01

ρ(x25e, sol_BLDC01_dense) <= 4e-3 # BLDC01 - BDU01

ρ(x25, sol_BLDF01_discrete(20.0))

ρ(x25e, sol_BLDC01_dense(20.0)) <= -0.78e-3 # BLDC01 - BDU02

sol_BLDC01_discrete = solve(prob_BLDC01, T=20.0, alg=LGG09(δ=0.01, vars=(25), n=49, approx_model=NoBloating()))

plot(sol_BLDC01_discrete, vars=(0, 25), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0, xlab="t", ylab="x25")

ρ(x25e, sol_BLDC01_discrete)

ρ(x25e, sol_BLDC01_discrete) <= 5.1e-3 # BLDC01 - BDS01

ρ(x25e, sol_BLDC01_discrete)

ρ(x25e, sol_BLDC01_discrete) <= 4e-3 # BLDC01 - BDU01

ρ(x25e, sol_BLDC01_discrete(20.0))

ρ(x25e, sol_BLDC01_discrete(20.0)) <= -0.78e-3 # BLDC01 - BDU02


text:
using Nullables

#############
# 2D scenes #
#############

struct Point
    x::Float64
    y::Float64
end

Point(arr::Vector{Float64}) = Point(arr[1], arr[2])

array(point::Point) = [point.x, point.y]

function dist(a::Point, b::Point)
    dx = a.x - b.x
    dy = a.y - b.y
    sqrt(dx * dx + dy * dy)
end

Base.:+(a::Point, b::Point) = Point(a.x + b.x, a.y + b.y)
Base.:*(a::Real, b::Point) = Point(b.x * a, b.y * a)

abstract type Obstacle end

"""
A 2D canvas with a set of obstacles
"""
mutable struct Scene
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
    obstacles::Array{Obstacle,1}
end

function Scene(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
    Scene(xmin, xmax, ymin, ymax, Obstacle[])
end

function add!(scene::Scene, object::Obstacle)
    push!(scene.obstacles, object)
end

function ccw(a::Point, b::Point, c::Point)
    (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y)
end

function line_intersects_line(a1::Point, a2::Point, b1::Point, b2::Point)
        # http://algs4.cs.princeton.edu/91primitives/
        if ccw(a1, a2, b1) * ccw(a1, a2, b2) > 0
            return false
        end
        if ccw(b1, b2, a1) * ccw(b1, b2, a2) > 0
            return false
        end
        return true
end

function line_of_site(scene::Scene, a::Point, b::Point)
    for obstacle in scene.obstacles
        if intersects_path(obstacle, a, b)
            return false
        end
    end
    return true
end

"""
An arbitrary polygon
"""
struct Polygon <: Obstacle
    vertices::Array{Point,1}
end

function intersects_path(poly::Polygon, path_start::Point, path_end::Point)
    n = length(poly.vertices)
    for start_vertex_idx=1:n
        end_vertex_idx = start_vertex_idx % n + 1 # loop over to 1
        v1 = poly.vertices[start_vertex_idx]
        v2 = poly.vertices[end_vertex_idx]
        if line_intersects_line(v1, v2, path_start, path_end)
            return true
        end
    end
    return false
end

"""
A line with thickness
"""
struct Wall <: Obstacle
    start::Point
    orientation::Int # x is 1, y is 2
    length::Float64
    thickness::Float64
    height::Float64
    poly::Polygon
    function Wall(start::Point, orientation::Int, length::Float64,
                  thickness::Float64, height::Float64)
        if orientation != 1 && orientation != 2
            error("orientation must be either 1 (x) or 2 (y)")
        end
        vertices = Vector{Point}(undef, 4)
        vertices[1] = start
        dx = orientation == 1 ? length : thickness
        dy = orientation == 2 ? length : thickness
        vertices[2] = Point(start.x + dx, start.y)
        vertices[3] = Point(start.x + dx, start.y + dy)
        vertices[4] = Point(start.x, start.y + dy)
        poly = Polygon(vertices)
        new(start, orientation, length, thickness, height, poly)
    end
end

function intersects_path(wall::Wall, path_start::Point, path_end::Point)
    intersects_path(wall.poly, path_start, path_end)
end

"""
A square
"""
struct Tree <: Obstacle
    center::Point
    size::Float64
    poly::Polygon
    function Tree(center::Point; size::Float64=10.)
        vertices = Vector{Point}(undef, 4)
        vertices[1] = Point(center.x - size/2, center.y - size/2)
        vertices[2] = Point(center.x + size/2, center.y - size/2)
        vertices[3] = Point(center.x + size/2, center.y + size/2)
        vertices[4] = Point(center.x - size/2, center.y + size/2)
        poly = Polygon(vertices)
        new(center, size, poly)
    end
end

function intersects_path(tree::Tree, path_start::Point, path_end::Point)
    intersects_path(tree.poly, path_start, path_end)
end

##########
# render #
##########

using PyPlot

const patches = PyPlot.matplotlib[:patches]

function render(polygon::Polygon, ax)
    num_sides = length(polygon.vertices)
    vertices = Matrix{Float64}(undef, num_sides, 2)
    for (i, pt) in enumerate(polygon.vertices)
        vertices[i,:] = [pt.x, pt.y]
    end
    poly = patches[:Polygon](vertices, true, facecolor="black")
    ax[:add_patch](poly)
end

render(wall::Wall, ax) = render(wall.poly, ax)
render(tree::Tree, ax) = render(tree.poly, ax)

function render(scene::Scene, ax)
    ax[:set_xlim]((scene.xmin, scene.xmax))
    ax[:set_ylim]((scene.ymin, scene.ymax))
    for obstacle in scene.obstacles
        render(obstacle, ax)
    end
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

# This tests are supposed to test queries on the circuits
@testset "Logistic Circuit Class Conditional" begin
    # Uses a Logistic Circuit with 4 variables, and tests 3 of the configurations to
    # match with python version.

    EPS = 1e-7;
    my_opts = (max_factors= 2,
            compact⋀=false,
            compact⋁=false)

    logistic_circuit = zoo_lc("little_4var.circuit", 2);
    @test logistic_circuit isa Vector{<:LogisticΔNode};

    flow_circuit = FlowΔ(logistic_circuit, 16, Float64, my_opts)
    @test flow_circuit isa Vector{<:FlowΔNode};

    # Step 1. Check Probabilities for 3 samples
    data = XData(Bool.([0 0 0 0; 0 1 1 0; 0 0 1 1]));

    true_prob = [3.43147972 4.66740416;
                4.27595352 2.83503504;
                3.67415087 4.93793472]

    CLASSES = 2
    calc_prob = class_conditional_likelihood_per_instance(flow_circuit, CLASSES, data)

    for i = 1:3
        for j = 1:2
            @test true_prob[i,j] ≈ calc_prob[i,j] atol= EPS;
        end
    end

    # 2. Testing different API
    fc2, calc_prob2 = class_conditional_likelihood_per_instance(logistic_circuit, CLASSES, data)
    for i = 1:3
        for j = 1:2
            @test true_prob[i,j] ≈ calc_prob2[i,j] atol= EPS;
        end
    end

end


text:
println("Loading Packges for Cytof5 test...")
using Cytof5
using Test
using LinearAlgebra
using Distributions
using RCall
import Random
println("Starting Tests for Cytof5 test...")

include("Util_tests.jl")
include("VB_tests.jl")
include("Model_tests.jl")
include("MCMC_tests.jl")


text:
module ImageDipoles

using LinearAlgebra
using StaticArrays

export SV3

const SV3 = SVector{3,Float64};

"Calculates the electric potential at point r of a charge q₁ and permanent dipole d₁ at position r₁.
   Distances are measured in nm, charge in units of e, potential in V.
   The full expression is:
     ϕ = (q₁ + d₁⋅rᵣ/Rᵣ²) / (4πϵ₀ Rᵣ)
   where rᵣ = r₁-r and Rᵣ = |rᵣ|
   In our units, 1/(4πϵ₀) = 1.4399644149416115 V nm/e"
@inline function ϕ_chargedipole(r₁,q₁,d₁,r)
    rᵣ = r - r₁
    Rᵣ2 = sum(abs2,rᵣ) + 1e-200
    ϕRᵣ = q₁ + d₁⋅rᵣ/Rᵣ2
    return 1.4399644149416115 * ϕRᵣ / √Rᵣ2
end

"Calculates the electric field at point r of a charge q₁ and permanent dipole d₁ at position r₁.
   Distances are measured in nm, charge in units of e, potential in V.
   The full expression is:
     E = ((q₁ + 3d₁⋅rᵣ/Rᵣ²) rᵣ - d₁) / (4πϵ₀ Rᵣ³)
   where rᵣ = r₁-r and Rᵣ = |rᵣ|
   In our units, 1/(4πϵ₀) = 1.4399644149416115 V nm/e"
@inline function E_chargedipole(r₁,q₁,d₁,r)
    rᵣ = r - r₁
    Rᵣ2 = sum(abs2,rᵣ) + 1e-200
    ERᵣ3 = (q₁ + 3d₁⋅rᵣ/Rᵣ2) * rᵣ - d₁
    return 1.4399644149416115 * ERᵣ3 / (Rᵣ2*√Rᵣ2)
end

@inline E_dipole(r1,d1,r) = E_chargedipole(r1,0.,d1,r)
@inline V_chargedipole_dipole(r1,q1,d1,r2,d2) = -dot(E_chargedipole(r1,q1,d1,r2),d2);
@inline V_dipole_dipole(r1,d1,r2,d2) = -dot(E_dipole(r1,d1,r2),d2);

function V_dipoles_images(rs, ds, rims, qims, dims)
    N = length(rs)
    M = length(rims)
    length(rs)==length(ds)==N || error("Must have length(rs)==length(ds)==N!")
    length(rims)==length(dims)==length(qims)==M || error("Must have length(rims)==length(dims)==length(qims)==M!")
    # array of dipole-dipole interactions
    Vdd = zeros(N,N)
    # array of dipole-image interactions (include factor 1/2)
    Vid = zeros(M,N)
    @views for ii = 1:N
        r1 = rs[ii]
        d1 = ds[ii]
        for jj = ii+1:N
            Vdd[jj,ii] = V_dipole_dipole(r1,d1,rs[jj],ds[jj])
            Vdd[ii,jj] = Vdd[jj,ii]
        end
        for jj = 1:M
            Vid[jj,ii] = 0.5 * V_chargedipole_dipole(rims[jj],qims[jj],dims[jj],r1,d1)
        end
    end
    return Vdd, Vid
end

V_dipoles(rs,ds) = V_dipoles_images(rs, ds, (), (), ())[1];

mutable struct ChargesDipoles
    rs::Vector{SV3}
    qs::Vector{Float64}
    ds::Vector{SV3}
end
ϕfun(cd::ChargesDipoles) = r -> sum(ϕ_chargedipole.(cd.rs,cd.qs,cd.ds,(r,)));
Efun(cd::ChargesDipoles) = r -> sum(E_chargedipole.(cd.rs,cd.qs,cd.ds,(r,)));
ϕfun(cd1::ChargesDipoles,cd2::ChargesDipoles) = r -> sum(ϕ_chargedipole.(cd1.rs,cd1.qs,cd1.ds,(r,))) + sum(ϕ_chargedipole.(cd2.rs,cd2.qs,cd2.ds,(r,)));
Efun(cd1::ChargesDipoles,cd2::ChargesDipoles) = r -> sum(E_chargedipole.(cd1.rs,cd1.qs,cd1.ds,(r,))) + sum(E_chargedipole.(cd2.rs,cd2.qs,cd2.ds,(r,)));

@inline function plate_dipole_image(r::T,d::T,zplate) where T
    zrel = r[3] - zplate
    # position is reflected in z
    rp = T( r[1],  r[2], r[3]-2zrel)
    # dipole moment is reflected in x and y
    dp = T(-d[1], -d[2], d[3])
    return rp, dp
end

@inline function sphere_chargedipole_image(r,q,d,R)
    ar2 = sum(abs2,r)
    ar = sqrt(ar2)
    ar3 = ar2*ar
    rp = r * R^2/ar2
    rdotd = dot(r,d)
    dp = R^3/ar3 * (2r*rdotd/ar2 - d)
    qp = R*(rdotd/ar3 - q/ar)
    return rp, qp, dp
end
@inline function sphere_chargedipole_image(r,q,d,R,R0)
    rp, qp, dp = sphere_chargedipole_image(r-R0,q,d,R)
    return rp + R0, qp, dp
end

function twosphere_chargedipole_images(r0,q0,d0,Rs,R0s,tol=1e-6)
    rrs, qqs, dds = [r0], [q0], [d0]
    tol = tol*(norm(d0)+abs(q0))
    @views for ii=0:1
        # do not get images starting in sphere if initial position is inside sphere
        norm(r0-R0s[ii%2+1])<=Rs[ii%2+1] && continue
        rs, qs, ds = [r0], [q0], [d0]
        while norm(ds[end])>tol || abs(qs[end])>tol
            push!.((rs, qs, ds), sphere_chargedipole_image(rs[end],qs[end],ds[end],Rs[ii%2+1],R0s[ii%2+1]))
            ii += 1
        end
        # take 2:end-1 because the first one is the original dipole and the last one is below the cutoff
        rrs = [rrs;rs[2:end-1]]
        qqs = [qqs;qs[2:end-1]]
        dds = [dds;ds[2:end-1]]
    end
    return ChargesDipoles(rrs,qqs,dds)
end

# stateful iterator over image chargedipoles
mutable struct Twosphere_Chargedipole_ImageIter
    Rs::Tuple{Float64,Float64}
    R0s::Tuple{SV3,SV3}
    r::SV3
    q::Float64
    d::SV3
    nextsphere::Int
end
function Base.iterate(iter::Twosphere_Chargedipole_ImageIter,state=nothing)
    iter.r, iter.q, iter.d = sphere_chargedipole_image(iter.r,iter.q,iter.d,iter.Rs[iter.nextsphere],iter.R0s[iter.nextsphere])
    # 1 -> 2, 2 -> 1
    iter.nextsphere = 3 - iter.nextsphere
    return (iter.r, iter.q, iter.d), nothing
end

# stateful iterator over image dipoles of two parallel plates at zs
mutable struct Twoplate_Dipole_ImageIter
    zs::Tuple{Float64,Float64}
    r::SV3
    d::SV3
    nextplate::Int
end
function Base.iterate(iter::Twoplate_Dipole_ImageIter,state=nothing)
    iter.r, iter.d = plate_dipole_image(iter.r,iter.d,iter.zs[iter.nextplate])
    # 1 -> 2, 2 -> 1
    iter.nextplate = 3 - iter.nextplate
    return (iter.r, iter.d), nothing
end

function twoplate_V_dipole_dipole(zs,r1,d1,r2,d2,tol=1e-6)
    # interaction of dipole d1 at r1 with dipole r2 d2 and its images
    Vdd = V_dipole_dipole(r2,d2,r1,d1)
    Vdi = 0.
    Vmax = r1==r2 ? 0. : abs(Vdd)
    for iplate = 1:2
        ii = 0
        ismall = 0
        for (r,d) in Twoplate_Dipole_ImageIter(zs,r2,d2,iplate)
            ii += 1
            Vc = V_dipole_dipole(r,d,r1,d1)
            Vdi += Vc
            Vmax = max(Vmax,abs(Vc))
            ismall = abs(Vc)/abs(Vmax) < tol ? ismall+1 : 0
            ismall > 2 && break
        end
        #println("did isphere = $isphere, ii = $ii")
    end
    Vdd, Vdi
end

function twosphere_V_dipole_dipole(Rs,R0s,r1,d1,r2,d2,tol=1e-6)
    # interaction of dipole d1 at r1 with dipole r2 d2 and its images
    Vdd = V_dipole_dipole(r2,d2,r1,d1)
    Vdi = 0.
    Vmax = r1==r2 ? 0. : abs(Vdd)
    for isphere = 1:2
        ii = 0
        ismall = 0
        for (r,q,d) in Twosphere_Chargedipole_ImageIter(Rs,R0s,r2,0.,d2,isphere)
            ii += 1
            Vc = V_chargedipole_dipole(r,q,d,r1,d1)
            Vdi += Vc
            Vmax = max(Vmax,abs(Vc))
            ismall = abs(Vc)/abs(Vmax) < tol ? ismall+1 : 0
            ismall > 2 && break
        end
        #println("did isphere = $isphere, ii = $ii")
    end
    Vdd, Vdi
end

function environ_n_dipole_interactions(Vdidfun, rs, ds)
    N = length(rs)
    N == length(ds) || error("length(rs) = $(length(rs)) must be equal to length(ds) = $(length(ds))!")
    Vdd = Array{Float64}(undef,N,N)
    Vid = similar(Vdd)
    for ii = 1:length(rs)
        for jj = 1:ii
            Vdd[ii,jj], Vid[ii,jj] = Vdidfun(rs[ii],ds[ii],rs[jj],ds[jj])
            Vdd[jj,ii], Vid[jj,ii] = Vdd[ii,jj], Vid[ii,jj]
        end
        Vdd[ii,ii] = 0.
    end
    return Vdd,Vid
end

function twoplate_n_dipole_interactions(zs, rs, ds, tol=1e-6)
    Vdidfun(r1,d1,r2,d2) = twoplate_V_dipole_dipole(zs,r1,d1,r2,d2,tol)
    environ_n_dipole_interactions(Vdidfun, rs, ds)
end

function twosphere_n_dipole_interactions(Rs, R0s, rs, ds, tol=1e-6)
    Vdidfun(r1,d1,r2,d2) = twosphere_V_dipole_dipole(Rs,R0s,r1,d1,r2,d2,tol)
    environ_n_dipole_interactions(Vdidfun, rs, ds)
end

end # module


text:
# # GCM EKE Spectra

# ## Load modules
using Plots
using NCDatasets
using Statistics: mean
using DelimitedFiles
using PrettyTables
using PaddedViews
using Dates
using Printf

# ## Setup run-time enviromnent
ENV["GKSwstype"] = "100"

# ## Extract data

# Specify location of .nc files
CLIMA_NETCDF = "../netcdf/";

# Get all netcdf file names in the CLIMA_NETCDF directory and select the last one
fnames = filter(x -> occursin(".nc", x), readdir( CLIMA_NETCDF ) );
filename = "$CLIMA_NETCDF"fnames[end]

# Get the spectra data
ds = NCDataset(filename, "r")
spec1 = ds["spectrum_1d"][:] # m, lat, lev, time,
spec2 = ds["spectrum_2d"][:] # m_t, n, lev, time,
time = ds["time"][:]
lev = ds["level"][:]
lat = ds["lat"][:]
m = ds["m"][:]
n = ds["n"][:]
m_t = ds["m_t"][:]
close(ds)

# ## Plot the power spectrum
t_i = length(time)
l_1 = 1
plot(n, sum(spec2[:,:,l_1,t_i], dims = 1)[1,:], title = "EKE power spec at time $(time[t_i])", xlabel = "n", ylabel = "(m2/s2)" )


text:
macro v(val)
    esc(:(view([$val],1)))
end

@testset "_gettarget" begin
    @test_throws UndefVarError _gettarget
    @test_throws UndefVarError _gettarget(X)
    @test typeof(MLDataPattern._gettarget) <: Function

    @testset "Any" begin
        @test_throws MethodError MLDataPattern._gettarget(:a)
        @test_throws MethodError MLDataPattern._gettarget("test")
        @test_throws MethodError MLDataPattern._gettarget(3.0)
        @test_throws MethodError MLDataPattern._gettarget(2)
        @test_throws MethodError MLDataPattern._gettarget(X)
        @test_throws MethodError MLDataPattern._gettarget(uppercase, "test")
        @test @inferred(MLDataPattern._gettarget(identity, @v(1))) === 1
        @test @inferred(MLDataPattern._gettarget(x->x.+1,[1,2])) == [2,3]
        @test @inferred(MLDataPattern._gettarget(identity, X)) === X
        @test @inferred(MLDataPattern._gettarget(identity, y)) === y
        @test @inferred(MLDataPattern._gettarget(identity, yv)) === yv
        @test @inferred(MLDataPattern._gettarget(identity, ys)) === ys
        # TODO 0.6: test that first parameter need not be
        #           "typeof(...) <: Function", just callable.
    end

    @testset "Array" begin
        @test @inferred(MLDataPattern._gettarget(identity, y)) === y
        @test @inferred(MLDataPattern._gettarget(identity, X)) === X
        @test @inferred(MLDataPattern._gettarget(x->x[1]>x[2], [2,1]))
        @test @inferred(MLDataPattern._gettarget(argmax, [2,3,1])) === 2
    end

    @testset "SubArray" begin
        tmp = [2,3,1]
        tmpv = view(tmp, :)
        @test @inferred(MLDataPattern._gettarget(argmax, tmpv)) === 2
        @test @inferred(MLDataPattern._gettarget(identity, yv)) === yv
    end

    @testset "AbstractSparseArray" begin
        tmp = [0,3,0,1]
        tmps = sparse(tmp)
        @test @inferred(MLDataPattern._gettarget(nnz, tmps)) === 2
        @test @inferred(MLDataPattern._gettarget(identity, ys)) === ys
    end

    @testset "DataSubset" begin
        @test_throws MethodError MLDataPattern._gettarget(DataSubset(CustomType()))
        @test_throws MethodError MLDataPattern._gettarget(identity, DataSubset(CustomObs(1)))
        @test @inferred(MLDataPattern._gettarget(identity, DataSubset(CustomType()))) == collect(1:100)
    end

    @testset "Tuple" begin
        @test_throws MethodError MLDataPattern._gettarget(("test",))
        @test_throws MethodError MLDataPattern._gettarget((1,))
        @test_throws MethodError MLDataPattern._gettarget((1,2.0))
        @test_throws MethodError MLDataPattern._gettarget(identity, ("test",))
        @test @inferred(MLDataPattern._gettarget(identity, (@v("test"),))) == "test"
        @test @inferred(MLDataPattern._gettarget(uppercase, (@v("test"),))) == "TEST"
        @test @inferred(MLDataPattern._gettarget(identity, (@v(1),))) === 1
        @test @inferred(MLDataPattern._gettarget(identity, (@v(1),@v(2.0)))) === 2.0
        @test @inferred(MLDataPattern._gettarget(x->x+2,(@v(1),@v(2.0)))) === 4.0
        @test @inferred(MLDataPattern._gettarget(identity, (@v(1),@v(2.0),@v(:a)))) === :a

        # nobs not checked
        @test @inferred(MLDataPattern._gettarget(identity, (X,@v(2)))) === 2
        @test @inferred(MLDataPattern._gettarget(identity, (X,Y))) === Y

        @test @inferred(MLDataPattern._gettarget(identity, (y,))) === y
        @test @inferred(MLDataPattern._gettarget(identity, (ys,))) === ys
        @test @inferred(MLDataPattern._gettarget(identity, (yv,))) === yv
        @test @inferred(MLDataPattern._gettarget(identity, (ys,yv))) === yv
        @test @inferred(MLDataPattern._gettarget(identity, (XX,X,y))) === y
        @test @inferred(MLDataPattern._gettarget(identity, (XX,X,Y))) === Y

        # nested tuples
        @test @inferred(MLDataPattern._gettarget(identity, (y,(y,Y)))) === (y,Y)
        @test @inferred(MLDataPattern._gettarget(identity, (y,(yv,)))) === (yv,)
        @test @inferred(MLDataPattern._gettarget(x->map(round,x),(1,(@v(3.1),@v(4.9))))) === (3.,5.)
    end

    @testset "Custom Types" begin
        @test_throws MethodError MLDataPattern._gettarget(identity, EmptyType())
        @test_throws MethodError MLDataPattern._gettarget(identity, (EmptyType(),))
        @test @inferred(MLDataPattern._gettarget(identity, CustomType())) == collect(1:100)
        @test @inferred(MLDataPattern._gettarget(identity, (X,CustomType()))) == collect(1:100)
        @test @inferred(MLDataPattern._gettarget(identity, (y,DataSubset(CustomType())))) == collect(1:100)

        @test @inferred(MLDataPattern._gettarget("test", DataSubset(CustomStorage(),1))) == "test - obs 1"
        @test @inferred(MLDataPattern._gettarget(identity, DataSubset(CustomStorage(),2))) == "obs 2"
        @test @inferred(MLDataPattern._gettarget(identity, (1.0,DataSubset(CustomStorage(),2)))) == "obs 2"
        @test @inferred(MLDataPattern._gettarget("a", (1.0,DataSubset(CustomStorage(),2)))) == "a - obs 2"
    end
end

println("<HEARTBEAT>")

@testset "gettarget" begin
    @test_throws UndefVarError gettarget
    @test_throws UndefVarError gettarget(X)
    @test typeof(MLDataPattern.gettarget) <: Function

    @testset "Any" begin
        @test @inferred(MLDataPattern.gettarget(:a)) == :a
        @test @inferred(MLDataPattern.gettarget("test")) == "test"
        @test @inferred(MLDataPattern.gettarget(uppercase, "test")) == "TEST"
        @test @inferred(MLDataPattern.gettarget(3.0)) === 3.0
        @test @inferred(MLDataPattern.gettarget(2)) === 2
        @test @inferred(MLDataPattern.gettarget(x->x+1,2)) === 3
    end

    @testset "Array" begin
        @test @inferred(MLDataPattern.gettarget(y)) === y
        @test @inferred(MLDataPattern.gettarget(identity, y)) === y
        @test @inferred(MLDataPattern.gettarget(X)) === X
        @test @inferred(MLDataPattern.gettarget(identity, X)) === X
        @test @inferred(MLDataPattern.gettarget(x->x[1]>x[2], [2,1]))
        @test @inferred(MLDataPattern.gettarget(argmax, [2,3,1])) === 2
    end

    @testset "SubArray" begin
        tmp = [2,3,1]
        tmpv = view(tmp, :)
        @test @inferred(MLDataPattern.gettarget(argmax, tmpv)) === 2
        @test @inferred(MLDataPattern.gettarget(yv)) === yv
        @test @inferred(MLDataPattern.gettarget(identity, yv)) === yv
    end

    @testset "AbstractSparseArray" begin
        tmp = [0,3,0,1]
        tmps = sparse(tmp)
        @test @inferred(MLDataPattern.gettarget(nnz, tmps)) === 2
        @test @inferred(MLDataPattern.gettarget(ys)) === ys
        @test @inferred(MLDataPattern.gettarget(identity, ys)) === ys
    end

    @testset "DataSubset" begin
        @test_throws MethodError MLDataPattern.gettarget(identity, DataSubset(CustomObs(1)))
        @test @inferred(MLDataPattern.gettarget(DataSubset(CustomType()))) == DataSubset(CustomType())
        @test @inferred(MLDataPattern.gettarget(identity, DataSubset(CustomType()))) == DataSubset(CustomType())
    end

    @testset "Tuple" begin
        @test @inferred(MLDataPattern.gettarget(("test",))) == ("test",)
        @test_throws MethodError MLDataPattern.gettarget(uppercase, ("test",))
        @test @inferred(MLDataPattern.gettarget(x->map(uppercase,x), ("test",))) == ("TEST",)
        @test @inferred(MLDataPattern.gettarget((1,2.0))) === (1,2.0)
        @test_throws MethodError MLDataPattern.gettarget(x->x+2,(1,2.0))
        @test @inferred(MLDataPattern.gettarget(x->map(y->y+2,x),(1,2.0))) === (3,4.0)
        @test @inferred(MLDataPattern.gettarget((1,2.0,:a))) === (1,2.0,:a)
        @test @inferred(MLDataPattern.gettarget((1,(2.0,:a)))) === (1,(2.0,:a))
        @test @inferred(MLDataPattern.gettarget((y,))) === (y,)
        @test @inferred(MLDataPattern.gettarget((Y,))) === (Y,)
        @test @inferred(MLDataPattern.gettarget((yv,))) == (yv,)
        @test @inferred(MLDataPattern.gettarget((y,(y,Y)))) === (y,(y,Y))
        @test_throws MethodError MLDataPattern.gettarget(x->map(round,x),(1,(3.1,4.9)))
    end

    @testset "Custom Types" begin
        @test @inferred(MLDataPattern.gettarget(EmptyType())) === EmptyType()
        @test @inferred(MLDataPattern.gettarget(identity, EmptyType())) === EmptyType()
        @test @inferred(MLDataPattern.gettarget(CustomType())) === CustomType()
        @test @inferred(MLDataPattern.gettarget(x->x,CustomType())) === CustomType()
        @test @inferred(MLDataPattern.gettarget("a",CustomObs(4))) == "a - obs 4"
        @test @inferred(MLDataPattern.gettarget(CustomObs(5))) == "obs 5"

        @test @inferred(MLDataPattern.gettarget(DataSubset(CustomStorage(),2))) === DataSubset(CustomStorage(),2)
        @test @inferred(MLDataPattern.gettarget(identity, DataSubset(CustomStorage(),2))) == DataSubset(CustomStorage(),2)

        @test @inferred(MLDataPattern.gettarget((1:100,DataSubset(CustomType())))) == (1:100,DataSubset(CustomType()))
        @test_throws MethodError MLDataPattern.gettarget("a",(CustomObs(2),CustomObs(2)))
    end
end

println("<HEARTBEAT>")

@testset "targets" begin
    @test typeof(targets) <: Function
    @test targets === MLDataPattern.targets

    @testset "questionable results to unusual parameters" begin
        @test_throws MethodError targets(x->x+1, 1)
        @test_throws MethodError targets(uppercase, "test")
        @test_throws MethodError targets((1,1:3))
        @test_throws MethodError targets(x->x, (1,1:3))
        @test_throws MethodError targets((1:3,1))
        @test_throws MethodError targets(3.)
        @test_throws MethodError targets(identity, 3.)
    end

    @testset "nobs mismatch" begin
        @test_throws DimensionMismatch targets((1:3,1:4))
        @test_throws DimensionMismatch targets((1:4,1:3))
        @test_throws DimensionMismatch targets((X,Y), ObsDim.First())
        @test_throws DimensionMismatch targets((X,Yt))
        @test_throws DimensionMismatch targets((y,Yt))
        @test_throws DimensionMismatch targets((X,y), (ObsDim.First(),ObsDim.Last()))
        @test_throws DimensionMismatch targets((X,y), (ObsDim.First(),))
        @test_throws DimensionMismatch targets(identity, (X,y), (ObsDim.First(),))
    end

    @testset "Arrays" begin
        @test @inferred(targets(y)) === y
        @test @inferred(targets(x->x=="setosa", y)) == vcat(trues(50), falses(100))
        @test @inferred(targets(Y)) === Y
        @test targets(Y, obsdim=:last) === Y
        @test @inferred(targets(Y, ObsDim.Last())) === Y
        @test @inferred(targets(x->x[1]==x[2], Y)) == trues(nobs(Y))
        @test targets(x->x[1]==x[2], Yt, obsdim=1) == trues(nobs(Y))
        @test @inferred(targets(x->x[1]==x[2], Yt, ObsDim.First())) == trues(nobs(Y))
        @test @inferred(targets(yv)) !== yv
        @test @inferred(targets(yv)) == yv
        @test @inferred(targets(ys)) === ys
        # TODO: inference broken in test modus but not otherwise :-(
        @test targets(uppercase, ["a","b","c"]) == ["A","B","C"]
    end

    @testset "Tuples" begin
        @test @inferred(targets((X,yv))) !== yv
        @test @inferred(targets((X,yv))) == yv
        @test @inferred(targets((yv,))) !== yv
        @test @inferred(targets((yv,))) == yv
        @test @inferred(targets((ys,))) === ys
        @test @inferred(targets((Y,))) === Y

        @test targets(x->map(uppercase,x), ["a" "b" "c"; "d" "e" "f"], obsdim=2) == [["A", "D"], ["B", "E"], ["C", "F"]]
        @test targets(x->map(uppercase,x), ["a" "b" "c"; "d" "e" "f"], ObsDim.Last()) == [["A", "D"], ["B", "E"], ["C", "F"]]
        @test targets((y,Y),obsdim=:last) === Y
        @test targets((y,Yt),obsdim=:first) === Yt
        @test targets((Yt,yt),obsdim=(:first,:last)) === yt
        @test targets((Y,Yt),obsdim=(:last,:first)) === Yt

        @test @inferred(targets((y,Y),ObsDim.Last())) === Y
        @test @inferred(targets((y,Yt),ObsDim.First())) === Yt
        @test @inferred(targets((Yt,yt),(ObsDim.First(),ObsDim.Last()))) === yt
        @test @inferred(targets((Y,Yt),(ObsDim.Last(),ObsDim.First()))) === Yt

        @test @inferred(targets(identity,(y,Y),ObsDim.Last())) === Y
        @test @inferred(targets(identity,(y,Yt),ObsDim.First())) === Yt
        @test @inferred(targets(identity,(Yt,yt),(ObsDim.First(),ObsDim.Last()))) === yt
        @test @inferred(targets(identity,(Y,Yt),(ObsDim.Last(),ObsDim.First()))) === Yt

        @test @inferred(targets(([:a,:b,:c],1:3))) == [1,2,3]
        @test @inferred(targets(x->x^2,([:a,:b,:c],1:3))) == [1,4,9]
        # nested tuples
        @test @inferred(targets(([:a,:b,:c],(1:3,4:6)))) === (1:3,4:6)
        @test @inferred(targets(x->x[1]*x[2],([:a,:b,:c],(1:3,4:6)))) == [4,10,18]
        @test @inferred(targets(DataSubset(([:a,:b,:c],(1:3,4:6)), [3,2,1]))) == ([3,2,1],[6,5,4])
        @test @inferred(targets(datasubset(([:a,:b,:c],(1:3,4:6)), [3,2,1]))) == ([3,2,1],[6,5,4])
        @test @inferred(targets(DataSubset(([:a,:b,:c],(1:3,4:6))))) === (1:3,4:6)
        @test @inferred(targets(datasubset(([:a,:b,:c],(1:3,4:6))))) == ([1,2,3],[4,5,6])
    end

    @testset "BatchView" begin
        ft = rand(6)
        str = ["a", "b", "b", "b", "b", "a"]
        bv = batchview(str)
        # TODO: inference broken in test modus but not otherwise :-(
        @test targets(bv) == [["a", "b"], ["b", "b"], ["b", "a"]]
        @test targets(uppercase, bv) == [["A", "B"], ["B", "B"], ["B", "A"]]
        bv = batchview((ft,str))
        # TODO: inference broken in test modus but not otherwise :-(
        @test targets(bv) == [["a", "b"], ["b", "b"], ["b", "a"]]
        @test targets(uppercase, bv) == [["A", "B"], ["B", "B"], ["B", "A"]]
    end

    @testset "ObsView" begin
        ft = rand(6)
        str = ["a", "b", "b", "b", "b", "a"]
        ov = obsview(str)
        # TODO: inference broken in test modus but not otherwise :-(
        @test targets(ov) == ["a", "b", "b", "b", "b", "a"]
        @test targets(uppercase, ov) == ["A", "B", "B", "B", "B", "A"]
        ov = obsview((ft,str))
        # TODO: inference broken in test modus but not otherwise :-(
        @test targets(ov) == ["a", "b", "b", "b", "b", "a"]
        @test targets(uppercase, ov) == ["A", "B", "B", "B", "B", "A"]
    end

    @testset "CustomType" begin
        @test @inferred(targets(CustomType())) == "batch 1:100"
        @test @inferred(targets(identity,CustomType())) == "batch 1:100"
        @test @inferred(targets(datasubset(CustomType()))) == "batch 1:100"
        @test @inferred(targets(datasubset(CustomType(),2:5))) == "batch 2:5"
        @test @inferred(targets(datasubset(CustomType(),7))) == "obs 7"
        @test @inferred(targets(x->x,CustomType())) == collect(1:100)
        @test @inferred(targets(x->x,datasubset(CustomType()))) == collect(1:100)
        @test @inferred(targets(x->x,datasubset(CustomType(),2:5))) == collect(2:5)
        @test @inferred(targets(x->x,datasubset(CustomType(),3))) == [3]
    end

    @testset "CustomStorage" begin
        @test @inferred(targets(CustomStorage())) == ["obs 1", "obs 2"]
        @test @inferred(targets(identity,CustomStorage())) == ["obs 1", "obs 2"]
        @test @inferred(targets(x->x,CustomStorage())) == ["obs 1", "obs 2"]
        @test @inferred(targets(datasubset(CustomStorage()))) == ["obs 1", "obs 2"]
        @test @inferred(targets(datasubset(CustomStorage(),[2,1]))) == ["obs 2", "obs 1"]
        @test @inferred(targets(datasubset(CustomStorage(),2))) == "obs 2"
        @test @inferred(targets(x->x,datasubset(CustomStorage()))) == ["obs 1", "obs 2"]
        @test @inferred(targets(x->x,datasubset(CustomStorage(),[2,1]))) == ["obs 2", "obs 1"]
        @test @inferred(targets(x->x,datasubset(CustomStorage(),2))) == ["obs 2"]
        @test @inferred(targets("a",CustomStorage())) == ["a - obs 1", "a - obs 2"]
        @test @inferred(targets("a",datasubset(CustomStorage()))) == ["a - obs 1", "a - obs 2"]
        @test @inferred(targets("a",datasubset(CustomStorage(),[2,1]))) == ["a - obs 2", "a - obs 1"]
        @test @inferred(targets("a",datasubset(CustomStorage(),2))) == ["a - obs 2"]
    end

    @testset "MetaDataStorage" begin
        @test @inferred(targets(MetaDataStorage())) == "full"
        @test @inferred(targets(identity,MetaDataStorage())) == "full"
        @test @inferred(targets(datasubset(MetaDataStorage()))) == "batch 1:3"
        @test @inferred(targets(datasubset(MetaDataStorage(),2:3))) == "batch 2:3"
        @test @inferred(targets(datasubset(MetaDataStorage(),2))) == "obs 2"
        @test_throws ObsDimTriggeredException targets(x->x,MetaDataStorage())
        @test_throws ObsDimTriggeredException targets(x->x,datasubset(MetaDataStorage()))
        @test_throws ObsDimTriggeredException targets(x->x,datasubset(MetaDataStorage(),2:3))
        @test_throws ObsDimTriggeredException targets(x->x,datasubset(MetaDataStorage(),2))
    end
end

println("<HEARTBEAT>")

@testset "eachtarget" begin
    @test typeof(eachtarget) <: Function
    @test eachtarget === MLDataPattern.eachtarget

    @testset "questionable results to unusual parameters" begin
        @test_throws MethodError eachtarget(x->x+1, 1)
        @test_throws MethodError eachtarget(uppercase, "test")
        @test_throws MethodError eachtarget((1,1:3))
        @test_throws MethodError eachtarget((1:3,1))
        @test_throws MethodError eachtarget(identity, 3.)
    end

    @testset "nobs mismatch" begin
        @test_throws DimensionMismatch eachtarget((1:3,1:4))
        @test_throws DimensionMismatch eachtarget((1:4,1:3))
        @test_throws DimensionMismatch eachtarget((X,Yt))
        @test_throws DimensionMismatch eachtarget((y,Yt))
        @test_throws DimensionMismatch eachtarget((X,y), (ObsDim.First(),))
        @test_throws DimensionMismatch eachtarget(identity, (X,y), (ObsDim.First(),))
    end

    @testset "Arrays" begin
        @test eltype(collect(@inferred(eachtarget(Y)))) <: Vector
        @test eltype(collect(@inferred(eachtarget(x->x,Y)))) <: SubArray
        @test collect(@inferred(eachtarget(y))) == y
        @test collect(@inferred(eachtarget(x->x=="setosa", y))) == vcat(trues(50), falses(100))
        @test collect(@inferred(eachtarget(Y))) == obsview(Y)
        @test collect(eachtarget(Y, obsdim=:last)) == obsview(Y)
        @test collect(@inferred(eachtarget(Y, ObsDim.Last()))) == obsview(Y)
        @test collect(@inferred(eachtarget(x->x[1]==x[2], Y))) == trues(nobs(Y))
        @test collect(eachtarget(x->x[1]==x[2], Yt, obsdim=1)) == trues(nobs(Y))
        @test collect(@inferred(eachtarget(x->x[1]==x[2], Yt, ObsDim.First()))) == trues(nobs(Y))
        @test collect(@inferred(eachtarget(yv))) == y
        @test collect(@inferred(eachtarget(ys))) == collect(ys)
        # TODO: inference broken in test modus but not otherwise :-(
        @test collect(eachtarget(uppercase, ["a","b","c"])) == ["A","B","C"]
    end

    @testset "Tuples" begin
        @test collect(@inferred(eachtarget((X,yv)))) == y
        @test collect(@inferred(eachtarget((yv,)))) == y
        @test collect(@inferred(eachtarget((ys,)))) == collect(ys)
        @test collect(@inferred(eachtarget((Y,)))) == obsview(Y)

        @test collect(eachtarget(x->map(uppercase,x), ["a" "b" "c"; "d" "e" "f"], obsdim=2)) == [["A", "D"], ["B", "E"], ["C", "F"]]
        @test collect(eachtarget(x->map(uppercase,x), ["a" "b" "c"; "d" "e" "f"], ObsDim.Last())) == [["A", "D"], ["B", "E"], ["C", "F"]]
        @test collect(eachtarget((y,Y),obsdim=:last)) == obsview(Y)
        @test collect(eachtarget((y,Yt),obsdim=:first)) == obsview(Y)
        @test collect(eachtarget((Yt,y),obsdim=(:first,:last))) == y
        @test collect(eachtarget((Y,Yt),obsdim=(:last,:first))) == obsview(Y)

        @test collect(@inferred(eachtarget((y,Y),ObsDim.Last()))) == obsview(Y)
        @test collect(@inferred(eachtarget((y,Yt),ObsDim.First()))) == obsview(Y)
        @test collect(@inferred(eachtarget((Yt,y),(ObsDim.First(),ObsDim.Last())))) == y
        @test collect(@inferred(eachtarget((Y,Yt),(ObsDim.Last(),ObsDim.First())))) == obsview(Y)

        @test collect(@inferred(eachtarget(identity,(y,Y),ObsDim.Last()))) == obsview(Y)
        @test collect(@inferred(eachtarget(identity,(y,Yt),ObsDim.First()))) == obsview(Y)
        @test collect(@inferred(eachtarget(identity,(Yt,y),(ObsDim.First(),ObsDim.Last())))) == y
        @test collect(@inferred(eachtarget(identity,(Y,Yt),(ObsDim.Last(),ObsDim.First())))) == obsview(Y)

        @test collect(@inferred(eachtarget(([:a,:b,:c],1:3)))) == [1,2,3]
        @test collect(@inferred(eachtarget(x->x^2,([:a,:b,:c],1:3)))) == [1,4,9]
        # nested tuples
        @test collect(@inferred(eachtarget(([:a,:b,:c],(1:3,4:6))))) == getobs(obsview((1:3,4:6)))
        @test collect(@inferred(eachtarget(x->x[1]*x[2],([:a,:b,:c],(1:3,4:6))))) == [4,10,18]
    end

    println("<HEARTBEAT>")

    @testset "ObsView" begin
        ft = rand(6)
        str = ["a", "b", "b", "b", "b", "a"]
        ov = obsview(str)
        @inferred eachtarget(ov)
        @inferred eachtarget(uppercase, ov)
        @test typeof(eachtarget(ov)) <: Base.Generator
        @test typeof(eachtarget(uppercase, ov)) <: Base.Generator
        @test collect(eachtarget(ov)) == ["a", "b", "b", "b", "b", "a"]
        @test collect(eachtarget(uppercase, ov)) == ["A", "B", "B", "B", "B", "A"]
        ov = obsview((ft,str))
        @inferred eachtarget(ov)
        @inferred eachtarget(uppercase, ov)
        @test typeof(eachtarget(ov)) <: Base.Generator
        @test typeof(eachtarget(uppercase, ov)) <: Base.Generator
        @test collect(eachtarget(ov)) == ["a", "b", "b", "b", "b", "a"]
        @test collect(eachtarget(uppercase, ov)) == ["A", "B", "B", "B", "B", "A"]
    end

    @testset "CustomType" begin
        @test collect(@inferred(eachtarget(CustomType()))) == collect("obs $i" for i in 1:100)
        @test collect(@inferred(eachtarget(identity,CustomType()))) == collect("obs $i" for i in 1:100)
        @test collect(@inferred(eachtarget(datasubset(CustomType())))) == collect("obs $i" for i in 1:100)
        @test collect(@inferred(eachtarget(datasubset(CustomType(),2:5)))) == collect("obs $i" for i in 2:5)
        @test collect(@inferred(eachtarget(datasubset(CustomType(),7)))) == ["obs 7"]
        @test collect(@inferred(eachtarget(x->x,CustomType()))) == collect(1:100)
        @test collect(@inferred(eachtarget(x->x,datasubset(CustomType())))) == collect(1:100)
        @test collect(@inferred(eachtarget(x->x,datasubset(CustomType(),2:5)))) == collect(2:5)
        @test collect(@inferred(eachtarget(x->x,datasubset(CustomType(),3)))) == [3]
    end

    @testset "CustomStorage" begin
        @test collect(@inferred(eachtarget(CustomStorage()))) == ["obs 1", "obs 2"]
        @test collect(@inferred(eachtarget(identity,CustomStorage()))) == ["obs 1", "obs 2"]
        @test collect(@inferred(eachtarget(x->x,CustomStorage()))) == ["obs 1", "obs 2"]
        @test collect(@inferred(eachtarget(datasubset(CustomStorage())))) == ["obs 1", "obs 2"]
        @test collect(@inferred(eachtarget(datasubset(CustomStorage(),[2,1])))) == ["obs 2", "obs 1"]
        @test collect(@inferred(eachtarget(datasubset(CustomStorage(),2)))) == ["obs 2"]
        @test collect(@inferred(eachtarget(x->x,datasubset(CustomStorage())))) == ["obs 1", "obs 2"]
        @test collect(@inferred(eachtarget(x->x,datasubset(CustomStorage(),[2,1])))) == ["obs 2", "obs 1"]
        @test collect(@inferred(eachtarget(x->x,datasubset(CustomStorage(),2)))) == ["obs 2"]
        @test collect(@inferred(eachtarget("a",CustomStorage()))) == ["a - obs 1", "a - obs 2"]
        @test collect(@inferred(eachtarget("a",datasubset(CustomStorage())))) == ["a - obs 1", "a - obs 2"]
        @test collect(@inferred(eachtarget("a",datasubset(CustomStorage(),[2,1])))) == ["a - obs 2", "a - obs 1"]
        @test collect(@inferred(eachtarget("a",datasubset(CustomStorage(),2)))) == ["a - obs 2"]
    end

    @testset "MetaDataStorage" begin
        @test collect(@inferred(eachtarget(MetaDataStorage()))) == ["obs 1", "obs 2", "obs 3"]
        @test collect(@inferred(eachtarget(identity,MetaDataStorage()))) == ["obs 1", "obs 2", "obs 3"]
        @test collect(@inferred(eachtarget(datasubset(MetaDataStorage())))) == ["obs 1", "obs 2", "obs 3"]
        @test collect(@inferred(eachtarget(datasubset(MetaDataStorage(),2:3)))) == ["obs 2", "obs 3"]
        @test collect(@inferred(eachtarget(datasubset(MetaDataStorage(),[3,2])))) == ["obs 3", "obs 2"]
        @test collect(@inferred(eachtarget(datasubset(MetaDataStorage(),2)))) == ["obs 2"]
        @test typeof(@inferred(eachtarget(x->x,MetaDataStorage()))) <: Base.Generator
        @test_throws ObsDimTriggeredException collect(eachtarget(x->x,MetaDataStorage()))
        @test_throws ObsDimTriggeredException collect(eachtarget(x->x,datasubset(MetaDataStorage())))
        @test_throws ObsDimTriggeredException collect(eachtarget(x->x,datasubset(MetaDataStorage(),2:3)))
        @test_throws ObsDimTriggeredException collect(eachtarget(x->x,datasubset(MetaDataStorage(),2)))
    end
end


text:
# Load GANTR code
include( "gantr.jl" )
# Load simulated data generation code
include( "binary_sampler.jl" )

using JLD2, LaTeXStrings

compriskDelta0pt1_p = sim_study_sampler_comprisk( 200, 0.1, 1.0, 2.2, 1.75, 1.2, 2.2, 5.3, 1.5)

# Save and load simulated data
#@save "compriskDelta0pt1_p.txt" compriskDelta0pt1_p
#@load "compriskDelta0pt1_p.txt"

# Define surial times and covariables
T, Z = compriskDelta0pt1

# No censoring
C = repeat( [1], inner=length(T))
Z_ar = [ [Z[i]] for i in 1:length(Z) ]

# Define data struct
data = Data( T, C, Z_ar )

# Define regression functions
function f_1( z::Array{Float64,1}, β::Array{Float64,1})
    return 1.0*(z[1] <= β[1])  # male white
end

function f_2( z::Array{Float64,1}, β::Array{Float64,1})
    return 1.0*(z[1] > β[1])  # male white
end

# Define model struct
model =  Model( data, 1., [ f_1, f_2 ])

# Precompile loglikelihood
model_loglikelihood( [rand(Uniform())], 1.0, 1.0, [ clayton_corm_lognorm_2dim_sampler( 0.2, 0.1) for j in 1:1000 ], model)

# Precompile MCMC
MCMC_comprisk_uniform( 3, 3, 3, 1.0, 1.0, model, 0.1, [1.0], 0.5, 0.0, 2.0)
# Run MCMC
run_uniform = MCMC_comprisk_uniform( 5000, 1, 1000, 1.0, 1.0, model, 0.1, [1.0], 0.5, quantile(Z,0.25), quantile(Z,0.75))

# Save and load MCMC run
#@save "run_uniform.txt" run_uniform
#@load "run_uniform.txt"

# Precompile MCMC
MCMC_comprisk_randwalk( 3, 1, 3, 1.0, 1.0, model, 0.1, 0.1, 0.1, [1.0], 0.5, quantile(Z)[2], quantile(Z)[4] )
# Run MCMC
run_rw = MCMC_comprisk_randwalk( 10000, 1, 1000, 1.0, 1.0, model, 0.1, 0.2, 0.1, [0.7], 0.5, quantile(Z)[2], quantile(Z)[4] )

# Save and load MCMC run
#@save "run_rw.txt" run_rw
#@load "run_rw.txt"

surv_vec_delta( [ 1. ], [1.], 1., 1., 1., 0.1, 3, 10., 2., model)

# surv_pop_1_map_unif = surv_vec_delta( [ 0.66 ], run_uniform[4], 1., 1., run_uniform[3], 0.1, 1000, 5., 50., model)
# surv_pop_2_map_unif = surv_vec_delta( [ 1.44 ], run_uniform[4], 1., 1., run_uniform[3], 0.1, 1000, 5., 50., model)
# surv_pop_1_mean_unif = surv_vec_delta( [ 0.66 ],  mean( run_uniform[2][500:5001] ), 1., 1., mean( run_uniform[1][500:5001] ), 0.1, 1000, 5., 50., model)
# surv_pop_2_mean_unif = surv_vec_delta( [ 1.44 ], mean( run_uniform_[2][500:5001] ), 1., 1., mean( run_uniform[1][500:5001] ), 0.1, 1000, 5., 50., model)

surv_pop_1_map_rw = surv_vec_delta( [ 0.66 ], run_rw_p[4], 1., 1., run_rw_p[3], 0.1, 1000, 5., 50., model)
surv_pop_2_map_rw = surv_vec_delta( [ 1.44 ], run_rw_p[4], 1., 1., run_rw_p[3], 0.1, 1000, 5., 50., model)
surv_pop_1_mean_rw = surv_vec_delta( [ 0.66 ],  mean( run_rw_p[2][500:5001] ), 1., 1., mean( run_rw_p[1][500:5001] ), 0.1, 1000, 5., 50., model)
surv_pop_2_mean_rw = surv_vec_delta( [ 1.44 ], mean( run_rw_p[2][500:5001] ), 1., 1., mean( run_rw_p[1][500:5001] ), 0.1, 1000, 5., 50., model)

using PyPlot

S_pop1 = zeros(Float64,100,251)
for i in 1:100
 S_pop1[i,:] = surv_vec_delta( [ 0.66 ],  run_rw_p[2][1001+90*i], 1., 1., run_rw_p[1][1001+90*i], 0.1, 1000, 5., 50., model)
end

S_pop2 = zeros(Float64,100,251)
for i in 1:100
 S_pop2[i,:] = surv_vec_delta( [ 1.44 ],  run_rw_p[2][1001+90*i], 1., 1., run_rw_p[1][1001+90*i], 0.1, 1000, 5., 50., model)
end

S2_pop1 = cred_band_mat(S_pop1,5)
S2_pop2 = cred_band_mat(S_pop2,5)

surv_pop_1_band_rw_u = [ findmax(S2_pop1[:,i])[1] for i in 1:251]
surv_pop_1_band_rw_d = [ findmin(S2_pop1[:,i])[1] for i in 1:251]
surv_pop_2_band_rw_u = [ findmax(S2_pop2[:,i])[1] for i in 1:251]
surv_pop_2_band_rw_d = [ findmin(S2_pop2[:,i])[1] for i in 1:251]

################################################################################

using Optim

x0 = rand(Uniform(),2)
lower = [ 0.00001, 0.00001 ]
upper = [ 0.99999, 1.99999 ]

f_opt(x::Array{Float64,1}) = - posterior(x,1.0,1.0,model,0.1,1000)

op = optimize( f_opt, lower, upper, x0, Fminbox(LBFGS()), Optim.Options(iterations = 20) )

op.minimizer

# Save and load op
#@save "op.txt" op
#@load "op.txt"


text:
"Holds positional values (not tuples) in relation to cardinal directions."
struct Direction
    val::String
end

"Returns an array of the cardinal Directions."
get_all_cardinals() = [Direction("North"), Direction("South"), Direction("East"), Direction("West")]

"Choose a random Direction from [North, South, East, West]."
function random_direction()
    rand(get_all_cardinals())
end

"Converts from a Direction value (not tuple) string to the engine's command notation."
function convert(direction::Direction)
    if direction == Direction("North") return Commands.NORTH
    elseif direction == Direction("South") return Commands.SOUTH
    elseif direction == Direction("East") return Commands.WEST
    elseif direction == Direction("West") return Commands.EAST
    elseif direction == Direction("Still") return Commands.STAY_STILL
    else throw(DomainError(direction))
    end
 end

"Returns the opposite cardinal direction given a direction."
function invert(direction::Direction)
   if direction == Direction("North") return Direction("South")
   elseif direction == Direction("South") return Direction("North")
   elseif direction == Direction("East") return Direction("West")
   elseif direction == Direction("West") return Direction("East")
   elseif direction == Direction("Still") return Direction("Still")
   elseif direction == nothing return nothing
   else throw(DomainError(direction))
   end
end

struct Position
    x::Int
    y::Int
end

Base.abs(position::Position) = Position(abs(position.x), abs(position.y))
Base.:+(pos_1::Position, pos_2::Position) = Position(pos_1.x + pos_2.x, pos_1.y + pos_2.y)
Base.:-(pos_1::Position, pos_2::Position) = Position(pos_1.x - pos_2.x, pos_1.y - pos_2.y)

"Returns a new position moved in that direction. The position considering a Direction cardinal value."
function directional_offset(position::Position, direction::Direction)
    if direction == Direction("North")
        offset = Position(0, -1)
    elseif direction == Direction("South")
        offset = Position(0, 1)
    elseif direction == Direction("East")
        offset = Position(1, 0)
    elseif direction == Direction("West")
        offset = Position(-1, 0)
    elseif direction == Direction("Still")
        offset = Position(0, 0)
    else
        throw(DomainError(direction))
    end
    position + offset
end

"Returns an array of all positions around this specific position in each cardinal direction."
function get_surrounding_cardinals(position::Position)
    [directional_offset(position, direction) for direction in get_all_cardinals()]
end

max_stars_repo_name: ewirkerman/Halite-III

max_issues_repo_name: ewirkerman/Halite-III

max_forks_repo_name: ewirkerman/Halite-III


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 47f98a52-f2b2-11ea-0ac6-9b4b74b7b15d
using Plots, PlutoUI

# ╔═╡ 5ac4989c-f2b3-11ea-3e8a-a79a4117e48c
using SymEngine

# ╔═╡ 6fe82e5a-f2b8-11ea-101d-d1ee5f9bd5b1
using ForwardDiff

# ╔═╡ 963c4910-f2b8-11ea-12bb-bb42539b1c93
using Zygote

# ╔═╡ 08cf2a80-f2b2-11ea-132d-2f26c934bbb2
md"""
# Numeric and automatic differentiation

STMO

**Michiel Stock**

## Motivation

Up to now, we confidently assumed that we would always be able to compute the derivative or gradient of any function. Despite differentiation being a relatively easy operation, it is frequenty not feasible (or desirable) to compute this by hand. *Numerical differentiation* can provide approximations of th derivate or gradient at a particular point. *Automatic differentiation* directly manipulates the computational graph to generate a function that computes the (exact) derivate. Such methods have advanced greatly in the last years and it is no exageration that their easy use in popular software libraries such as TenserFlow and PyTorch are a cornerstone of deep learning and other machine learning and scientific computing fields.

"""

# ╔═╡ 5154a672-f2b2-11ea-1280-3ba5dadff7f9
md"""
## Definition of a derivative

$$\frac{\text{d}f(x)}{\text{d}x} = f'(x) = \lim _{h\to 0}{\frac {f(x+h)-f(x)}{h}}.$$

Derivation is in essence a mechanical process, following the rules below.

![](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/derivatives.jpeg)

When we work with function of several variables, we use *partial derivatives* (e.g. $\frac{\partial f(x, y)}{\partial x}$), indicating we keep all variables but $x$ fixed.
"""

# ╔═╡ f53fbf24-f2b2-11ea-2f59-d943884fee9d
md"""
Our running example:

$$f(x) = \log x + \frac{\sin x}{x}$$
"""

# ╔═╡ e6a55456-fe6f-11ea-328c-47fedb388194
md"**Assignment**: implement this function."

# ╔═╡ 2be8656c-f2b3-11ea-2244-5740c807deff
f(x) = log(x)+ sin(x)/x

# ╔═╡ 2e00ef04-f2b3-11ea-1960-8136321c71db
plot(f, 1, 5, label="f(x)")

# ╔═╡ f7cd052e-f2bb-11ea-239a-47966e96e909
md"We will evaluate this function in `a=2`."

# ╔═╡ 7244eba2-f2b3-11ea-045f-ffd208852eda
a = 2.0

# ╔═╡ 401df7a4-f2b3-11ea-12a1-e9b46bcfcf1b
md"""
## Symbolic differentiation

Computing derivatives, as you have seen in basic calculus courses.

By hand or automatically:
- Maple
- Sympy (python)
- Mathematica
- Maxima

Differentiation is *easy* compared to *integration* or *sampling*.

Advantages:
- exact derivatives!
- gives the formula for different evaluations.
    - insight in the system
    - in some cases, closed-form solution extrema by solving $\frac{\text{d}f(x)}{\text{d}x}=0$
- no hyperparameters or tweaking: just works!

Disadvantages:
- some software not flexible enough (gradients, arrays, for-loops,...)
- sometimes explosion of terms: *expression swell*
- not always numerically optimal!
"""

# ╔═╡ 5db7350a-f2b3-11ea-293d-ffc572e0dbe8
@vars x  # define variable + the @ turns it into a macro

# ╔═╡ bbe02ec4-120d-11eb-39fb-8545d56de2ad
f(x)

# ╔═╡ 64b28e4a-f2b3-11ea-0fc0-7d37b4794686
df = diff(f(x), x)

# ╔═╡ 684ff7f4-f2b3-11ea-14f7-7d3f329ecd79
df(a)

# ╔═╡ 7a6a18de-f2b3-11ea-1b3d-ef7a2b1b0434
begin
plot(f, 1, 5, label="\$f(x)\$", xlabel="\$x\$", lw=2, color=:green)
plot!(df, 1, 5, label="\$f'(x)\$", lw=2, color=:orange)
end

# ╔═╡ 913234a2-f2b3-11ea-3501-e3bfc54f91fe
md"""
# Numerical differentiation

Finite difference approximation of the derivative/gradient based on a number of function evaluations.

Often based on the limit definition of a derivative. Theoretical analysis using Taylor approximation:

$$f(x + h) = f(x) + \frac{h}{1!}f'(x) + \frac{h^2}{2!}f''(x) + \frac{h^3}{3!}f^{(3)}(x)+\ldots$$

**Forward difference**

$$f'(x)\approx \frac{f(x+h) - f(x)}{h}$$

**Central difference**

$$f'(x)\approx \frac{f(x+h) - f(x-h)}{2h}$$

**Complex step method**

$$f'(x)\approx \frac{\text{Im}(f(x +ih))}{h}$$
"""

# ╔═╡ 0238579a-fe70-11ea-0968-0f11ae9557db
md"**Assignment**: Implement these functions."

# ╔═╡ eef066e0-f2b3-11ea-3cfa-5f35ebf13f76
diff_fordiff(f, x; h=1e-10) = (f(x+h) - f(x)) / h;

# ╔═╡ f2158710-f2b3-11ea-16c9-dded5e1747d4
diff_centrdiff(f, x; h=1e-10) = (f(x+h)-f(x-h))/ 2h;

# ╔═╡ f54ee2b2-f2b3-11ea-24e8-f556ab9e3b85
diff_complstep(f, x; h=1e-10) = imag(f(x+im*h)) / h;

# ╔═╡ 6d18ecec-07a7-11eb-0a75-b1e47ef58f95
abs(df(a) - diff_fordiff(f, a; h=1e-20))

# ╔═╡ b23725a0-07a7-11eb-2813-054ec2b6ae5d
diff_fordiff(f, a)

# ╔═╡ 1055a996-f2b4-11ea-14d2-6116cbbaa6ba
diff_centrdiff(f, a)

# ╔═╡ 13db037a-f2b4-11ea-2a38-bd7c1c60f7d9
diff_complstep(f, a)

# ╔═╡ 4f42234e-f2b4-11ea-2c86-93449b8dc60b
md"""
## Intermezzo: floats

Real numbers are always represented as floating point numbers in a computer.

![Encoding of a real number using a `Float32`.](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/floats.png)

By default, Julia uses double precision floats (`Float64`). For brevity, let us take a look at the bit representation of a float. We use `Float32` for brevity's sake.
"""

# ╔═╡ 6adbc7fe-f2b4-11ea-1861-b398d517edfa
num = Float32(10.789)

# ╔═╡ 6cd1accc-f2b4-11ea-0077-e7685b77b1be
bitstring(num)

# ╔═╡ 77d4fd72-f2b4-11ea-0bd6-97fc2299ca2f
md"The first bit encodes the *sign*, here positive."

# ╔═╡ 773af248-f2b4-11ea-154a-cb67beebed3d
sign(num)

# ╔═╡ 8a5ee764-f2b4-11ea-23f3-27664d910c16
md"The next eight bits specify the *exponent*, the magnitude of the number."

# ╔═╡ 8ded9128-f2b4-11ea-0f44-9321077205de
exponent(num)

# ╔═╡ c3a1abd8-f2b4-11ea-0509-a5296ae8b361
md"While the final 23 bits specify the *mantissa*, a number between $[1,2]$ representing the precision."

# ╔═╡ c9f92e66-f2b4-11ea-322f-7778f149453d
significand(num)

# ╔═╡ cfb07cc6-f2b4-11ea-08a7-dd0b13d15b5e
md"These can be used to reconstrunct the number."

# ╔═╡ d56b74f2-f2b4-11ea-12b4-5f54d4999b38
significand(num) * 2^exponent(num)

# ╔═╡ d9f2d114-f2b4-11ea-32b0-43fa1a872a97
md"The *machine precision* of a number can be retained using `eps`. This is the relative error."

# ╔═╡ de1fd700-f2b4-11ea-36dd-4bfba5ecc134
eps(1.2)

# ╔═╡ df14ce72-f2b4-11ea-307b-8fac2889681b
eps(1.2e10)

# ╔═╡ e584f976-f2b4-11ea-0245-67f6bcc4be0e
eps(1.2e-10)

# ╔═╡ 4eff1680-1210-11eb-0afb-99f8fcc72aad
2+2

# ╔═╡ eb68674c-f2b4-11ea-0cc2-f792d17a5005
md"""
This brings us with numerical issues we might encounter using numerical differentiation.

**First sin of numerical analysis**:

> *thou shalt not add small numbers to big numbers*

**second sin of numerical analysis**:

> *thou shalt not subtract numbers which are approximately equal*
"""

# ╔═╡ f4edad2c-f2b4-11ea-05fb-3d43c18e9d5c
md"### Back to numerical differentiation"

# ╔═╡ fd9ed266-f2b4-11ea-3c53-3bd2fe30a7f1
begin
	fexamp(x) = 64x*(1-x)*(1-2x)^2*(1-8x+8x^2)^2
	dfexamp = diff(fexamp(x), x)
	error(diff, h; x=1.0) = max(abs(Float64(dfexamp(x)) - diff(fexamp, x, h=h)), 1e-50)
	stepsizes = map(t->10.0^t, -20:0.1:-1);
	plot(stepsizes, error.(diff_fordiff, stepsizes), label="forward difference",
    xscale=:log10, yscale=:log10, lw=2, legend=:bottomright, color=:blue)
	plot!(stepsizes, error.(diff_centrdiff, stepsizes), label="central difference", 		lw=2, color=:red)
	plot!(stepsizes, error.(diff_complstep, stepsizes), label="complex step", lw=2,
            color=:yellow)
	xlabel!("\$h\$")
	ylabel!("absolute error")
end

# ╔═╡ 4e6b2b04-f2b5-11ea-0c52-5f32438e8c68
md"""
Advantages of numerical differentiation:
- easy to implement
- general, no assumptions needed

Disadvantages:
- not numerically stable (round-off errors)
- not efficient for gradients ($\mathcal{O}(n)$ evaluations for $n$-dimensional vectors)
"""

# ╔═╡ 55ff7230-f2b5-11ea-2898-db5e21c764fd
md"""## Approximations of multiplications with gradients

**Gradient-vector approximation**

$$\nabla f(\mathbf{x})^\intercal \mathbf{d} \approx \frac{f(\mathbf{x}+h\cdot\mathbf{d}) - f(\mathbf{x}-h\cdot\mathbf{d})}{2h}$$

**Hessian-vector approximation**

$$\nabla^2 f(\mathbf{x}) \mathbf{d} \approx \frac{\nabla f(\mathbf{x}+h\cdot\mathbf{d}) - \nabla f(\mathbf{x}-h\cdot\mathbf{d})}{2h}$$
"""

# ╔═╡ 688a1f40-f2b5-11ea-3fa8-ebb735323ddd
grad_vect(f, x, d; h=1e-10) = (f(x + h * d) - f(x - h * d)) / (2h)

# ╔═╡ 6d6da8ba-f2b5-11ea-1064-df27ecea7d13
dvect = randn(10) / 10

# ╔═╡ 6fd08e42-f2b5-11ea-36e6-d7ff10412e2b
xvect = 2rand(10)

# ╔═╡ 9a87095e-f2b5-11ea-334e-57a637e20c43
A = randn(10, 10) |> A -> A * A' / 100

# ╔═╡ 9345b744-ff03-11ea-1c19-1913536154bb
md"$$g(\mathbf{x}) = \exp(-\mathbf{x}^\intercal A\mathbf{x})$$"

# ╔═╡ 75ae48fe-f2b5-11ea-0614-912091cbfb6a
g(x) = exp(- sum(x .* (A * x)))

# ╔═╡ 8a738272-f2b5-11ea-2de6-f9ae6b30a466
md"Correct gradient and Hessian (by hand)"

# ╔═╡ 7b95d73c-f2b5-11ea-1b29-858f7b6c99bf
Dg(x) = -2g(x) * A * x

# ╔═╡ 836ca94a-f2b5-11ea-07e0-bd830f810fba
D²g(x) = -2g(x) * A - 2A * x * Dg(x)'

# ╔═╡ 8fe85430-f2b5-11ea-3fe7-5d9a0db03932
g(xvect)

# ╔═╡ b372c25c-f2b5-11ea-1c55-e5045e0dc5b7
Dg(xvect)

# ╔═╡ b6e48b12-f2b5-11ea-3bf2-81e5b2a3dc98
Dg(xvect)' * dvect

# ╔═╡ c0960f82-f2b5-11ea-0df4-c15d19a698f4
h = 1e-10  # small stepsize

# ╔═╡ bebbd674-f2b5-11ea-19c6-4fe1257b0277
(Dg(xvect + h * dvect) - Dg(xvect - h * dvect)) / 2h

# ╔═╡ d6c6d4ea-f2b7-11ea-0b41-d7f629cdcd2f
md"""
## Forward differentiation

Accumulation of the gradients along the *computational graph*.

![](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/forwarddiff.png)

Forward differentiation computes the gradient from the inputs to the outputs.

### Differentiation rules

**Sum rule**:

$$\frac{\partial (f(x)+g(x))}{\partial x} =  \frac{\partial f(x)}{\partial x} + \frac{\partial f(x)}{\partial x}$$

**Product rule**:

$$\frac{\partial (f(x)g(x))}{\partial x} =  f(x)\frac{\partial g(x)}{\partial x} + g(x)\frac{\partial f(x)}{\partial x}$$

**Chain rule**:

$$\frac{\partial (g(f(x))}{\partial x} =  \frac{\partial g(u)}{\partial u}\mid_{u=f(x)} \frac{\partial f(x)}{\partial x}$$

## Dual numbers

Forward differentiation can be viewed as evaluating function using *dual numbers*, which can be viewed as truncated Taylor series:

$$v + \dot{v}\epsilon\,,$$

where $v,\dot{v}\in\mathbb{R}$ and $\epsilon$ a nilpotent number, i.e. $\epsilon^2=0$. For example, we have

$$(v + \dot{v}\epsilon) + (u + \dot{u}\epsilon) = (v+u) + (\dot{v} +\dot{u})\epsilon$$

$$(v + \dot{v}\epsilon)(u + \dot{u}\epsilon) = (vu) + (v\dot{u} +\dot{v}u)\epsilon\,.$$

These dual numbers can be used as

$$f(v+\dot{v}\epsilon) = f(v) + f'(v)\dot{v}\epsilon\,.$$
"""

# ╔═╡ fe18341e-f2b7-11ea-0236-03b44eebd328
struct Dual{T}
    v::T
    vdot::T
end

# ╔═╡ 0aba93a6-f2b8-11ea-1461-d96878b7a5c1
md"Let's implement some basic rules showing linearity."

# ╔═╡ 014ec92c-f2b8-11ea-3986-8d74a8c2458d
begin
	Base.:+(a::Dual, b::Dual) = Dual(a.v + b.v, a.vdot + b.vdot)
	Base.:*(a::Dual, b::Dual) = Dual(a.v * b.v, a.v * b.vdot + b.v * a.vdot)
	Base.:+(c::Real, b::Dual) = Dual(c + b.v, b.vdot)
	Base.:*(v::Real, b::Dual) = Dual(v, 0.0) * b
	Base.:*(a::Dual, b::Real) = v * a
end

# ╔═╡ 3f5a58da-f2b8-11ea-3ab0-5db5e9199bb1
md"And some more advanced ones, based on differentiation."

# ╔═╡ 48e5da96-f2b8-11ea-0108-2104f2c8ac24
Base.:/(a::Dual, b::Dual) = Dual(a.v / b.v, (a.vdot * b.v - a.v * b.vdot) / b.v^2)

# ╔═╡ 263e18dc-fe70-11ea-03a0-8fd5e031b06f
md"**Assignment**: complete this code."

# ╔═╡ 2f18c362-f2b8-11ea-1beb-719ff5cd5f8d
Base.:sin(a::Dual) = Dual(sin(a.v), cos(a.v)*a.vdot)

# ╔═╡ 393f92a6-fe70-11ea-024e-ff30dcfc7103
Base.:cos(a::Dual) = Dual(cos(a.v), -sin(a.v)*a.vdot)

# ╔═╡ 329cb67e-f2b8-11ea-07b1-e13fd8ddbe42
Base.:exp(a::Dual) = Dual(exp(a.v), exp(a.v)*a.vdot)

# ╔═╡ 36a2ca10-f2b8-11ea-0fae-53b552c76115
Base.:log(a::Dual) = Dual(log(a.v), 1/ a.v * a.vdot)

# ╔═╡ 4ee45daa-f2b8-11ea-10ea-f929853879b1
f(Dual(a, 1.0))

# ╔═╡ 555c8d88-f2b8-11ea-34d8-a17ce87ccdc1
md"This directly works for vectors!"

# ╔═╡ 59eee760-f2b8-11ea-29f6-732292e2f95d
q(x) = 10.0 * x[1] * x[2] + x[1] * x[1] + sin(x[1]) / x[2]

# ╔═╡ 5ccd625e-f2b8-11ea-38c1-098e2f31ce53
q([1, 2])

# ╔═╡ 60b52f6e-f2b8-11ea-2afc-4be2d4ae3e19
q(Dual.([1, 2], [1, 0]))  # partial wrt x1

# ╔═╡ 63d5e0be-f2b8-11ea-16e8-d9528b4e8945
q(Dual.([1, 2], [0, 1]))  # partial wrt x2

# ╔═╡ 6971e7e6-f2b8-11ea-2a76-9dda6b1455b6
md"In practice, we prefer to use a package to do this."

# ╔═╡ 723fe4ca-f2b8-11ea-1440-a9e1cd213b66
ForwardDiff.derivative(f, a)

# ╔═╡ 76284520-f2b8-11ea-077e-43167cc7484c
ForwardDiff.gradient(g, xvect)

# ╔═╡ 85a5c518-f2b8-11ea-2e12-25d16a12e7c8
ForwardDiff.gradient(q, [1, 2])

# ╔═╡ 8a254852-f2b8-11ea-142d-71bb270a5c01
md"""
Forward differentiation:

- exact gradients!
- computational complexity scales with **number of inputs**
- used when you have more outputs than inputs
"""

# ╔═╡ 91d0e03e-f2b8-11ea-07e5-db056636ac21
md"""
# Reverse differentiation

Compute the gradient from the output toward the inputs using the chain rule.

![](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/reversediff.png)

Reverse differentiation:

- also exact!
- main workhorse for training artificial neural networks.
- efficient when more inputs than outputs (machine learning: thousands of parameters vs. one loss)
"""

# ╔═╡ 3c39d114-1214-11eb-07a5-81df5b18bd3f
f(x)

# ╔═╡ a041fbf8-f2b8-11ea-28ba-0772497c649f
md"Constructing the derivative or gradient can be done by appending `'` to the function."

# ╔═╡ 9b5db398-f2b8-11ea-00bc-77de99792925
f'(a)  # that's it

# ╔═╡ bdce5964-f2b8-11ea-03d1-57705d295040
md"More verbose, but exactly the same:"

# ╔═╡ c6728784-f2b8-11ea-12b5-d7eb039d77ce
Zygote.gradient(f, a)  # returns a tuple, since you can differentiate wrt multiple arguments

# ╔═╡ e1057afe-fe6e-11ea-27fe-6def4d8e1f2d
md"This works for function with multiple inputs."

# ╔═╡ eb76e252-fe6e-11ea-0939-c51ff3c73827
Zygote.gradient((x1, x2) -> (x1^2 + x2^2 - 0.1x1*x2)/(x1+1.0), 0.2, 0.3)

# ╔═╡ d6d310ee-f2b8-11ea-24e1-3b405f071545
md"Fuctions with vectorial input."

# ╔═╡ d3ae2c14-f2b8-11ea-1ac2-fb5e0b4b8500
g'(xvect)

# ╔═╡ dbb1e68c-f2b8-11ea-1220-6f5715822e3d
md"Finding the Hessian:"

# ╔═╡ df557ed2-f2b8-11ea-22ee-e1c90a2655bf
Zygote.hessian(g, xvect)

# ╔═╡ ee52d2b8-f2b8-11ea-15c4-19e9f2ad10e2
md"""
## Artificial neural networks

Multi-layer perceptron.

![](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/ANN_example.png)

Forward differentiation.

![](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/Forwardprop.png)

Reverse differentation or backpropagation.

![](https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/03.AutoDiff/Figures/Backprop.png)

Returns effect of changing layer output on the loss. Can be related directly to the parameters!
"""

# ╔═╡ 213fc38a-f380-11ea-135e-6f72924c19aa
md"""
## Differentiating complex objects

Automatic differentiation can be used beyond machine learning and optimization:

- [physical engines](https://arxiv.org/abs/1611.01652) to learn robot control
- differentiating [protein](https://github.com/lupoglaz/TorchProteinLibrary) [structures](https://www.cell.com/cell-systems/fulltext/S2405-4712(19)30076-6)
- Sinkhorn algorithm
- [dynamic programming](https://arxiv.org/abs/1802.03676)
- [differential equations](https://julialang.org/blog/2019/01/fluxdiffeq)

Everything is computed by some straightforward and differentiable functions!
"""

# ╔═╡ 79221eee-f2be-11ea-3e43-4596f5ce5694
md"""
## Application: inverse kinematics of a robot arm

As an illustration of using gradients, let us study a [simple robot arm](https://appliedgo.net/roboticarm/). This arm consists of three joints that can be moved in the 2D plane. These lengths of the three segments are, respectivly 2, 1 and 1.5 meter.

The forward kinetics (i.e., going from the angles of the joints to the position of the arms) can be computed by the following functions.
"""

# ╔═╡ fb91a9c4-f2c0-11ea-38c4-5b554446e997
pos_second_joint((θ₁, θ₂, θ₃)) = [2cos(θ₁), 2sin(θ₁)]

# ╔═╡ 07c45d78-f2c3-11ea-131e-b9f0b67ca3cd
pos_third_joint((θ₁, θ₂, θ₃)) = pos_second_joint((θ₁, θ₂, θ₃)) + [cos( θ₂), sin( θ₂)]

# ╔═╡ 33023d10-f2c1-11ea-3375-e9ac23172ab6
pos_hand((θ₁, θ₂, θ₃)) = pos_third_joint((θ₁, θ₂, θ₃)) .+ [1.5cos(θ₃), 1.5sin(θ₃)]

# ╔═╡ 31cdefd2-f381-11ea-3037-e9c023ebc1e9
md"Suppose that there is a teapot at a certain location we want to grasp with the arm."

# ╔═╡ 4def0626-f2c1-11ea-2f01-31a8c42b121c
position_teapot = [-1, 2.3]

# ╔═╡ 46c54108-f381-11ea-141f-31a696a7a8ef
md"Below is a plot of the arm and the teapot. Can you turn the angles to bring the hand to the position of the teapot?"

# ╔═╡ 8049560e-f383-11ea-1c8e-39c629d20982
@bind θ₁ Slider(0:0.01:2π, show_value=true)

# ╔═╡ c7b5a1aa-f383-11ea-2a84-8da9c9e60dbc
@bind θ₂ Slider(0:0.01:2π, show_value=true)

# ╔═╡ cfad1a8c-f383-11ea-3dad-214f3da9485d
@bind θ₃ Slider(0:0.01:2π, show_value=true)

# ╔═╡ 83ddbe6a-f381-11ea-2415-255cbfb2590b
md"First complete a function to determine the distance between the arm and the teapot."

# ╔═╡ 30e9891c-f2c5-11ea-022d-5dbf9c064e43
distance_to_teapot(θ) = (position_teapot .- pos_hand(θ)).^2 |> sum |> sqrt

# ╔═╡ 5812f84c-f2c1-11ea-0f6c-31b9fedc3704
begin
	psj = pos_second_joint((θ₁, θ₂, θ₃))
	ptj = pos_third_joint((θ₁, θ₂, θ₃))
	phand = pos_hand((θ₁, θ₂, θ₃))
	plot([0, psj[1], ptj[1], phand[1]], [0, psj[2],ptj[2], phand[2]], label="Arm")
	scatter!([0], [0], label="First joint (fixed)")
	scatter!([psj[1]], [psj[2]], label="Second joint")
	scatter!([ptj[1]], [ptj[2]], label="Second joint")
	scatter!([phand[1]], [phand[2]], label="Hand")
	scatter!([position_teapot[1]], [position_teapot[2]], m=:star, label="teapot")
	title!("Distance to teapot = $(distance_to_teapot((θ₁, θ₂, θ₃)))")
	xlims!(-5, 5)
	ylims!(-5, 5)
end

# ╔═╡ b07cb958-f381-11ea-3378-5517440464ca
md"Can you use gradients to position the arm?"

# ╔═╡ eaac2514-f8fe-11ea-2229-81f5b59f58fe
distance_to_teapot([θ₁, θ₂, θ₃])

# ╔═╡ fb044216-f381-11ea-2faf-fd53dad74950
distance_to_teapot'([θ₁, θ₂, θ₃])

# ╔═╡ 675b222c-f382-11ea-257a-afcb5a5afdd8
  # use simple gradient descent to move the arm to the teapot
begin
	#initial value of 0
	θ = [0.0, 0.0, 0.0]
	for i in 1:100
		θ .-=0.05distance_to_teapot'(θ)
	end
end

# ╔═╡ ccd482ae-1215-11eb-3c4b-a342a810ec9f
θ

# ╔═╡ 881dedd4-f385-11ea-342e-a3700e90ac66
md"""
## Exercise

Consider the *Wheeler's Ridge* function:

$$f(\mathbf{x}) = -\exp(-(x_1 x_2 - a)^2 -(x_2 -a)^2)\,,$$

at the point $\mathbf{x}_0=[1.5, 1]^T$. We set $a=1.5$.

Implement this function.
"""

# ╔═╡ c1e19e26-f385-11ea-17b6-6d2d2e5cb7c0
fwr(x) = missing

# ╔═╡ fe18a0a0-f8fe-11ea-1595-4d95256dabca
fwr([1,2]) isa Missing || contour(-2:0.1:2, -2:0.1:2, (x1, x2) -> fwr([x1, x2]))

# ╔═╡ df5f20d8-f385-11ea-104d-9b9e7b094eb7
md"""

**Assignments**

1. Compute the gradient by hand.
2. Find the gradient and Hessian at $\mathbf{x}_0$ by numerical differentiation.
3. Compute the gradient and Hessian at $\mathbf{x}_0$ using automatic differentiation.
4. (optional) Use the function `quiver!` to draw the gradient as a vector field on the contour plot.

"""

# ╔═╡ e4a088f0-f385-11ea-00da-b12e40e4357f
x₀ = [1.5, 1.0]  # btw written as x\_0<TAB>

# ╔═╡ 8f3fecac-f8fd-11ea-136d-1513f60ee3ed

# ╔═╡ f8863a18-f385-11ea-0965-0f074e58d9b2

# ╔═╡ 4c751e90-f387-11ea-2358-273423b85c00

# ╔═╡ 4d2c34ae-f387-11ea-1649-19c3a4520f74

# ╔═╡ 34d48ec0-f2b9-11ea-147e-d3554a9fcc04
md"""
# References

- Gunes et. al. (2015) *Automatic differentiation in machine learning: a survey*
- Kochenderfer, M. J. and Wheeler, T., '*Algorithms for Optimization*'. MIT Press (2019)
"""

# ╔═╡ Cell order:
# ╟─08cf2a80-f2b2-11ea-132d-2f26c934bbb2
# ╠═47f98a52-f2b2-11ea-0ac6-9b4b74b7b15d
# ╟─5154a672-f2b2-11ea-1280-3ba5dadff7f9
# ╟─f53fbf24-f2b2-11ea-2f59-d943884fee9d
# ╟─e6a55456-fe6f-11ea-328c-47fedb388194
# ╠═2be8656c-f2b3-11ea-2244-5740c807deff
# ╟─2e00ef04-f2b3-11ea-1960-8136321c71db
# ╟─f7cd052e-f2bb-11ea-239a-47966e96e909
# ╠═7244eba2-f2b3-11ea-045f-ffd208852eda
# ╟─401df7a4-f2b3-11ea-12a1-e9b46bcfcf1b
# ╠═5ac4989c-f2b3-11ea-3e8a-a79a4117e48c
# ╠═5db7350a-f2b3-11ea-293d-ffc572e0dbe8
# ╠═bbe02ec4-120d-11eb-39fb-8545d56de2ad
# ╠═64b28e4a-f2b3-11ea-0fc0-7d37b4794686
# ╠═684ff7f4-f2b3-11ea-14f7-7d3f329ecd79
# ╟─7a6a18de-f2b3-11ea-1b3d-ef7a2b1b0434
# ╟─913234a2-f2b3-11ea-3501-e3bfc54f91fe
# ╟─0238579a-fe70-11ea-0968-0f11ae9557db
# ╠═eef066e0-f2b3-11ea-3cfa-5f35ebf13f76
# ╠═f2158710-f2b3-11ea-16c9-dded5e1747d4
# ╠═f54ee2b2-f2b3-11ea-24e8-f556ab9e3b85
# ╠═6d18ecec-07a7-11eb-0a75-b1e47ef58f95
# ╠═b23725a0-07a7-11eb-2813-054ec2b6ae5d
# ╠═1055a996-f2b4-11ea-14d2-6116cbbaa6ba
# ╠═13db037a-f2b4-11ea-2a38-bd7c1c60f7d9
# ╟─4f42234e-f2b4-11ea-2c86-93449b8dc60b
# ╠═6adbc7fe-f2b4-11ea-1861-b398d517edfa
# ╠═6cd1accc-f2b4-11ea-0077-e7685b77b1be
# ╟─77d4fd72-f2b4-11ea-0bd6-97fc2299ca2f
# ╠═773af248-f2b4-11ea-154a-cb67beebed3d
# ╟─8a5ee764-f2b4-11ea-23f3-27664d910c16
# ╠═8ded9128-f2b4-11ea-0f44-9321077205de
# ╟─c3a1abd8-f2b4-11ea-0509-a5296ae8b361
# ╠═c9f92e66-f2b4-11ea-322f-7778f149453d
# ╟─cfb07cc6-f2b4-11ea-08a7-dd0b13d15b5e
# ╠═d56b74f2-f2b4-11ea-12b4-5f54d4999b38
# ╟─d9f2d114-f2b4-11ea-32b0-43fa1a872a97
# ╠═de1fd700-f2b4-11ea-36dd-4bfba5ecc134
# ╠═df14ce72-f2b4-11ea-307b-8fac2889681b
# ╠═e584f976-f2b4-11ea-0245-67f6bcc4be0e
# ╠═4eff1680-1210-11eb-0afb-99f8fcc72aad
# ╟─eb68674c-f2b4-11ea-0cc2-f792d17a5005
# ╟─f4edad2c-f2b4-11ea-05fb-3d43c18e9d5c
# ╟─fd9ed266-f2b4-11ea-3c53-3bd2fe30a7f1
# ╟─4e6b2b04-f2b5-11ea-0c52-5f32438e8c68
# ╟─55ff7230-f2b5-11ea-2898-db5e21c764fd
# ╠═688a1f40-f2b5-11ea-3fa8-ebb735323ddd
# ╠═6d6da8ba-f2b5-11ea-1064-df27ecea7d13
# ╠═6fd08e42-f2b5-11ea-36e6-d7ff10412e2b
# ╠═9a87095e-f2b5-11ea-334e-57a637e20c43
# ╟─9345b744-ff03-11ea-1c19-1913536154bb
# ╠═75ae48fe-f2b5-11ea-0614-912091cbfb6a
# ╟─8a738272-f2b5-11ea-2de6-f9ae6b30a466
# ╠═7b95d73c-f2b5-11ea-1b29-858f7b6c99bf
# ╠═836ca94a-f2b5-11ea-07e0-bd830f810fba
# ╠═8fe85430-f2b5-11ea-3fe7-5d9a0db03932
# ╠═b372c25c-f2b5-11ea-1c55-e5045e0dc5b7
# ╠═b6e48b12-f2b5-11ea-3bf2-81e5b2a3dc98
# ╠═c0960f82-f2b5-11ea-0df4-c15d19a698f4
# ╠═bebbd674-f2b5-11ea-19c6-4fe1257b0277
# ╟─d6c6d4ea-f2b7-11ea-0b41-d7f629cdcd2f
# ╠═fe18341e-f2b7-11ea-0236-03b44eebd328
# ╟─0aba93a6-f2b8-11ea-1461-d96878b7a5c1
# ╠═014ec92c-f2b8-11ea-3986-8d74a8c2458d
# ╟─3f5a58da-f2b8-11ea-3ab0-5db5e9199bb1
# ╠═48e5da96-f2b8-11ea-0108-2104f2c8ac24
# ╟─263e18dc-fe70-11ea-03a0-8fd5e031b06f
# ╠═2f18c362-f2b8-11ea-1beb-719ff5cd5f8d
# ╠═393f92a6-fe70-11ea-024e-ff30dcfc7103
# ╠═329cb67e-f2b8-11ea-07b1-e13fd8ddbe42
# ╠═36a2ca10-f2b8-11ea-0fae-53b552c76115
# ╠═4ee45daa-f2b8-11ea-10ea-f929853879b1
# ╟─555c8d88-f2b8-11ea-34d8-a17ce87ccdc1
# ╠═59eee760-f2b8-11ea-29f6-732292e2f95d
# ╠═5ccd625e-f2b8-11ea-38c1-098e2f31ce53
# ╠═60b52f6e-f2b8-11ea-2afc-4be2d4ae3e19
# ╠═63d5e0be-f2b8-11ea-16e8-d9528b4e8945
# ╟─6971e7e6-f2b8-11ea-2a76-9dda6b1455b6
# ╠═6fe82e5a-f2b8-11ea-101d-d1ee5f9bd5b1
# ╠═723fe4ca-f2b8-11ea-1440-a9e1cd213b66
# ╠═76284520-f2b8-11ea-077e-43167cc7484c
# ╠═85a5c518-f2b8-11ea-2e12-25d16a12e7c8
# ╟─8a254852-f2b8-11ea-142d-71bb270a5c01
# ╟─91d0e03e-f2b8-11ea-07e5-db056636ac21
# ╠═963c4910-f2b8-11ea-12bb-bb42539b1c93
# ╠═3c39d114-1214-11eb-07a5-81df5b18bd3f
# ╟─a041fbf8-f2b8-11ea-28ba-0772497c649f
# ╠═9b5db398-f2b8-11ea-00bc-77de99792925
# ╟─bdce5964-f2b8-11ea-03d1-57705d295040
# ╠═c6728784-f2b8-11ea-12b5-d7eb039d77ce
# ╟─e1057afe-fe6e-11ea-27fe-6def4d8e1f2d
# ╠═eb76e252-fe6e-11ea-0939-c51ff3c73827
# ╟─d6d310ee-f2b8-11ea-24e1-3b405f071545
# ╠═d3ae2c14-f2b8-11ea-1ac2-fb5e0b4b8500
# ╟─dbb1e68c-f2b8-11ea-1220-6f5715822e3d
# ╠═df557ed2-f2b8-11ea-22ee-e1c90a2655bf
# ╟─ee52d2b8-f2b8-11ea-15c4-19e9f2ad10e2
# ╟─213fc38a-f380-11ea-135e-6f72924c19aa
# ╟─79221eee-f2be-11ea-3e43-4596f5ce5694
# ╠═fb91a9c4-f2c0-11ea-38c4-5b554446e997
# ╠═07c45d78-f2c3-11ea-131e-b9f0b67ca3cd
# ╠═33023d10-f2c1-11ea-3375-e9ac23172ab6
# ╟─31cdefd2-f381-11ea-3037-e9c023ebc1e9
# ╠═4def0626-f2c1-11ea-2f01-31a8c42b121c
# ╟─46c54108-f381-11ea-141f-31a696a7a8ef
# ╟─5812f84c-f2c1-11ea-0f6c-31b9fedc3704
# ╠═8049560e-f383-11ea-1c8e-39c629d20982
# ╠═c7b5a1aa-f383-11ea-2a84-8da9c9e60dbc
# ╠═cfad1a8c-f383-11ea-3dad-214f3da9485d
# ╟─83ddbe6a-f381-11ea-2415-255cbfb2590b
# ╠═30e9891c-f2c5-11ea-022d-5dbf9c064e43
# ╟─b07cb958-f381-11ea-3378-5517440464ca
# ╠═eaac2514-f8fe-11ea-2229-81f5b59f58fe
# ╠═fb044216-f381-11ea-2faf-fd53dad74950
# ╠═675b222c-f382-11ea-257a-afcb5a5afdd8
# ╠═ccd482ae-1215-11eb-3c4b-a342a810ec9f
# ╟─881dedd4-f385-11ea-342e-a3700e90ac66
# ╠═c1e19e26-f385-11ea-17b6-6d2d2e5cb7c0
# ╟─fe18a0a0-f8fe-11ea-1595-4d95256dabca
# ╟─df5f20d8-f385-11ea-104d-9b9e7b094eb7
# ╠═e4a088f0-f385-11ea-00da-b12e40e4357f
# ╠═8f3fecac-f8fd-11ea-136d-1513f60ee3ed
# ╠═f8863a18-f385-11ea-0965-0f074e58d9b2
# ╠═4c751e90-f387-11ea-2358-273423b85c00
# ╠═4d2c34ae-f387-11ea-1649-19c3a4520f74
# ╟─34d48ec0-f2b9-11ea-147e-d3554a9fcc04

max_stars_repo_name: tjvneste/STMO

max_issues_repo_name: tjvneste/STMO

max_forks_repo_name: tjvneste/STMO


text:
module Intcode

Base.Experimental.@optlevel 3

export IntcodeMachine
export run_intcode!, reset!
export set!

import Base: show, fetch, length

const AbstractTape = AbstractVector{<:Integer}

abstract type AbstractMode end
struct PositionMode <: AbstractMode end
struct ImmediateMode <: AbstractMode end

const _m_position = PositionMode()
const _m_immediate = ImmediateMode()

show(io::IO, ::ImmediateMode) = print(io, "MI")
show(io::IO, ::PositionMode) = print(io, "MP")

abstract type AbstractOp end
struct Add <: AbstractOp end; length(::Add) = 4
struct Mul <: AbstractOp end; length(::Mul) = 4
struct Input <: AbstractOp end; length(::Input) = 2
struct Output <: AbstractOp end; length(::Output) = 2
struct JumpTrue <: AbstractOp end; length(::JumpTrue) = 3
struct JumpFalse <: AbstractOp end; length(::JumpFalse) = 3
struct Less <: AbstractOp end; length(::Less) = 4
struct Equals <: AbstractOp end; length(::Equals) = 4
struct Halt <: AbstractOp end; length(::Halt) = 1

struct JumpValue
    value::Int
end

show(io::IO, jp::JumpValue) = print(io, "Jump($(jp.value))")

const _op_add = Add()
const _op_mul = Mul()
const _op_in = Input()
const _op_out = Output()
const _op_jpt = JumpTrue()
const _op_jpf = JumpFalse()
const _op_ls = Less()
const _op_eq = Equals()
const _op_halt = Halt()

const _opcode_map = Dict(
    [1,0] => _op_add,
    [2,0] => _op_mul,
    [3,0] => _op_in,
    [4,0] => _op_out,
    [5,0] => _op_jpt,
    [6,0] => _op_jpf,
    [7,0] => _op_ls,
    [8,0] => _op_eq,
    [9,9] => _op_halt,
)

show(io::IO, ::Add) = print(io, "A")
show(io::IO, ::Mul) = print(io, "M")
show(io::IO, ::Input) = print(io, "I")
show(io::IO, ::Output) = print(io, "O")
show(io::IO, ::JumpTrue) = print(io, "JT")
show(io::IO, ::JumpFalse) = print(io, "JF")
show(io::IO, ::Less) = print(io, "LS")
show(io::IO, ::Equals) = print(io, "EQ")
show(io::IO, ::Halt) = print(io, "H")

struct IntcodeInstruction{
    M1 <: AbstractMode,
    M2 <: AbstractMode,
    M3 <: AbstractMode,
    Opcode <: AbstractOp,
    Tape <: AbstractTape,
}
    mode::Tuple{M1,M2,M3}
    opcode::Opcode
    data::Tape
end

show(io::IO, i::IntcodeInstruction) = print(io, "<$(i.opcode)$(i.mode): $(i.data)>")
length(i::IntcodeInstruction) = length(i.opcode)

abstract type AbstractState end
struct Running <: AbstractState end
struct Idle <: AbstractState end
struct Halted <: AbstractState end

show(io::IO, ::Running) = print(io, "RUNNING")
show(io::IO, ::Idle) = print(io, "IDLE")
show(io::IO, ::Halted) = print(io, "HALTED")

const _s_run = Running()
const _s_idle = Idle()
const _s_halt = Halted()

mutable struct IntcodeMachine{Tape <: AbstractTape}
    state::Union{Nothing,<:AbstractState}
    pc::Int
    tape::Tape
    init::Tape
    IntcodeMachine(tape::Tape) where {Tape <: AbstractTape} =
        new{Tape}(_s_idle, 0, copy(tape), copy(tape))
end

show(io::IO, ic::IntcodeMachine) = print(io, "Intcode ($(ic.state)): pc=$(ic.pc)")

@inline function opcode(ds::AbstractTape)::Union{Nothing,<:AbstractOp}
    1 <= length(ds) <= 2 || return nothing
    op = [0,0]
    r = 1:length(ds)
    op[r] .= ds
    res = get(_opcode_map, op, nothing)
    if isnothing(res)
        @debug "Opcode '$ds' unknown"
        return nothing
    end
    return res
end

@inline function modes(d::Int)::Union{Nothing,<:AbstractMode}
    d == 0 && return _m_position
    d == 1 && return _m_immediate
    @debug "Mode '$d' unknown"
    return nothing
end

@noinline function modes(ds::AbstractTape)::Union{Nothing,Tuple{<:AbstractMode,<:AbstractMode,<:AbstractMode}}
    cs = [0,0,0]
    cs[1:length(ds)] .= ds
    ms = modes.(cs)
    any(isnothing, ms) && return nothing
    tuple(ms...)
end

@inline function _opmode(ds::AbstractTape)
    op = opcode(view(ds, 1:min(2, length(ds))))
    isnothing(op) && return nothing
    ms = modes(view(ds, 3:length(ds)))
    isnothing(ms) && return nothing
    (ms, op)
end

function opmode(tape::AbstractTape)
    _opmode(digits(tape[1]))
end

@inline function data_seg(op::AbstractOp, tape::AbstractTape)::Union{Nothing,<:AbstractTape}
    length(op) == 1 && return Int[]
    length(op) == 2 && return view(tape, 2:2)
    length(op) == 3 && return view(tape, 2:3)
    length(op) == 4 && return view(tape, 2:4)
    nothing
end

@inline function decode(tape::AbstractTape)::Union{Nothing,IntcodeInstruction}
    decodes = opmode(tape)
    isnothing(decodes) && return nothing
    ms, op = decodes
    data = data_seg(op, tape)
    isnothing(data) && return nothing
    IntcodeInstruction(ms, op, data)
end

function (i::IntcodeInstruction)(tape::AbstractTape)::Union{Nothing,Some{<:Any}}
    i.opcode(tape, i.mode, i.data)
end

function fetch_data(tape::AbstractTape, ::PositionMode, addr::Integer)
    tape[addr + 1]
end

function fetch_data(::AbstractTape, ::ImmediateMode, data::Integer)
    data
end

function set_data!(tape::AbstractTape, addr::Integer, v::Integer)
    tape[addr + 1] = v
end

@inline function (op::Add)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    x = fetch_data(tape, mode[1], data[1])
    y = fetch_data(tape, mode[2], data[2])
    res = x + y
    set_data!(tape, data[3], res) |> Some
end

@inline function (op::Mul)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    x = fetch_data(tape, mode[1], data[1])
    y = fetch_data(tape, mode[2], data[2])
    res = x * y
    set_data!(tape, data[3], res) |> Some
end

function (op::Input)(
    tape::AbstractTape,
    ::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    input = nothing
    try
        input = parse(Int, readline())
    catch
        @error "Cannot parse input as Int: $input"
        return nothing
    end
    isnothing(input) && return nothing
    set_data!(tape, data[1], input) |> Some
end

function (op::Output)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    val = fetch_data(tape, mode[1], data[1])
    println(val)
    Some(val)
end

function (op::JumpTrue)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    fetch_data(tape, mode[1], data[1]) != 0 || return Some(nothing)
    fetch_data(tape, mode[2], data[2]) + 1 |> JumpValue |> Some
end

function (op::JumpFalse)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    fetch_data(tape, mode[1], data[1]) == 0 || return Some(nothing)
    fetch_data(tape, mode[2], data[2]) + 1 |> JumpValue |> Some
end

function (op::Less)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    x = fetch_data(tape, mode[1], data[1])
    y = fetch_data(tape, mode[2], data[2])
    set_data!(tape, data[3], x < y) |> Some
end

function (op::Equals)(
    tape::AbstractTape,
    mode::Tuple{M1,M2,M3},
    data::AbstractTape,
) where {M1 <: AbstractMode,M2 <: AbstractMode,M3 <: AbstractMode,}
    x = fetch_data(tape, mode[1], data[1])
    y = fetch_data(tape, mode[2], data[2])
    set_data!(tape, data[3], x == y) |> Some
end

function fetch_intcode!(ic::IntcodeMachine)::Union{Nothing,IntcodeInstruction}
    i = decode(view(ic.tape, ic.pc:length(ic.tape)))
    if isnothing(i)
        ic.state = nothing
        return nothing
    else
        return i
    end
end

should_halt(::Halt) = true
should_halt(::AbstractOp) = false
should_halt(i::IntcodeInstruction) = should_halt(i.opcode)
should_jump(res::Some{<:Any}) = something(res) isa JumpValue
should_jump(::Nothing) = false

function halt!(ic::IntcodeMachine)
    ic.state = _s_halt
    return Some(_s_halt)
end

function jump!(ic::IntcodeMachine, jp::JumpValue)
    ic.pc = jp.value
end

function exec!(ic::IntcodeMachine, i::IntcodeInstruction)::Union{Nothing,Some{<:Any}}
    should_halt(i) && return halt!(ic)
    res = i(ic.tape)
    if isnothing(res)
        ic.state = nothing
    elseif should_jump(res)
        jump!(ic, something(res))
    end
    return res
end

function increment!(ic::IntcodeMachine, i::IntcodeInstruction)
    if ic.pc > length(ic.tape)
        if ic.debug; debug_failure(ic, "pc past code") end
        return false
    end
    ic.pc += length(i)
    return true
end

isrunning(ic::IntcodeMachine) = ic.state === _s_run
issuccess(ic::IntcodeMachine) = !isnothing(ic.state)
isidle(ic::IntcodeMachine) = ic.state === _s_idle

function reset!(ic::IntcodeMachine)
    ic.state = _s_idle
    ic.pc = 0
    ic.tape = copy(ic.init)
end

function start!(ic::IntcodeMachine)
    if !isidle(ic)
        debug_failure(ic, "Intcode is not ready, please reset")
        return false
    end
    ic.state = _s_run
    ic.pc = 1
    return true
end

function run_intcode!(ic::IntcodeMachine)::Union{Nothing,Some{<:AbstractState}}
    start!(ic) || return nothing
    while isrunning(ic)
        i = fetch_intcode!(ic)
        issuccess(ic) || return debug_failure(ic, i)
        res = exec!(ic, i)
        issuccess(ic) || return debug_failure(ic, i, res)
        debug_info(ic, i, res)
        should_jump(res) || increment!(ic, i) || return debug_failure(ic, i, res)
    end
    issuccess(ic) || return nothing
    Some(ic.state)
end

function set!(ic::IntcodeMachine, noun::Integer, verb::Integer)
    ic.tape[2] = noun
    ic.tape[3] = verb
end

function fetch(ic::IntcodeMachine)::Union{Nothing,Some{<:Integer}}
    issuccess(ic) || return nothing
    ic.tape |> first |> Some
end

@inline function debug_info(
    ic::IntcodeMachine,
    i::Union{Nothing,AbstractOp,IntcodeInstruction},
    res
)
    @debug "$ic; i=$i; res=$res"
end

@inline function debug_failure(
    ic::IntcodeMachine,
    i::Union{Nothing,AbstractOp,IntcodeInstruction},
    res
)
    @debug "[FAILED] $ic; i=$i; res=$res"
end

@inline function debug_failure(
    ic::IntcodeMachine,
    i::Union{Nothing,AbstractOp,IntcodeInstruction}
)
    @debug "[FAILED] $ic; i=$i"
end

@inline function debug_failure(ic::IntcodeMachine, msg::AbstractString)
    @error "[FAILED] $msg" ic
end

end # module


text:
"""
Simplex Crossover (SPX).

`ϵ>0` controls how the original simplex is inflated, `ϵ=1` means no
inflation.

See Tsutsui, Yamamura & Higuchi "Multi-parent recombination with simplex crossover in real coded genetic algorithms", 1999,
Proc. of the Genetic and Evolutionary Computation Conference
"""
struct SimplexCrossover{NP} <: CrossoverOperator{NP,1}
    ϵ::Float64      # inflation rate

    SimplexCrossover{NP}(ϵ::Number = Math.sqrt(NP + 1)) where NP = new{NP}(ϵ)
    SimplexCrossover{NP}(params::Parameters) where NP = new{NP}(params[:SPX_ϵ])
end

const SPX_DefaultOptions = ParamsDict(
    :SPX_ϵ => 1.1 # how much to inflate
)

#masscenter(pop, parentIndices) = mapslices(mean, pop[:, parentIndices], 2)

function apply!(xover::SimplexCrossover,
                target::AbstractIndividual, targetIndex::Int,
                pop, parentIndices)
    @assert length(parentIndices) == numparents(xover)

    offset = fill!(similar(target), 0)
    parentSum = copy(viewer(pop, parentIndices[1]))
    for i in 1:(numparents(xover)-1)
        s = rand()^(1/i)
        parent = viewer(pop, parentIndices[i])
        parentNext = viewer(pop, parentIndices[i+1])
        @inbounds offset .= s * (xover.ϵ * (parent .- parentNext) .+ offset)
        @inbounds parentSum .+= parentNext
    end
    s = (1-xover.ϵ)/numparents(xover)
    parentLast = viewer(pop, parentIndices[end])
    @inbounds target .= s .* parentSum + xover.ϵ .* parentLast .+ offset
    return target
end


text:
using DifferentialEquations
using SimpleDiffEq
using Distributions
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_sde!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    ifrac = β*c*I/N*S*δt
    rfrac = γ*I*δt
    ifrac_noise = sqrt(ifrac)*rand(Normal(0,1))
    rfrac_noise = sqrt(rfrac)*rand(Normal(0,1))
    @inbounds begin
        du[1] = S-(ifrac+ifrac_noise)
        du[2] = I+(ifrac+ifrac_noise) - (rfrac + rfrac_noise)
        du[3] = R+(rfrac+rfrac_noise)
    end
    for i in 1:3
        if du[i] < 0 du[i]=0 end
    end
    nothing
end;

δt = 0.1
nsteps = 400
tmax = nsteps*δt
tspan = (0.0,nsteps)
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0]; # S,I,R

p = [0.05,10.0,0.25,δt]; # β,c,γ,δt

Random.seed!(1234);

prob_sde = DiscreteProblem(sir_sde!,u0,tspan,p)

sol_sde = solve(prob_sde,solver=FunctionMap);

df_sde = DataFrame(sol_sde')
df_sde[!,:t] = t;

@df df_sde plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_sde,solver=FunctionMap)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


text:
function genball(R)
    sphere_side = 2R + 1
    sphere = zeros(Bool, (sphere_side, sphere_side, sphere_side))

    for i in -R:R
        for j in -R:R
            for k in -R:R
                dist = i^2 + j^2 + k^2
                if dist < R^2
                    sphere[k+R+1, j+R+1, i+R+1] = 1
                end
            end
        end
    end

    return sphere
end

function gendisk(R)
    sphere  = zeros(Bool, (2R + 1, 2R + 1))

    for i in -R:R
        for j in -R:R
            dist = i^2 + j^2
            if dist < R^2
                sphere[j+R+1, i+R+1] = 1
            end
        end
    end

    return sphere
end

function gentorus(side, r, R)
    a = zeros(Bool, (side, side, side))
    for i in 1:side
        for j in 1:side
            for k in 1:side
                x = i - side÷2
                y = j - side÷2
                z = k - side÷2
                if (x^2 + y^2 + z^2 + R^2 - r^2)^2 < 4R^2*(x^2 + y^2)
                    a[k,j,i] = 1
                end
            end
        end
    end

    return a
end

@test euler_characteristic(genball(100)) == 1
@test euler_characteristic(gendisk(100)) == 1
@test euler_characteristic(ones(Bool, (100, 100, 100))) == 1
@test euler_characteristic(hcat(ones(Bool, (100, 100, 100)),
                                zeros(Bool, (100, 100, 100)),
                                ones(Bool, (100, 100, 100)))) == 2
let a = ones(Bool, (100, 100));
    a[50:60, 50:60] .= 0;
    @test euler_characteristic(a) == 0
end

let a = gentorus(200, 10, 90);
    b = cat(a[:,:,begin:end-5], a[:,:,begin+5:end]; dims = [3])
    @test euler_characteristic(b) == -1
end


text:
@testset "tabulate" begin
    @testset "basic" begin
        @system STabulate(Controller) begin
            T => [
                # a b
                  0 4 ; # A
                  1 5 ; # B
                  2 6 ; # C
                  3 7 ; # D
            ] ~ tabulate(rows=(:A, :B, :C, :D), columns=(:a, :b))
            Aa(T) => T.A.a ~ preserve
            Ba(T) => T.B[:a] ~ preserve
            Ca(T) => T[:C].a ~ preserve
            Da(T) => T[:D][:a] ~ preserve
            Ab(T.A.b) ~ preserve
            Bb(T.B[:b]) ~ preserve
            Cb(T[:C].b) ~ preserve
            Db(T[:D][:b]) ~ preserve
        end
        s = instance(STabulate)
        @test s.Aa' == 0
        @test s.Ba' == 1
        @test s.Ca' == 2
        @test s.Da' == 3
        @test s.Ab' == 4
        @test s.Bb' == 5
        @test s.Cb' == 6
        @test s.Db' == 7
    end

    @testset "parameter" begin
        @system STabulateParameter(Controller) begin
            T ~ tabulate(rows=(:A, :B, :C), columns=(:a, :b), parameter)
            Aa(T.A.a) ~ preserve
            Ba(T.B.a) ~ preserve
            Ca(T.C.a) ~ preserve
            Ab(T.A.b) ~ preserve
            Bb(T.B.b) ~ preserve
            Cb(T.C.b) ~ preserve
        end
        o = STabulateParameter => :T => [
            # a b
              0 3 ; # A
              1 4 ; # B
              2 5 ; # C
        ]
        s = instance(STabulateParameter, config=o)
        @test s.Aa' == 0
        @test s.Ba' == 1
        @test s.Ca' == 2
        @test s.Ab' == 3
        @test s.Bb' == 4
        @test s.Cb' == 5
    end

    @testset "default columns" begin
        @system STabulateDefaultColumns(Controller) begin
            T => [
                # a b
                  0 2 ; # a
                  1 3 ; # b
            ] ~ tabulate(rows=(:a, :b))
            aa(T) => T.a.a ~ preserve
            ba(T) => T.b.a ~ preserve
            ab(T) => T.a.b ~ preserve
            bb(T) => T.b.b ~ preserve
        end
        s = instance(STabulateDefaultColumns)
        @test getfield(s.T, :rows) == getfield(s.T, :columns) == (:a, :b)
        @test s.aa' == 0
        @test s.ba' == 1
        @test s.ab' == 2
        @test s.bb' == 3
    end

    @testset "arg" begin
        @system STabulateArg(Controller) begin
            x => 0 ~ preserve(parameter)
            T(x) => [
                # a b
                  x   x+2 ; # A
                  x+1 x+3 ; # B
            ] ~ tabulate(rows=(:A, :B), columns=(:a, :b))
        end
        s = instance(STabulateArg)
        @test s.T.A.a == 0
        @test s.T.B.a == 1
        @test s.T.A.b == 2
        @test s.T.B.b == 3
    end
end