text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Generate some fake data into datasets for grammar experiments
"""
module FakeDataGen

export generate_fake_data

import Compat.ASCIIString
using RLESUtils, DataFrameSets
using DataFrames
using Datasets

const DATAPATH = Pkg.dir("Datasets/data")
const BIN_SYNTH_NAME = "bin_synth"
const BIN_TS_SYNTH_NAME = "bin_ts_synth"
const VHDIST_NAME = "vhdist"
const VDIST_NAME = "vdist"

function generate_fake_data(dataset::AbstractString; kwargs...)
    generate_fake_data(Val{Symbol(dataset)}; kwargs...)
end

"""
bin_synth dataset
synthetic data with
static binary features
"""
function generate_fake_data(::Type{Val{Symbol(BIN_SYNTH_NAME)}};
    n_feats::Int64=20,
    n_samples::Int64=5000)

    mkpath(joinpath(DATAPATH, BIN_SYNTH_NAME))
    D = convert(Array{Int64}, rand(Bool, n_samples, n_feats))
    colnames = [Symbol("x$i") for i = 1:n_feats]
    feats = DataFrame(D)
    names!(feats, colnames)
    filename = joinpath(DATAPATH, BIN_SYNTH_NAME, "feats.csv.gz")
    writetable(filename, feats)

    #labels
    labels = DataFrame()
    labels[:x1] = feats[:x1] #x1
    labels[:x2] = feats[:x2] #x2
    labels[:x1_and_x3] = feats[:x1] & feats[:x3] #x1 and x3
    labels[:x2_or_x4] = feats[:x2] | feats[:x4] #x2 or x4
    labels[:x1_and_x4] = feats[:x1] & feats[:x4] #x1 and x4
    labels[:x2_or_x3] = feats[:x2] | feats[:x3] #x2 or x3
    filename = joinpath(DATAPATH, BIN_SYNTH_NAME, "labels.csv.gz")
    writetable(filename, labels)
end

"""
bin_ts_synth dataset
synthetic data with
binary time series features
"""
function generate_fake_data(::Type{Val{Symbol(BIN_TS_SYNTH_NAME)}};
    n_feats::Int64=20,
    n_time::Int64=50,
    n_samples::Int64=5000)

    outdir = joinpath(DATAPATH, BIN_TS_SYNTH_NAME)
    mkpath(outdir)

    records = DataFrame[]
    colnames = [Symbol("x$i") for i = 1:n_feats]
    for j = 1:n_samples
        B = rand(Bool, n_time, n_feats)
        for i = 1:size(B, 2)
            r = rand(1:3)
            if r == 1
                B[:,i] = false
            elseif r == 2
                B[:,i] = true
            end
        end
        A = convert(Array{Int64}, B)
        D = DataFrame(A)
        names!(D, colnames)
        push!(records, D)
    end
    Ds = DFSet(records)
    meta = getmeta(Ds)

    #labels, saved into meta df
    meta[:F_x1] = [reduce(|, r[:x1]) for r in records] #F(x1)
    meta[:G_x1] = [reduce(&, r[:x1]) for r in records] #G(x1)
    meta[:F_x2] = [reduce(|, r[:x2]) for r in records] #F(x2)
    meta[:G_x2] = [reduce(&, r[:x2]) for r in records] #G(x2)
    meta[:F_x3] = [reduce(|, r[:x3]) for r in records] #F(x3)
    meta[:G_x3] = [reduce(&, r[:x3]) for r in records] #G(x3)
    meta[:F_x4] = [reduce(|, r[:x4]) for r in records] #F(x4)
    meta[:G_x4] = [reduce(&, r[:x4]) for r in records] #G(x4)
    meta[:F_x5] = [reduce(|, r[:x5]) for r in records] #F(x5)
    meta[:G_x5] = [reduce(&, r[:x5]) for r in records] #G(x5)

    #F(x1 & x3)
    meta[:F_x1_and_x3] = [reduce(|, r[:x1] & r[:x3]) for r in records]
    #F(x1 | x3)
    meta[:F_x1_or_x3] = [reduce(|, r[:x1] | r[:x3]) for r in records]
    #G(x1 & x3)
    meta[:G_x1_and_x3] = [reduce(&, r[:x1] & r[:x3]) for r in records]
    #G(x1 | x3)
    meta[:G_x1_or_x3] = [reduce(&, r[:x1] | r[:x3]) for r in records]
    #F(x2 & x4)
    meta[:F_x2_and_x4] = [reduce(|, r[:x2] & r[:x4]) for r in records]
    #F(x2 | x4)
    meta[:F_x2_or_x4] = [reduce(|, r[:x2] | r[:x4]) for r in records]
    #G(x2 & x4)
    meta[:G_x2_and_x4] = [reduce(&, r[:x2] & r[:x4]) for r in records]
    #G(x2 | x4)
    meta[:G_x2_or_x4] = [reduce(&, r[:x2] | r[:x4]) for r in records]
    #F(x1 & x4)
    meta[:F_x1_and_x4] = [reduce(|, r[:x1] & r[:x4]) for r in records]
    #F(x1 | x4)
    meta[:F_x1_or_x4] = [reduce(|, r[:x1] | r[:x4]) for r in records]
    #G(x1 & x4)
    meta[:G_x1_and_x4] = [reduce(&, r[:x1] & r[:x4]) for r in records]
    #G(x1 | x4)
    meta[:G_x1_or_x4] = [reduce(&, r[:x1] | r[:x4]) for r in records]

    save_csvs(outdir, Ds)
    Ds
end

"""
vertical and horizontal miss distances from ACASX dataset
real-valued time series features
actually not fake data...
"""
function generate_fake_data(::Type{Val{Symbol(VHDIST_NAME)}};
    cas_data_name::ASCIIString="dasc",
    vdist_name::Symbol=:abs_alt_diff,
    hdist_name::Symbol=:horizontal_range,
    nmac_name::Symbol=:nmac,
    outname::ASCIIString=VHDIST_NAME
    )

    outdir = joinpath(DATAPATH, outname)
    mkpath(outdir)

    cas_data = dataset(cas_data_name)

    records = DataFrame[]
    colnames = [:vdist, :hdist]
    for r in getrecords(cas_data)
        D = r[[vdist_name, hdist_name]]
        names!(D, colnames)
        push!(records, D)
    end

    Ds = DFSet(getmeta(cas_data), records)
    save_csvs(outdir, Ds)
end

"""
vertical distances from ACASX dataset
real-valued time series features
actually not fake data...
"""
function generate_fake_data(::Type{Val{Symbol(VDIST_NAME)}};
    cas_data_name::ASCIIString="dasc",
    vdist_name::Symbol=:abs_alt_diff,
    nmac_name::Symbol=:nmac,
    nmac_thresh::Float64=100.0,
    outname::ASCIIString=VDIST_NAME
    )

    outdir = joinpath(DATAPATH, outname)
    mkpath(outdir)

    cas_data = dataset(cas_data_name)
    meta = metadf(length(cas_data))

    records = DataFrame[]
    nmac = Bool[]
    for r in getrecords(cas_data)
        D = DataFrame()
        D[:vdist] = r[vdist_name]
        push!(records, D)
        push!(nmac, any(D[:vdist] .< nmac_thresh))
    end

    meta[nmac_name] = nmac
    Ds = DFSet(meta, records)
    save_csvs(outdir, Ds)
end

end #module


text:
#using PlotlyJS
using Plots
using UnicodePlots
using Crayons
using ProgressBars

# Arranged (y, x) not (x, y)
"""
Gets all the moves the player (`to_move`) can make in a `position`.
"""
function get_moves(position::Array, to_move::Int)
    moves = []

    for y in 1:length(position)
        for x in 1:length(position[y])

            if position[y][x] == to_move

                if isassigned(position, y + to_move)
                    if position[y + to_move][x] == 0
                        push!(moves, [[y, x], [y + to_move, x]])
                    end # if

                    if isassigned(position[y+to_move], x + 1)
                        if position[y + to_move][x + 1] == to_move * -1
                            push!(moves, [[y, x], [y + to_move, x + 1]])
                        end # if
                    end # if
                    if isassigned(position[y+to_move], x - 1)
                        if position[y + to_move][x - 1] == to_move * -1
                            push!(moves, [[y, x], [y + to_move, x - 1]])
                        end # if
                    end # if
                end # if

            end # if

        end #for
    end # for

    return moves

end # function

"""
Returns a position which is `pos` with the `move` made. The returned array is a deepcopy.
"""
function new_p(pos::Array, move::Array)
    _pos = deepcopy(pos)
    _pos[move[2][1]][move[2][2]] = _pos[move[1][1]][move[1][2]]
    _pos[move[1][1]][move[1][2]] = 0

    return _pos
end # function

"""
Gets the index of some `value` in an `array`. Returns `0` if the value does not exist in the array.
"""
function extract(value::Any, array::Array)
    for i in 1:length(array)
        if array[i] == value
            return i
        end # if
    end # for
    return 0
end # function

"""
Checks to see if anyone has won or if the game is draw. Returns `"live"` if the game is ongoing, and `0` if the game is a draw.
"""
function check_victors(pos::Array)

    if -1 in pos[1]
        return -1
    elseif 1 in pos[end]
        return 1
    elseif length(get_moves(pos, 1)) == 0 && length(get_moves(pos, -1)) == 0
        return 0
    elseif length(get_moves(pos, 1)) == 0
        return -1
    elseif length(get_moves(pos, -1)) == 0
        return 1
    end # if

    return "live"

end

function reflect(pos::Array)
    r = [Int[], Int[], Int[]]
    for row in 1:length(pos)
        r[row] = [pos[row][end - sq + 1] for sq in 1:length(pos[row])]
    end # for
end # function

"""
Use recursive DFS to create a set of links and positions which covers every board state possible in hexapawn and links them together.
"""
function dfs!(pos::Array, links::Array, seen::Array; player = 1)

    for m in get_moves(pos, player)

        if extract(new_p(pos, m), seen) == 0
            # println("Checking: " * string(new_p(pos, move)))
            # println("Links: " * string(links))

            append!(seen, [new_p(pos, m)])
            push!(links, Int[])

            if check_victors(new_p(pos, m)) == "live"
                dfs!(new_p(pos, m), links, seen, player=player*-1)
            end # If
        end # if

        append!(links[extract(pos, seen)], extract(new_p(pos, m), seen))

    end # for

end # function

#== SELF PLAY ==#

"""
Using `links` to get a new `pos` using the weights from `w`.
"""
function move(pos::Int, w::Array, links::Array)

    _move = [0, 0]
    if sum(w[pos]) > 0
        bead = rand(1:sum(w[pos]))
        i = 1
        for m in w[pos]
            if bead > m
                bead -= m
            else
                _move[1] = pos
                pos = links[pos][i]
                _move[2] = i
                break
            end # if
            i += 1
        end # for
    else
        pos = links[pos][1]
    end # if

    return [pos, _move]

end # function

"""
Pits `w` against `b` for one game. Uses `links` for connections. Returns `[w_moves, b_moves, victor]` (victor according to check_victors()).
If the game goes over 200 moves the function returns `"ERROR"`. Real games of hexapawn do not last very long.
"""
function pit(w::Array, b::Array, links::Array)

    current_pos = 1
    w_moves = []
    b_moves = []

    for i in 1:200

        if check_victors(positions[current_pos]) != "live"
            return [w_moves, b_moves, check_victors(positions[current_pos])]
        end # if

        out = move(current_pos, w, links)
        current_pos = out[1]
        push!(w_moves, out[2])

        if check_victors(positions[current_pos]) != "live"
            return [w_moves, b_moves, check_victors(positions[current_pos])]
        end # if

        out = move(current_pos, b, links)
        current_pos = out[1]
        push!(b_moves, out[2])

    end # for

    println("Got an error")
    println("Position #: " * string(current_pos))
    println("Position: " * string(positions[current_pos]))
    println("Links: " * string(links[current_pos]))
    return "ERROR"

end # function

"""
Trains the weights in `w` and `b` for `games` games. Bonuses and punishments can be set using `win`, `lose`, and `draw`.
"""
function train!(w::Array, b::Array; games = 200, win = 3, lose = -1, draw = 1, do_plot = false)

    data_to_plot = [[] for l in w]
    println("Training")

    for i in ProgressBar(1:games)

        out = pit(w, b, links)

        if out != "ERROR"

            w_moves = out[1]
            b_moves = out[2]
            result = out[3]

            # edit weights based on results
            if result == 1 # W WINS
                for m in w_moves

                    w[m[1]][m[2]] += win

                end # for
                for m in b_moves

                    b[m[1]][m[2]] += lose

                end # for
            elseif result == -1 # B WINS
                for m in w_moves

                    w[m[1]][m[2]] += lose

                end # for
                for m in b_moves

                    b[m[1]][m[2]] += win

                end # for
            else # DRAWS
                for m in w_moves

                    w[m[1]][m[2]] += draw

                end # for
                for m in b_moves

                    b[m[1]][m[2]] += draw

                end # for
            end # if

            # make sure there is at least one bead in the box
            for m in w_moves
                if sum(w[m[1]]) <= 0
                    w[m[1]][1] == 1
                end # if
            end # for
            for m in b_moves
                if sum(b[m[1]]) <= 0
                    b[m[1]][1] == 1
                end # if
            end # for

            if do_plot

                for d in w
                    append!(data_to_plot[extract(d, w)], sum(d))
                end # for

            end # if

        end # if

    end # for

    if do_plot

        # plt = lineplot([i for i in 1:length(data_to_plot[1])], [x for x in data_to_plot[1]], name = "box 1", xlabel = "Games", ylabel = "Beads")
        # println(plt)

        filtered = sort([[extract(u, w), sum(u)] for u in w if sum(u) / length(u) != 5 && sum(u) != 0], by= x -> x[2], rev = true)[1:20]
        plt = barplot([string(i[1]) for i in filtered], [x[2] for x in filtered], title = "Sum in box")
        println(plt)

        plt = scatter([y for y in data_to_plot[1]])
        Plots.savefig(plt, "box1")

        #plt = lineplot([i for i in 1:length(data_to_plot[2])], [x for x in data_to_plot[2]], name = "box 2", xlabel = "Games")
        #for box in filtered[1:10]
        #    lineplot!(plt, [i for i in 1:length(data_to_plot[box[1]])], data_to_plot[box[1]], name = "box " * string(box[1]))
        #end # for
        #println(plt)

    end # if

end # function

#== ALLOW THE HUMAN TO PLAY IT ==#

"""
Takes a move typed by someone in the format: `Y X, Y X` where 1 1 is the top left.
"""
function human_move(raw::String)

    move = [[0, 0], [0, 0]]
    move[1][1] = parse(Int, string(raw[1]))
    move[1][2] = parse(Int, string(raw[3]))
    move[2][1] = parse(Int, string(raw[6]))
    move[2][2] = parse(Int, string(raw[8]))

    return move

end # function

"""
Prints out the current position without all those ugly brackets.
"""
function format(pos::Array)

    formatted = [[] for i in 1:length(pos)]
    for row in 1:length(pos)
        for square in pos[row]
            if square == 0
                append!(formatted[row], "-")
            elseif square == 1
                append!(formatted[row], "W")
            elseif square == -1
                append!(formatted[row], "B")
            end # if
        end # for
    end # for

    for row in formatted
        for square in row
            print(Crayon(reset = true), square * " ")
        end # for
        println()
    end # for

end # function

"""
Does all the pretty human stuff to let humans play a game against the alg.
"""
function play(w::Array, links::Array, positions::Array; base = [[1, 1, 1], [0, 0, 0], [-1, -1, -1]])

    println("You play as W/B? ")
    human = uppercase(readline())
    pos = deepcopy(base)

    blue = (150, 150, 255)
    green = (100, 255, 100)

    for i in 1:100

        println()
        println(Crayon(foreground = green, bold = true), "Move " * string(i))

        format(pos)

        if human == "W"

            println(Crayon(foreground = blue, bold = false), "Your move? ")
            _move = human_move(readline())
            println(_move)
            println()
            pos = new_p(pos, _move)
            format(pos)

            if check_victors(pos) != "live"
                println(Crayon(foreground = green, bold = true), "The game has ended")
                break
            end # if

        end # if

        println()
        println("ML to move...")
        c_move = move(extract(pos, positions), w, links)
        pos = positions[c_move[1]]
        format(pos)

        if check_victors(pos) != "live"
            println(Crayon(foreground = green, bold = true), "The game has ended")
            break
        end # if

        if human == "B"

            println(Crayon(foreground = blue, bold = false), "Your move? ")
            _move = human_move(readline())
            println(_move)
            println()
            pos = new_p(pos, _move)
            format(pos)

            if check_victors(pos) != "live"
                println(Crayon(foreground = green, bold = true), "The game has ended")
                break
            end # if

        end # if

    end # for

end # function

#== If you wanted to fiddle with the additional params, chnaces are you only need to touch what's below this line ==#

# Setup
println("Dimension to train for (max 9)? [y,x]")
dims = readline()

println("How many games to self-play?")
games_to = parse(Int, readline())

base = []
for i in 1:parse(Int, dims[1])
    push!(base, Int[])
    for x in 1:parse(Int, dims[3])
        if i == 1
            push!(base[i], 1)
        elseif i == parse(Int, dims[1])
            push!(base[i], -1)
        else
            push!(base[i], 0)
        end # if
    end # for
end # for
println(base)

player = 1
positions = [base]
links = [[]]

# DFS
dfs!(base, links, positions)
println("Got " * string(length(positions)) * " positions")

#==
for l in 1:length(links)
    println(string(l) * " : " * string(links[l]))
end # for
==#

# Setup weights
weights = [Int[2*length(base) for m in l] for l in links]

# Self-play
train!(weights, weights, games=games_to, do_plot=true)

# Human play
while true
    play(weights, links, positions, base=base)
    println(Crayon(reset = true), "")
end # while

ext: jl

max_stars_repo_name: Squalm/MENACE

max_issues_repo_name: Squalm/MENACE

max_forks_repo_name: Squalm/MENACE


text:
abstract type AbstractVgg end
struct Vgg19  <: AbstractVgg end
struct Vgg416 <: AbstractVgg end

vgg19() = UnstableChain([
        # First Layer
        Conv((3,3), 3 => 64, relu; pad = 1),
        Conv((3,3), 64 => 64, relu; pad = 1),
        MaxPool((2,2)),
        # Second Layer
        Conv((3,3), 64 => 128, relu; pad = 1),
        Conv((3,3), 128 => 128, relu; pad = 1),
        MaxPool((2,2)),
        # Third Layer
        Conv((3,3), 128 => 256, relu; pad = 1),
        Conv((3,3), 256 => 256, relu; pad = 1),
        Conv((3,3), 256 => 256, relu; pad = 1),
        Conv((3,3), 256 => 256, relu; pad = 1),
        MaxPool((2,2)),
        # Fourth Layer
        Conv((3,3), 256 => 512, relu; pad = 1),
        Conv((3,3), 512 => 512, relu; pad = 1),
        Conv((3,3), 512 => 512, relu; pad = 1),
        Conv((3,3), 512 => 512, relu; pad = 1),
        MaxPool((2,2)),
        # Fifth Layer
        Conv((3,3), 512 => 512, relu; pad = 1),
        Conv((3,3), 512 => 512, relu; pad = 1),
        Conv((3,3), 512 => 512, relu; pad = 1),
        Conv((3,3), 512 => 512, relu; pad = 1),
        MaxPool((2,2)),
        # Fully Connected Layers
        x -> reshape(x, :, size(x, 4)),
        Dense(25088, 4096, relu),
        Dense(4096, 4096, relu),
        Dense(4096, 1000, relu),
        # Truncate to a small positive value to avoid NaNs
        x -> log.(max.(x, Float32(1e-7))),
        softmax
    ])

function vgg416()
    loops = (80, 81, 82, 83, 83)
    #loops = (40, 41, 42, 43, 43)
    layers = []
    # First Layer
    push!(layers, Conv((3,3), 3 => 64, relu; pad = 1))
    for _ in 1:loops[1]
    push!(layers, Conv((3,3), 64 => 64, relu; pad = 1))
    end
    push!(layers,
        Conv((3,3), 64 => 128, relu; pad = 1),
        MaxPool((2,2))
    )

    # Second Layer
    for _ in 1:loops[2]
        push!(layers, Conv((3,3), 128 => 128, relu; pad = 1))
    end
    push!(layers,
        Conv((3,3), 128 => 256, relu; pad = 1),
        MaxPool((2,2))
    )
    # Third Layer
    for _ in 1:loops[3]
        push!(layers, Conv((3,3), 256 => 256, relu; pad = 1))
    end
    push!(layers,
        Conv((3,3), 256 => 512, relu; pad = 1),
        MaxPool((2,2)),
    )

    # Fourth Layer
    for _ in 1:loops[4]
        push!(layers, Conv((3,3), 512 => 512, relu; pad = 1))
    end
    push!(layers, MaxPool((2,2)))
    # Fifth Layer
    for _ in 1:loops[5]
        push!(layers, Conv((3,3), 512 => 512, relu; pad = 1))
    end
    push!(layers,
        MaxPool((2,2)),
        # Fully Connected Layers
        x -> reshape(x, :, size(x, 4)),
        Dense(25088, 4096, relu),
        Dense(4096, 4096, relu),
        Dense(4096, 1000),
        # Add a small positive value to avoid NaNs
        x -> max.(x, Float32(1e-7)),
        softmax
    )

    # Return a vetor version of a chain
    return UnstableChain(layers)
end

function vgg19_inference(batchsize)
    x = rand(Float32, 224, 224, 3, batchsize)

    x = (x .- mean(x)) ./ std(x)

    forward = vgg19()
    return Actualizer(forward, x)
end

_forward(::Vgg19) = vgg19()
_forward(::Vgg416) = vgg416()

function vgg_training(vgg::T, batchsize) where {T <: AbstractVgg}
    X = randn(Float32, 224, 224, 3, batchsize)
    Y = zeros(Float32, 1000, batchsize)
    for col in 1:batchsize
        Y[rand(1:1000), col] = one(eltype(Y))
    end

    # Get the forward pass
    forward = _forward(vgg)

    f = ForwardLoss(_forward(vgg), Flux.crossentropy)
    return Actualizer(f, X, Y; optimizer = nGraph.SGD(Float32(0.05)))
end

function random_labels!(y::AbstractArray{T,2}) where {T}
    y .= zero(eltype(y))
    for j in 1:size(y, 2)
        y[rand(1:size(y, 1)), j] = one(eltype(y))
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
module Game2048

export Dirs, left, right, up, down, move!, move

@enum Dirs left right up down

include("uint_board/move-up.jl")
# include("gamesim/move-board.jl")
# include("gamesim/sim-game.jl")
# include("gamesim/play_game_with_policy.jl")

# include("identity-skip.jl")

include("rotate-mirror.jl")

include("bitboard/make_lookups.jl")
include("bitboard/bitboard.jl")
include("bitboard/simulate_bb.jl")
# include("bitboard/playahead.jl")
# include("bitboard/fingerprint.jl")

# include("monte-carlo-tree.jl")

# include("attempts2020/5_policy_gradient.jl")
end # module


text:
@testset "Quadtree.jl" begin
    width, height = 800, 800
    origin = (width/2, height/2)

    @init "Quadtree test 1" width height
    @layout aside(:v, 250)

    @create @options begin
        capacity = (4, 1:1000)
        search_radius = (200, 1:800)
        reset = Button("Reset", @hexpand)
    end

    Quadtrees.position(p::Point) = (p.x, p.y)

    pointlist = [randompoint(0, 0, width, height) for i in 1:500]

    mouse = Ref{Point}()

    onevent(:clicked, reset) do
        empty!(pointlist)
    end

    @use function update()
        k = floor(Int, capacity[])
        qtree = Quadtree{Point}(origin, width, height, k)
        append!(qtree, pointlist)

        background("black")
        sethue("white")
        draw(qtree, items = true)

        if !ismissing(mouse)
            x, y = mouse[]
            r = search_radius[]
            setcolor"yellowgreen"
            circle(x, y, r, :stroke)
            draw.(search(qtree, x, y, r), 4)
        end
    end

    @use function mousemove(event)
        mouse[] = event.pos
    end

    @use function mouse1motion((pos, x, y))
        push!(pointlist, pos)
    end

    start()
end


text:
module TestReduce
using DataArrays, Base.Test, StatsBase, Compat

srand(1337)

## extended test of sum

for skipna in (true, false)
    if VERSION < v"0.5-"
        @test sum(@data(Int8[]); skipna=skipna) === 0
        @test sum(@data(Int[]); skipna=skipna) === 0
        @test sum(@data(Float64[]); skipna=skipna) === 0.0

        @test sum(@data([@compat(Int8(3))]); skipna=skipna) === 3
        @test sum(@data([3]); skipna=skipna) === 3
        @test sum(@data([3.0]); skipna=skipna) === 3.0

        z = DataArray(reshape(1:16, (2,2,2,2)))
        fz = convert(DataArray{Float64}, z)
        bfz = convert(DataArray{BigFloat}, z)
        @test sum(z) === 136
        @test sum(fz) === 136.0
        @test sum(bfz) == 136
    else
        @test sum(@data(Int8[]); skipna=skipna) === Int32(0)
        @test sum(@data(Int[]); skipna=skipna) === 0
        @test sum(@data(Float64[]); skipna=skipna) === 0.0

        @test sum(@data([@compat(Int8(3))]); skipna=skipna) === Int32(3)
        @test sum(@data([3]); skipna=skipna) === 3
        @test sum(@data([3.0]); skipna=skipna) === 3.0

        z = DataArray(reshape(1:16, (2,2,2,2)))
        fz = convert(DataArray{Float64}, z)
        bfz = convert(DataArray{BigFloat}, z)
        @test sum(z) === 136
        @test sum(fz) === 136.0
        @test sum(bfz) == 136
    end
end

@test sum(@data(Int[NA])) === NA
@test sum(@data(Int[NA]); skipna=true) === 0
@test sum(@data(Int[NA, NA])) === NA
@test sum(@data(Int[NA, NA]); skipna=true) === 0
@test sum(@data(Int[NA, NA, 1]); skipna=true) === 1
@test sum(@data(Int[NA, NA, 1, 2]); skipna=true) === 3
@test sum(@data(Int[NA, 1, NA, 1, 2]); skipna=true) === 4

z = DataArray(reshape(1:16, (2,2,2,2)))
z[6] = NA
fz = convert(DataArray{Float64}, z)
bfz = convert(DataArray{BigFloat}, z)
@test isna(sum(z))
@test isna(sum(fz))
@test isna(sum(bfz))
@test sum(z; skipna=true) === 130
@test sum(fz; skipna=true) === 130.0
@test sum(bfz; skipna=true) == 130

bs = DataArrays.sum_pairwise_blocksize(@functorize(identity))
for n in [bs-64, bs-1, bs, bs+1, bs+2, 2*bs-2:2*bs+3..., 4*bs-2:4*bs+3...]
    da = DataArray(randn(n))
    s = sum(da.data)
    @test_approx_eq sum(da) s
    @test_approx_eq sum(da; skipna=true) s

    da2 = copy(da)
    da2[1:2:end] = NA
    @test isna(sum(da2))
    @test_approx_eq sum(da2; skipna=true) sum(dropna(da2))

    da2 = convert(DataArray{BigFloat}, da2)
    @test isna(sum(da2))
    @test_approx_eq sum(da2; skipna=true) sum(dropna(da2))

    da2 = copy(da)
    da2[2:2:end] = NA
    @test isna(sum(da2))
    @test_approx_eq sum(da2; skipna=true) sum(dropna(da2))

    da2 = convert(DataArray{BigFloat}, da2)
    @test isna(sum(da2))
    @test_approx_eq sum(da2; skipna=true) sum(dropna(da2))
end

## other reductions

macro same_behavior(ex1, ex2)
    quote
        v = try
            $ex2
        catch e
            e
        end
        isa(v, Exception) ? @test_throws(typeof(v), $ex1) : @test_approx_eq($ex1, v)
    end
end

_varuc(x; kw...) = var(x; corrected=false, kw...)
_varzm(x; kw...) = var(x; mean=0, kw...)
_varzmuc(x; kw...) = var(x; corrected=false, mean=0, kw...)
_var1m(x; kw...) = var(x; mean=1, kw...)
_var1muc(x; kw...) = var(x; corrected=false, mean=1, kw...)
_std1m(x; kw...) = stdm(x, 1; kw...)

for fn in (prod, minimum, maximum, mean,
           var, _varuc, _varzm, _varzmuc, _var1m, _var1muc,
           std, _std1m, Base.sumabs, Base.sumabs2)
    for n in [0, 1, 2, 62, 63, 64, 65, 66]
        da = DataArray(randn(n))
        @same_behavior fn(da) fn(da.data)
        @same_behavior fn(da; skipna=true) fn(da.data)

        da2 = copy(da)
        da2[1:2:end] = NA
        n > 0 && @test isna(fn(da2))
        @same_behavior fn(da2; skipna=true) fn(dropna(da2))

        da2 = convert(DataArray{BigFloat}, da2)
        n > 0 && @test isna(fn(da2))
        @same_behavior fn(da2; skipna=true) fn(dropna(da2))

        da2 = copy(da)
        da2[2:2:end] = NA
        n > 1 && @test isna(fn(da2))
        @same_behavior fn(da2; skipna=true) fn(dropna(da2))

        da2 = convert(DataArray{BigFloat}, da2)
        n > 1 && @test isna(fn(da2))
        @same_behavior fn(da2; skipna=true) fn(dropna(da2))
    end
end

## reduce and mapreduce drivers

for fn in (+, *, |, &)
    da = convert(DataArray, bitrand(10))

    s = mapreduce(@functorize(identity), fn, da.data)
    @test mapreduce(@functorize(identity), fn, da) == s
    @test mapreduce(@functorize(identity), fn, da; skipna=true) == s
    @test reduce(fn, da) == s
    @test reduce(fn, da; skipna=true) == s
end

# make sure reductions of & and | are still calling Base
@test isna(reduce(&, @data([true, NA])))
@test !reduce(&, @data([false, NA]))
@test reduce(|, @data([true, NA]))
@test isna(reduce(|, @data([false, NA])))

# weighted mean
da1 = DataArray(randn(128))
da2 = DataArray(randn(128))
@same_behavior mean(da1, weights(da2)) mean(da1.data, weights(da2.data))
@same_behavior mean(da1, weights(da2.data)) mean(da1.data, weights(da2.data))
@same_behavior mean(da1, weights(da2); skipna=true) mean(da1.data, weights(da2.data))
@same_behavior mean(da1, weights(da2.data); skipna=true) mean(da1.data, weights(da2.data))

da1[1:3:end] = NA
@same_behavior mean(da1, weights(da2); skipna=true) mean(dropna(da1), weights(da2.data[!da1.na]))
@same_behavior mean(da1, weights(da2.data); skipna=true) mean(dropna(da1), weights(da2.data[!da1.na]))

da2[1:2:end] = NA
keep = !da1.na & !da2.na
@test isna(mean(da1, weights(da2)))
@same_behavior mean(da1, weights(da2); skipna=true) mean(da1.data[keep], weights(da2.data[keep]))
end


text:
input_nums = zeros(Bool, 2021)

open("day01/day01.in") do io
    for line in readlines(io)
        num = parse(Int32, line)
        input_nums[parse(Int32, line)] = true
    end
end

for num1 in 1:1010
    for num2 in 1:num1
        num3 = 2020 - num1 - num2
        if num3 > 0 && input_nums[num1] && input_nums[num2] && input_nums[num3]
            println(num1 * num2 * num3)
        end
    end
end

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
# The MIT License (MIT)
#
# Copyright (c) 2016 Alexander Samoilov
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of
# this software and associated documentation files (the "Software"), to deal in
# the Software without restriction, including without limitation the rights to
# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
# the Software, and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
# COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Ld = 365.
NA = 0.42
Mask_L = 1 # um

x = y = flipdim((-256:255) * (Ld / (2.*NA)*31) / 512, 2)

tm = zeros(length(x))

max_stars_repo_name: alsam/julia-samples

max_issues_repo_name: alsam/julia-samples

max_forks_repo_name: alsam/julia-samples


text:
include("day5_solution.jl")
using Test

@test parse_rows("FBFBBFF") == 44
@test parse_columns("RLR") == 5

@test get_seat_id("BFFFBBFRRR") == 567
@test get_seat_id("FFFBBBFRRR") == 119
@test get_seat_id("BBFFBBFRLL") == 820

@test id_to_row_column(357) == (44, 5)


text:
using Yao, YaoExtensions
using YaoZX
using Test

@testset "decompose gates" begin
	@test YaoZX.getclocs(cnot(2,-2,1)) == (-2,)
	@test YaoZX.getclocs(control(4, (-2, 3), 1=>X)) == (-2, 3)
	for g in [control(5, (2, 1), 4=>X),
		 	put(5, 3=>ConstGate.T), put(5, 2=>ConstGate.Sdag),
			put(5, 4=>ConstGate.Tdag), put(5, 2=>ConstGate.S),
			control(5, 4, 3=>Y), put(5, 3=>shift(0.4)),
			control(5, 2, 3=>Ry(0.5)),
			control(5, 2, 3=>Rz(0.5)),
			control(5, 2, 3=>Rx(0.5)),
			control(5, -2, 3=>Rx(0.5)),
			cphase(5, 2, 4, 0.5)
			]
		@test operator_fidelity(decompose_zx(g), g) ≈ 1
	end
end


text:
type MidPointCdsEngine{P <: AbstractDefaultProbabilityTermStructure, Y <: YieldTermStructure} <: PricingEngine
  probability::P
  recoveryRate::Float64
  discountCurve::Y
end

function _calculate!(pe::MidPointCdsEngine, swap::CreditDefaultSwap)
  today_ = settings.evaluation_date
  settlementDate = reference_date(pe.discountCurve)

  # Upfront flow NPV, either we are on the run (no flow) or we are forward start
  upFPV01 = 0.0
  if has_occurred(swap.upfrontPayment, settlementDate)
    # date determining the probability survival so we have to pay the upfront (did not knock out)
    effectiveUpfrontDate = swap.protectionStart > reference_date(pe.probability) ? swap.protectionStart : reference_date(pe.probability)

    upFPV01 = survival_probability(pe.probability, effectiveUpfrontDate) * discount(pe.discountCurve, date(swap.upfrontPayment))
  end

  swap.results.upfrontNPV = upFPV01 * amount(swap.upfrontPayment)
  swap.results.couponLegNPV = 0.0
  swap.results.defaultLegNPV = 0.0

  for i in eachindex(swap.leg.coupons)
    if has_occurred(swap.leg.coupons[i], settlementDate)
      continue
    end

    coup = swap.leg.coupons[i]

    # first calculate NPV of both legs as positive, applying correct sign at the end
    paymentDate = date(coup)
    startDate = accrual_start_date(coup)
    endDate = accrual_end_date(coup)

    if i == 1
      startDate = swap.protectionStart
    end

    effectiveStartDate = (startDate <= today_ && today_ <= endDate) ? today_ : startDate
    defaultDate = effectiveStartDate + Dates.Day(floor(Int(endDate - effectiveStartDate) / 2)) # midpoint

    S = survival_probability(pe.probability, paymentDate)
    P = default_probability(pe.probability, effectiveStartDate, endDate)

    # on one side, we add the fixed rate payments in case of survival
    swap.results.couponLegNPV += S * amount(coup) * discount(pe.discountCurve, paymentDate)

    # possibly including accrual in case of default
    if swap.settlesAccrual
      if swap.paysAtDefaultTime
        swap.results.couponLegNPV += P * accrued_amount(coup, defaultDate) * discount(pe.discountCurve, defaultDate)
      else
        # pays at end
        swap.results.couponLegNPV += P * amount(coup) * discount(pe.discountCurve, paymentDate)
      end
    end

    # on the other side we add the payment in case of default
    claim = amount(swap.claim, defaultDate, swap.notional, pe.recoveryRate)

    if swap.paysAtDefaultTime
      swap.results.defaultLegNPV += P * claim * discount(pe.discountCurve, defaultDate)
    else
      swap.results.defaultLegNPV += P * claim * discount(pe.discountCurve, paymentDate)
    end

  end

  upfrontSign = 1.0
  if isa(swap.side, Seller)
    swap.results.defaultLegNPV *= -1.0
  else
    swap.results.couponLegNPV *= -1.0
    swap.results.upfrontNPV *= -1.0
    upfrontSign = -1.0
  end

  swap.results.value = swap.results.defaultLegNPV + swap.results.couponLegNPV + swap.results.upfrontNPV

  if swap.results.couponLegNPV != 0.0
    swap.results.fairSpread = -swap.results.defaultLegNPV * swap.spread / swap.results.couponLegNPV
  else
    swap.results.fairSpread = -1.0
  end

  upfrontSensitivity = upFPV01 * swap.notional
  if upfrontSensitivity != 0.0
    swap.results.fairUpfront = -upfrontSign * (swap.results.defaultLegNPV + swap.results.couponLegNPV) / upfrontSensitivity
  else
    swap.results.fairUpfront = -1.0
  end

  basisPoint = 1.0e-4
  if swap.spread != 0.0
    swap.results.couponLegBPS = swap.results.couponLegNPV * basisPoint / swap.spread
  else
    swap.results.couponLegBPS = -1.0
  end

  # need to factor in swap upfront
  swap.results.upfrontBPS = -1.0

  return swap
end


text:
struct UnionVector{
    T,
    ETS <: NTypes,
    TD <: AbstractVector,
    TM <: AbstractVector{UInt8},
    TV <: Tuple
} <: Abstract.UnionVector{T}

    types::ETS
    data::TD
    typeid::TM
    views::TV

    function UnionVector(
        types::ETS,
        data::TD,
        typeid::TM,
        views::TV,
    ) where {
        ETS <: NTypes,
        TD <: AbstractVector,
        TM <: AbstractVector{UInt8},
        TV <: Tuple,
    }
        A = new{asunion(types),ETS,TD,TM,TV}(types, data, typeid, views)
        @static if VERSION >= v"1.6-"
            return verify(A)
        else
            return A
        end
    end
end

if VERSION >= v"1.6-"
    Adapt.adapt_structure(to, A::UnionVector) =
        UnionVector(A.types, Adapt.adapt(to, A.data), Adapt.adapt(to, A.typeid))
else
    Adapt.adapt_structure(to, A::UnionVector) = UnionVector(
        A.types,
        Adapt.adapt(to, A.data),
        Adapt.adapt(to, A.typeid),
        Adapt.adapt(to, A.views),
    )
end

executor_type(A::UnionVector) = executor_type(A.data)

function UnionVector(ETS::ElTypeSpec, data::AbstractVector, typeid::AbstractVector{UInt8})
    types = asntypes(ETS)
    views = foldlargs((), types...) do views, v
        return (views..., reinterpret(ofsamesize(eltype(data), valueof(v)), data))
    end
    return UnionVector(types, data, typeid, views)
end

UnionArrays.buffereltypefor(ETS::ElTypeSpec) =
    foldlargs(Nothing, asntypes(ETS)...) do S, x
        Base.@_inline_meta
        T = valueof(x)
        if sizeof_aligned(S) < sizeof_aligned(T)
            T
        else
            S
        end
    end

UnionArrays.buffereltypeof(A::UnionVector) = eltype(A.data)

function UnionArrays.eltypebyid(A::UnionVector, ::Val{id}) where {id}
    id > length(A.types) && return Union{}
    id < 1 && return Union{}
    return valueof(A.types[id])
end

Base.size(A::UnionVector) = size(A.data)

UnionVector(undef::UndefInitializer, ETS::ElTypeSpec, n::Integer) =
    UnionVector(undef, Vector, ETS, n)

UnionVector(
    undef::UndefInitializer,
    VectorType::Type{<:AbstractVector},
    ETS::ElTypeSpec,
    n::Integer,
) = UnionVector(undef, VectorType, VectorType{UInt8}, ETS, n)

function UnionVector(
    ::UndefInitializer,
    DataVectorType::Type{<:AbstractVector},
    TypeTagVectorType::Type{<:AbstractVector{UInt8}},
    ETS::ElTypeSpec,
    n::Integer,
)
    typeid = TypeTagVectorType(undef, n)
    fill!(typeid, 0)
    BT = buffereltypefor(ETS)
    data = DataVectorType{BT}(undef, n)
    return UnionVector(ETS, data, typeid)
end

function UnionVector(ETS::ElTypeSpec, items::AbstractVector)
    A = UnionVector(undef, ETS, length(items))
    for i in axes(A, 1)
        # TODO: find the best way to support arrays with offset
        @inbounds A[i] = items[i - firstindex(A) + firstindex(items)]
    end
    return A
end

function UnionVector(data::AbstractVector)
    ETS = foldl(data, init=Union{}) do ETS, x
        Union{ETS,typeof(x)}
    end
    return UnionVector(ETS, data)
end

struct InvalidUnionVector{T<:UnionVector} <: Exception
    A::T
end

function Base.showerror(io::IO, err::InvalidUnionVector)
    A = err.A
    unmatches = filter!(collect(Pair{Int,Any}, pairs(A.views))) do (_, xs)
        UInt(pointer(A.data)) != UInt(pointer(xs))
    end
    print(io, "invalid UnionVector: ")
    println(io, length(unmatches), " incompatible views")
    println(io, "data pointer: ", pointer(A.data))
    for (i, xs) in unmatches
        println(io, i, "-th view pointer: ", pointer(xs))
    end
end

function verify(A::UnionVector)
    foldlargs(true, A.views...) do ok, xs
        Base.@_inline_meta
        ok && UInt(pointer(A.data)) == UInt(pointer(xs))
    end || throw(InvalidUnionVector(A))
    return A
end

Base.@propagate_inbounds typeat(A::UnionVector{<:Any, ETS}, i) where ETS =
    fieldtype(ETS, Int(A.typeid[i]))

struct TypeIDLookupFailed <: Exception
    id::Int
end

Base.showerror(io::IO, err::TypeIDLookupFailed) =
    print(io, "lookup of type id $(err.id) failed")

# Using CPS for begging the compiler to union split things:
@inline function view_by_id(f::F, A::UnionVector, id::Integer) where {F}
    @noinline unreachable() = throw(TypeIDLookupFailed(id))
    return terminating_foldlargs(unreachable, 1, A.views...) do j, xs
        Base.@_inline_meta
        if j == id
            Reduced(f(xs))
        else
            j + 1
        end
    end
end

@inline function Base.getindex(A::UnionVector, i::Int)
    @boundscheck checkbounds(A, i)
    return view_by_id(A, A.typeid[i]) do xs
        Base.@_inline_meta
        unpad(@inbounds xs[i])
    end
end

struct ElTypeLookupFailed{T} <: Exception end
Base.showerror(io::IO, ::ElTypeLookupFailed{T}) where {T} =
    print(io, "unsupported element type (type conversion not implemented yet); $T")

# Using CPS for begging the compiler to union split things:
# TODO: handle conversion
@inline function view_and_id(f::F, A::UnionVector, ::Type{T}) where {F,T}
    @noinline unreachable() = throw(ElTypeLookupFailed{T}())
    V = Val(T)
    return terminating_foldlargs(unreachable, 1, A.views...) do id, xs
        Base.@_inline_meta
        if paddedtype(eltype(xs)) === valueof(V)
            Reduced(f(xs, id))
        else
            id + 1
        end
    end
end

# Base.checkbounds(::Type{Bool}, A::UnionVector, i) =
#     checkbounds(Bool, A.data, i) && checkbounds(Bool, A.typeid, i)

@inline function Base.setindex!(A::UnionVector, v, i::Int)
    @boundscheck checkbounds(A, i)
    typeid = A.typeid
    view_and_id(A, typeof(v)) do xs, id
        Base.@_inline_meta
        @inbounds typeid[i] = id
        @inbounds xs[i] = v
        nothing
    end
end

function Base.fill!(A::UnionVector, x)
    v = convert(eltype(A), x)
    p = v  # TODO: handle padding
    view_and_id(A, typeof(v)) do xs, id
        fill!(xs, p)
        fill!(A.typeid, id)
        nothing
    end
    return A
end


text:
#Test Boundary abstraction
using CavChaos

bnd0 = CircularBnd();
bnd1 = EllipticalBnd();
bnd2 = PolygonalBnd(5,0.5);
bnd3 = FourierBnd();

#Boundary test function
function testbndfunc(bnd::Boundary,theta::Array{Float64,1})
    println("========== Testing $(summary(bnd)) ==========");

    #Scalar functions
    tic();
    r0 = similar(theta);
    r1 = similar(theta); dtheta_r1 = similar(theta);
    r2 = similar(theta); dtheta_r2 = similar(theta); dtheta2_r2 = similar(theta);
    for i=1:length(theta)
        r0[i] = rfunc(bnd,theta[i]);
        r1[i],dtheta_r1[i] = rderiv(bnd,theta[i]);
        r2[i],dtheta_r2[i],dtheta2_r2[i] = rderiv2(bnd,theta[i]);
    end
    println("Scalar functions compute in $(toq())s.");

    #Array functions
    tic();
    r0A = rfunc(bnd,theta);
    r1A,dtheta_r1A = rderiv(bnd,theta);
    r2A,dtheta_r2A,dtheta2_r2A = rderiv2(bnd,theta);
    println("Array functions compute in $(toq())s.");

    #Test agreement
    print("Testing function agreement... ");
    if !(r0 == r1 == r2 == r0A == r1A == r2A)
        println("r agreement failed!");
        println("Returning [theta r0 r1 r2 r0A r1A r2A].");
        where = find(!(r0 .== r1 .== r2 .== r0A .== r1A .== r2A));
        return [theta[where] r0[where] r1[where] r2[where] r0A[where] r1A[where] r2A[where]];
    end
    if !(dtheta_r1 == dtheta_r2 == dtheta_r1A == dtheta_r2A)
        println("dtheta_r agreement failed!");
        println("Returning [theta dtheta_r1 dtheta_r2 dtheta_r1A dtheta_r2A].");
        where = find(!(dtheta_r1 .== dtheta_r2 .== dtheta_r1A .== dtheta_r2A));
        return [theta[where] dtheta_r1[where] dtheta_r2[where] dtheta_r1A[where] dtheta_r2A[where]];
    end
    if !(dtheta2_r2 == dtheta2_r2A)
        println("dtheta2_r agreement failed.");
        println("Returning [theta dtheta2_r2 dtheta2_r2A].");
        where = find(!(dtheta2_r2 .== dtheta2_r2A));
        return [theta[where] dtheta2_r2[where] dtheta2_r2A[where]];
    end
    println("passed.");

    #Estimate derivatives
    dtheta = theta[2];
    dtheta_r0 = (r0[2:end]-r0[1:end-1])/dtheta; #length-1
    dtheta2_r0 = (dtheta_r0[2:end]-dtheta_r0[1:end-1])/dtheta; #length-2
    dtheta2_r1 = (dtheta_r1[2:end]-dtheta_r1[1:end-1])/dtheta; #length-1
    dtheta_r1m = 0.5*(dtheta_r1[1:end-1]+dtheta_r1[2:end]);
    dtheta2_r2m = 0.5*(dtheta2_r2[1:end-1]+dtheta2_r2[2:end]);
    if typeof(bnd) == PolygonalBnd
        #exclude corners for PolygonalBnd
        angdiv = 2*pi/bnd.sides;
        sides = (abs(mod(theta+0.5*angdiv,angdiv)-0.5*angdiv) .> dtheta)
        error_0 = maximum(abs(dtheta_r0 - dtheta_r1m)[sides[1:end-1]]);
        error2_0 = maximum(abs(dtheta2_r0 - dtheta2_r2[2:end-1])[sides[2:end-1]]);
        error2_1 = maximum(abs(dtheta2_r1 - dtheta2_r2m)[sides[1:end-1]]);
    else
        error_0 = maximum(abs(dtheta_r0 - dtheta_r1m));
        error2_0 = maximum(abs(dtheta2_r0 - dtheta2_r2[2:end-1]));
        error2_1 = maximum(abs(dtheta2_r1 - dtheta2_r2m));
    end

    #Test accuracy
    tol = 1e-4;
    print("Testing accuracy of derivatives... ");
    if error_0 > tol
        println("dtheta_r1 error is $(string(error_0))!");
        println("Returning [theta dtheta_r0 dtheta_r1m].");
        where = find(abs(dtheta_r0 - dtheta_r1m) .> tol);
        return [theta[where] dtheta_r0[where] dtheta_r1m[where]];
    end
    if error2_0 > tol
        println("dtheta2_r2 error is $(tring(error2_0))!");
        println("Returning [theta dtheta2_r0 dtheta2_r2].");
        where = find(abs(dtheta2_r0 - dtheta2_r2[2:end-1]) .> tol);
        return [theta[2:end-1][where] dtheta2_r0[where] dtheta2_r2[2:end-1][where]];
    end
    if error2_1 > tol
        println("dtheta2_r2 error is $(string(error2_1))!");
        println("Returning [theta dtheta2_r1 dtheta2_r2m].");
        where = find(abs(dtheta2_r1 - dtheta2_r2m) .> tol);
        return [theta[where] dtheta2_r1[where] dtheta2_r2m[where]];
    end
    println("passed.");

    return "passed";
end

function testbnds()
    theta = linspace(0,2*pi,1000);

    println("############### BOUNDARY TEST ###############");

    testinfo0 = testbndfunc(bnd0,theta);
    testinfo1 = testbndfunc(bnd1,theta);
    testinfo2 = testbndfunc(bnd2,theta);
    testinfo3 = testbndfunc(bnd3,theta);

    print("############### ");
    if testinfo0 == testinfo1 == testinfo2 == testinfo3 == "passed"
        print("PASSED ALL");
    else
        print("FAILED");
    end
    println(" ###############\n");

    return (testinfo0,testinfo1,testinfo2,testinfo3);
end

testbnds()


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

using Markdown
using InteractiveUtils

# ╔═╡ 194baef2-0417-11ec-05ab-4527ef614024
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 38c9ced5-dcf8-4e03-ac07-7c435687861b
using FlowOverCircle, Plots

# ╔═╡ 50ce80a3-a1e8-4ba9-a032-dad315bcb432
md"
# Markov Neural Operator

JingYu Ning
"

# ╔═╡ 59769504-ebd5-4c6f-981f-d03826d8e34a
md"
This demo trains a Markov neural operator (MNO) introduced by [Zongyi Li *et al.*](https://arxiv.org/abs/2106.06898) with only one time step information. Then composed the operator to a Markov chain and inference the Navier-Stokes equations."

# ╔═╡ 823b3547-6723-43cf-85e6-cc6eb44efea1
md"
## Generate data
"

# ╔═╡ 5268feee-bda2-4612-9d4c-a1db424a11c7
begin
	n = 20
	data = FlowOverCircle.gen_data(LinRange(100, 100+n-1, n))
end;

# ╔═╡ 9b02b6a2-33c3-4ca6-bfba-0bd74b664830
begin
	anim = @animate for i in 1:size(data)[end]
		heatmap(data[1, 2:end-1, 2:end-1, i]', color=:coolwarm, clim=(-1.5, 1.5))
		scatter!(
			[size(data, 3)÷2-2.5], [size(data, 3)÷2-2.25],
			markersize=40, color=:black, legend=false, ticks=false
		)
		annotate!(5, 5, text("i=$i", :left))
	end
	gif(anim, fps=2)
end

# ╔═╡ 55058635-c7e9-4ee3-81c2-0153e84f4c8e
md"
## Inference

Use the first data generated above as the initial state, and apply the operator recurrently.
"

# ╔═╡ fbc287b8-f232-4350-9948-2091908e5a30
begin
	m = FlowOverCircle.get_model()

	states = Array{Float32}(undef, size(data))
	states[:, :, :, 1] .= view(data, :, :, :, 1)
	for i in 2:size(data)[end]
		states[:, :, :, i:i] .= m(view(states, :, :, :, i-1:i-1))
	end
end

# ╔═╡ a0b5e94c-a839-4cc0-a325-1a4ac39fafbc
begin
	anim_model = @animate for i in 1:size(states)[end]
		heatmap(states[1, 2:end-1, 2:end-1, i]', color=:coolwarm, clim=(-1.5, 1.5))
		scatter!(
			[size(data, 3)÷2-2.5], [size(data, 3)÷2-2.25],
			markersize=40, color=:black, legend=false, ticks=false
		)
		annotate!(5, 5, text("i=$i", :left))
	end
	gif(anim_model, fps=2)
end

# ╔═╡ Cell order:
# ╟─50ce80a3-a1e8-4ba9-a032-dad315bcb432
# ╟─59769504-ebd5-4c6f-981f-d03826d8e34a
# ╟─194baef2-0417-11ec-05ab-4527ef614024
# ╠═38c9ced5-dcf8-4e03-ac07-7c435687861b
# ╟─823b3547-6723-43cf-85e6-cc6eb44efea1
# ╠═5268feee-bda2-4612-9d4c-a1db424a11c7
# ╟─9b02b6a2-33c3-4ca6-bfba-0bd74b664830
# ╟─55058635-c7e9-4ee3-81c2-0153e84f4c8e
# ╠═fbc287b8-f232-4350-9948-2091908e5a30
# ╟─a0b5e94c-a839-4cc0-a325-1a4ac39fafbc


text:
# using LarSurf
using DICOM
using Distributed
if nprocs() == 1
    addprocs(3)
end
include("dicom_support.jl")
# Dataset can be downloaded from here:
# https://www.ircad.fr/research/3d-ircadb-01/
# this is path where data on all my computers are stored
dicomdir = homedir() * "/data/medical/orig/3Dircadb1.1/MASKS_DICOM/portalvein/"

datap = read_dicom_dir(dicomdir, "image_")

segmentation = datap["data3d"]
voxelsize_mm = datap["voxelsize_mm"]
using LarSurf
tm = time()
@everywhere using LarSurf
@everywhere using Distributed
tm = time() - tm
@info "everywhere time: $tm"

# tm = time()
# @everywhere using LarSurf
# @everywhere using Distributed
# tm = time() - tm
# @info "everywhere time second: $tm"
#
block_size = [64, 64, 64]
setup_time = @elapsed LarSurf.lsp_setup(block_size)
println("setup time: $setup_time")

tmd = @timed larmodel = LarSurf.lsp_get_surface(segmentation; voxelsize=voxelsize_mm)
val, tm, mem, gc = tmd
println("Total time: $tm")
V, FV = larmodel
FVtri = LarSurf.triangulate_quads(FV)

@info "Loading ViewerGL"
using ViewerGL
ViewerGL.VIEW([
    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

FVtri = LarSurf.triangulate_quads(FV)
# tm1 = @elapsed EV1 = LarSurf.Smoothing.get_EV_quads(FVtri)
# tm2 = @elapsed EV2 = LarSurf.Smoothing.get_EV_quads2(FVtri)
# @info "time get EV" tm1 tm2
# Plasm.view(val)
ViewerGL.VIEW([
    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,0,1,1))
    # ViewerGL.GLGrid(V,EV1,ViewerGL.Point4d(.9,0,.9,0.9))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

iterations = 3
Vs = LarSurf.Smoothing.smoothing_FV(V, FVtri, 0.6, iterations)

ViewerGL.VIEW([
    ViewerGL.GLGrid(Vs,FVtri,ViewerGL.Point4d(1,0,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])

iterations = 3
Vtaubin = LarSurf.Smoothing.smoothing_FV_taubin(V, FVtri, 0.4, -0.35, iterations)

ViewerGL.VIEW([
    ViewerGL.GLGrid(Vtaubin,FVtri,ViewerGL.Point4d(1,0,1,0.1))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])


text:
# ------------------- Ch. 9: Online Planning --------------------

# Reachable state space usually much smaller than than full state space.
#
# Algorithms to compute approximately optimal policies starting
#   from the current state.

export  RolloutLookahead,
        randstep,
        rollout,
        ForwardSearch,
        forward_search,
        BranchAndBound,
        branch_and_bound,
        SparseSampling,
        sparse_sampling,
        MonteCarloTreeSearch,
        simulate!,
        bonus,
        explore,
        HeuristicSearch,
        LabeledHeuristicSearch,
        expand,
        label!

# -------- Rollout with Lookahead ---------

"""
    struct RolloutLookahead
        𝒫 # problem
        π # rollout policy
        d # depth
    end
"""
struct RolloutLookahead
    𝒫 # problem
    π # rollout policy
    d # depth
end

"""
    randstep(𝒫::MDP, s, a) = 𝒫.TR(s, a)
"""
randstep(𝒫::MDP, s, a) = 𝒫.TR(s, a)

"""
    function rollout(𝒫, s, π, d)

Complexity: 𝒪(d)
"""
function rollout(𝒫, s, π, d)
    if d ≤ 0
        return 0.0
    end
    a = π(s)
    s′, r = randstep(𝒫, s, a)
    return r + 𝒫.γ*rollout(𝒫, s′, π, d-1)
end

"""
    function (π::RolloutLookahead)(s)

Policy function: similar to ValueFunctionPolicy but using a single call to rollout for value function.

Complexity: O(|𝒜|×|𝒮| × d)

Consider variation that uses average of m rollouts to calculate U(s) (vs just 1 with this version).
"""
function (π::RolloutLookahead)(s)
    U(s) = rollout(π.𝒫, s, π.π, π.d)
    return greedy(π.𝒫, U, s).a
end

# -------- Forward Search ---------
"""
    struct ForwardSearch
        𝒫 # problem
        d # depth
        U # value function at depth d
    end
"""
struct ForwardSearch
    𝒫 # problem
    d # depth
    U # value function at depth d
end

"""
    function forward_search(𝒫, s, d, U)

Determine optimal action (and its value) to take from state `s` by expanding all possible transitions up to a depth `d` using depth-first search. `U(s)` used to evaluate the terminal (depth 0) value function.

Complexity: 𝒪( (|𝒜|×|𝒮|)ᵈ )
"""
function forward_search(𝒫, s, d, U)
    if d ≤ 0
        return (a=nothing, u=U(s))
    end

    best = (a=nothing, u=-Inf)
    U′(s) = forward_search(𝒫, s, d-1, U).u

    for a in 𝒫.𝒜
        u = lookahead(𝒫, U′, s, a)
        if u > best.u
            best = (a=a, u=u)
        end
    end

    return best
end

(π::ForwardSearch)(s) = forward_search(π.𝒫, s, π.d, π.U).a

# --------- Branch and Bound ---------

"""
    struct BranchAndBound
        𝒫 # problem
        d # depth
        Ulo # lower bound on value function at depth d Qhi # upper bound on action value function
    end
"""
struct BranchAndBound
    𝒫 # problem
    d # depth
    Ulo # lower bound on value function at depth d
    Qhi # upper bound on action value function
end

"""
    function branch_and_bound(𝒫, s, d, Ulo, Qhi)
"""
function branch_and_bound(𝒫, s, d, Ulo, Qhi)
    if d ≤ 0
        return (a=nothing, u=Ulo(s))
    end

    U′(s) = branch_and_bound(𝒫, s, d-1, Ulo, Qhi).u
    best = (a=nothing, u=-Inf)
    for a in sort(𝒫.𝒜, by=a->Qhi(s,a), rev=true)
        if Qhi(s, a) < best.u
            return best # safe to prune
        end
        u = lookahead(𝒫, U′, s, a)
        if u > best.u
            best = (a=a, u=u)
        end
    end
    return best
end

(π::BranchAndBound)(s) = branch_and_bound(π.𝒫, s, π.d, π.Ulo, π.Qhi).a

# --------- Sparse Sampling ---------

"""
    struct SparseSampling
        𝒫 # problem
        d # depth
        m # number of samples
        U # value function at depth d
    end
"""
struct SparseSampling
    𝒫 # problem
    d # depth
    m # number of samples
    U # value function at depth d
end

"""
    function sparse_sampling(𝒫, s, d, m, U)

Complexity: 𝒪( (|𝒜| × m)ᵈ )
"""
function sparse_sampling(𝒫, s, d, m, U)
    if d ≤ 0
        return (a=nothing, u=U(s))
    end

    best = (a=nothing, u=-Inf)

    for a in 𝒫.𝒜
        u = 0.0
        for i in 1:m
            s′, r = randstep(𝒫, s, a)
            a′, u′ = sparse_sampling(𝒫, s′, d-1, m, U)
            u += (r + 𝒫.γ*u′) / m  # update mean
        end

        if u > best.u
            best = (a=a, u=u)
        end
    end

    return best
end

(π::SparseSampling)(s) = sparse_sampling(π.𝒫, s, π.d, π.m, π.U).a

# --------- Monte Carlo Tree Search ---------

"""
    struct MonteCarloTreeSearch
        𝒫 # problem
        N # visit counts
        Q # action value estimates
        d # depth
        m # number of simulations
        c # exploration constant
        π # rollout policy
    end
"""
struct MonteCarloTreeSearch
    𝒫 # problem
    N # visit counts
    Q # action value estimates
    d # depth
    m # number of simulations
    c # exploration constant
    π # rollout policy
end

"""
    function (π::MonteCarloTreeSearch)(s)

Estimate optimal next action by first simulating m random trajectories.  Each simulation improves estimated action-value function.  The returned action is the optimal action over the estimated action-value function from the given state `s`.

Complexity: 𝒪(|𝒜| × d × m)
"""
function (π::MonteCarloTreeSearch)(s)
    # simulate m times (m trajectories)
    for k in 1:π.m
        simulate!(π, s)
    end
    # Return action with highest action-value
    return _argmax(a->π.Q[(s,a)], π.𝒫.𝒜)
end

"""
    function simulate!(π::MonteCarloTreeSearch, s, d=π.d)

Execute exploration-bonus policy for 1 random trajectory, updating system  counters, and return sample trajectory's value.

Complexity: `𝒪(|𝒜| × d)`
"""
function simulate!(π::MonteCarloTreeSearch, s, d=π.d)
    if d ≤ 0
        return 0.0
    end

    𝒫, N, Q, c = π.𝒫, π.N, π.Q, π.c
    𝒜, TR, γ = 𝒫.𝒜, 𝒫.TR, 𝒫.γ

    # if never visited state s, init visit count and action-value estimate.
    if !haskey(N, (s, first(𝒜)))
        for a in 𝒜
            N[(s,a)] = 0
            Q[(s,a)] = 0.0
        end
        # first time to state s: return sample trajectory value as utility est.
        return rollout(𝒫, s, π.π, d)
    end

    # find next action to try
    a = explore(π, s) # 𝒪(|𝒜|)
    # simulate 1 step
    s′, r = TR(s,a)
    # recurse remaining stochastic trajectory
    q = r + γ*simulate!(π, s′, d-1)
    # update visit count
    N[(s,a)] += 1
    # update running average of action-value
    Q[(s,a)] += (q-Q[(s,a)])/N[(s,a)]

    return q
end

"""
    bonus(Nsa, Ns)

Monte Carlo exploration bonus term helper function.

Complexity: 𝒪(1)
"""
bonus(Nsa, Ns) = Nsa == 0 ? Inf : sqrt(log(Ns)/Nsa)

"""
    function explore(π::MonteCarloTreeSearch, s)

Use monte carlo exploration heuristic to find next action to simulate.  Balances need to explore state-action space with
Complexity: 𝒪(|𝒜|)
"""
function explore(π::MonteCarloTreeSearch, s)

    𝒜, N, Q, c = π.𝒫.𝒜, π.N, π.Q, π.c

    Ns = sum(N[(s,a)] for a in 𝒜)

    return _argmax(a->Q[(s,a)] + c*bonus(N[(s,a)], Ns), 𝒜)
end

# --------- Heuristic Search ---------

# Use m simulations of a greedy policy with respect to value function U from state s.  U is initialized to an upperbound Ū (referred to as a heuristic).  Updates U with each lookahead step during simulation.  After simulations and value function estimate U improvement, return greedy action.
"""
    struct HeuristicSearch
        𝒫 # problem
        Uhi # upper bound on value function
        d # depth
        m # number of simulations
    end
"""
struct HeuristicSearch
    𝒫 # problem
    Uhi # upper bound on value function
    d # depth
    m # number of simulations
end

"""
    function simulate!(π::HeuristicSearch, U, s)

Simulate depth-d trajectory following greedy policy and updating U along the way.

Complexity: 𝒪(d × |𝒜|×|𝒮|)
"""
function simulate!(π::HeuristicSearch, U, s)
    𝒫, d = π.𝒫, π.d
    for d in 1:d                # 𝒪( d × ... (|𝒜|×|𝒮|))
        a, u = greedy(𝒫, U, s) # 𝒪(|𝒜|×|𝒮|)
        U[s] = u
        s = rand(𝒫.T(s, a))
    end
end

"""
    function (π::HeuristicSearch)(s)

Guranteed to converge to optimal value function iff `Uhi` is indeed an upperbound on `U`.

Complexity: 𝒪(m × d × |𝒮|×|𝒜|)
"""
function (π::HeuristicSearch)(s)
    # initialize U with upper bound
    U = [π.Uhi(s) for s in π.𝒫.𝒮]
    # run m random trajectory simulations to improve estimate on U.
    for i in 1:m
        simulate!(π, U, s)
    end
    # Return greedy action using improved U.
    return greedy(π.𝒫, U, s).a
end

# --------- Labled Heuristic Search ---------

"""
    struct LabeledHeuristicSearch
        𝒫 # problem
        Uhi # upper bound on value function
        d # depth
        δ # gap threshold
    end
"""
struct LabeledHeuristicSearch
    𝒫 # problem
    Uhi # upper bound on value function
    d # depth
    δ # gap threshold
end

"""
    function (π::LabeledHeuristicSearch)(s)
"""
function (π::LabeledHeuristicSearch)(s)
    U, solved = [π.Uhi(s) for s in 𝒫.𝒮], Set()
    while s ∉ solved
        simulate!(π, U, solved, s)
    end
    return greedy(π.𝒫, U, s).a
end

"""
    function simulate!(π::LabeledHeuristicSearch, U, solved, s)
"""
function simulate!(π::LabeledHeuristicSearch, U, solved, s)
    visited = []
    for d in 1:π.d
        if s ∈ solved
            break
        end

        push!(visited, s)
        a, u = greedy(π.𝒫, U, s)
        U[s] = u
        s = rand(π.𝒫.T(s, a))
    end

    while !isempty(visited)
        if label!(π, U, solved, pop!(visited))
            break
        end
    end
end

"""
    function expand(π::LabeledHeuristicSearch, U, solved, s)
"""
function expand(π::LabeledHeuristicSearch, U, solved, s)
    𝒫, δ = π.𝒫, π.δ
    𝒮, 𝒜, T = 𝒫.𝒮, 𝒫.𝒜, 𝒫.T

    found, toexpand, envelope = false, Set(s), []
    while !isempty(toexpand)
        s = pop!(toexpand)
        push!(envelope, s)
        a, u = greedy(𝒫, U, s)
        if abs(U[s] - u) > δ
            found = true
        else
            for s′ in 𝒮
                if T(s,a,s′) > 0 && s′ ∉ (solved ∪ envelope)
                    push!(toexpand, s′)
                end
            end
        end
    end

    return (found, envelope)
end

"""
    function label!(π::LabeledHeuristicSearch, U, solved, s)
"""
function label!(π::LabeledHeuristicSearch, U, solved, s)
    if s ∈ solved
        return false
    end
    found, envelope = expand(π, U, solved, s)
    if found
        for s ∈ reverse(envelope)
            U[s] = greedy(π.𝒫, U, s).u
        end
    else
        union!(solved, envelope)
    end

    return found
end


text:
function verifyQuadRule(qr::TriTetQuadRule{N_zeta, N_Point}) where {N_zeta, N_Point}
    (sum(qr.w) ≈ 1.0) || (return false)
    for k = 1:N_Point
        (sum(qr.zeta[k]) ≈ 1.0) || (return false)
    end
    c = zeros(SVector{N_zeta,Float64})
    for k = 1:N_Point
        c += qr.w[k] * qr.zeta[k]
    end
    @test all(c .≈ (1 / N_zeta))
    return true
end

@testset "quadrature" begin
    # Triangle
    for k = 1:5
        @test verifyQuadRule(getTriQuadRule(k))
    end
    qr_2 = getTriQuadRule(2)
    @test all(sort(qr_2.zeta[1]) .== sort(qr_2.zeta[2]) .== sort(qr_2.zeta[3]))

    # Tetrahedron
    for k = 1:5
        @test verifyQuadRule(getTetQuadRule(k))
    end
    @test_throws ErrorException getTriQuadRule((6))
    @test_throws ErrorException getTetQuadRule((6))
end


text:
const MBVec = [
"ADASYN",
"AllKNN",
"BorderlineSMOTE",
"CondensedNearestNeighbour",
"EditedNearestNeighbours",
"InstanceHardnessThreshold",
"RandomOverSampler",
"RepeatedEditedNearestNeighbours",
"SMOTE",
"SMOTEENN",
"SVMSMOTE"
]


text:
# runtests.jl

using Test: @test, @testset, detect_ambiguities
using ImagePhantoms

include("helper.jl")

include("core.jl")
include("ellipse.jl")
include("gauss2.jl")
include("rect.jl")
include("triangle.jl")
include("shepplogan.jl")
include("disk-phantom.jl")
include("focus-chart.jl")

@testset "mri-sense" begin
    include("mri-sense.jl")
end

@testset "ImagePhantoms" begin
    @test isempty(detect_ambiguities(ImagePhantoms))
end


text:
"""
    air.(broadcasting_expression) :: Broadcasted

Broadcast without materialization.

The idea is taken from @dawbarton's `_lazy` function:
https://discourse.julialang.org/t/19641/20.
"""
air(::Any) =
    error("Function `air` must be called with a dot as in: ",
          "`air.(YOUR .+ BROADCASTING .* EXPRESSION ./ HERE)`")

struct Aired{T}
    value::T
end
@inline Broadcast.broadcasted(::typeof(air), x) = Aired(x)
@inline Broadcast.materialize(x::Aired) = x.value

"""
    (@~ broadcasting_expression) :: Broadcasted
"""
macro ~(ex)
    esc(:($air.($ex)))
end


text:
using Dates
include("Network.jl")

"""
    flow(;apikey, bbox, s=1.0, m=1/60)

Returns real-time vehicle speed in every s seconds for m minutes in a bounding box using HERE Traffic API.

# Arguments
apikey::String   : API key to access HERE
bbox::String     : Bounding box defined as `"lat₁,long₁;lat₂,long₂"`
s::Float64=1.0   : gap between every API call
m::Float64=1/60  : time-period of analysis

For more information visit: https://developer.here.com/documentation/traffic/dev_guide/topics/guide.html

Example request: https://developer.here.com/documentation/traffic/dev_guide/topics_v6.1/example-flow-bounding-box.html
"""
function flow(;apikey, bbox, s=1.0, m=1/60)
    # Create request url
    baseurl = "https://traffic.ls.hereapi.com"
    path = "traffic/6.2"
    resource = "flow.xml"
    url = "$baseurl/$path/$resource?apiKey=$apikey&bbox=$bbox&responseattributes=shape"

    # Fetch network
    df = network(apikey=apikey, bbox=bbox)

    # Fetch real-time speeds
    t = 1
    while t ≤ m * 60
        Tₒ = Dates.format(now(), "HH:MM:SS:sss")
        df[!, "rts @$Tₒ"] = zeros(nrow(df))
        r = HTTP.request(:GET, url)
        f = parsexml(String(r.body))
        rootnode = root(f)
        k = 1
        for node in iterate(rootnode)
            if nodename(node) == "FI"
                childnode = lastelement(node)
                v = parse(Float64, childnode["SU"])
                df[k, "rts @$Tₒ"] = v
                k += 1
            end
        end
        t += s
        sleep(s)
    end

    return df
end


text:
#! /Users/malcolm/bin/julia
#
# Check on the number of arguments, print usage unless 1 or 2.
nargs = length(ARGS);
if nargs == 0 || nargs > 2
    println("usage: etl.jl infile [outfile]");
    exit();
end

# Assign first argument to input file
# If second argument this is the output file otherwise STDOUT
infile = ARGS[1];
if nargs == 2
   outfile = ARGS[2];
   try
      outf = open(outfile,"w");
   catch
      error("Can't create output file: ", outfile);
   end
else
   outf = STDOUT;
end

# One liner to double up single quotes
escticks(s) = replace(s,"'","''");

# Read all file into a matrix, first dimensions is number of lines
qq = readdlm("quotes.tsv",'\t');
n = size(qq)[1];

# Going to store all categories in a dictionary
j = 0;
cats = Dict(String,Int64);

# Main loop to load up the quotes table
for i = 1:n
  cat = qq[i,1];
  if haskey(cats,cat)
    jd = cats[cat];
  else
    j = j + 1; jd = j;
    cats[cat] = jd;
  end
  sql = "insert into quotes values($i,$jd,";
  if (length(qq[i,2]) > 0)
    sql *= string("'", escticks(qq[i,2]), "',");
  else
    sql *= string("null,");
  end
  sql *= string("'", escticks(qq[i,3]), "');");
  wrlte(outf,"$sql\n");
end

# Now dump the categories
for cat = keys(cats)
  jd = cats[cat];
  write(outf,"insert into categories($jd,'$cat');\n");
end

close(outf);	# Will have no effect if outf = STDOUT

max_stars_repo_path: Module 3/Chapter08/etl.jl

max_issues_repo_path: Module 3/Chapter08/etl.jl

max_forks_repo_path: Module 3/Chapter08/etl.jl


text:
# UDT decomposition (QR decomposition with U = Q and DT = R)
struct UDT{E,Er<:Real,M<:AbstractArray{E}} <: Factorization{E}
    U::M
    D::Vector{Er}
    T::M
    function UDT{E,Er,M}(U, D, T) where {E,Er,M<:AbstractArray{E}}
        new{E,Er,M}(U, D, T)
    end
end

UDT(U::AbstractArray{E}, D::Vector{Er}, T::AbstractArray{E}) where {E,Er<:Real} =
    UDT{E,Er,typeof(U)}(U, D, T)
function UDT{E}(U::AbstractArray, D::AbstractVector{Er}, T::AbstractArray) where {E,Er<:Real}
    UDT(convert(AbstractArray{E}, U),
        convert(Vector{Er}, D),
        convert(AbstractArray{E}, T))
end

# Conversion
Base.AbstractMatrix(F::UDT) = (F.U * Diagonal(F.D)) * F.T
Base.AbstractArray(F::UDT) = AbstractMatrix(F)
Base.Matrix(F::UDT) = Array(AbstractArray(F))
Base.Array(F::UDT) = Matrix(F)

function Matrix!(res::AbstractMatrix, F::UDT)
    tmp = F.U * Diagonal(F.D)
    mul!(res, tmp, F.T)
end

# iteration for destructuring into components
Base.iterate(S::UDT) = (S.U, Val(:D))
Base.iterate(S::UDT, ::Val{:D}) = (S.D, Val(:T))
Base.iterate(S::UDT, ::Val{:T}) = (S.T, Val(:done))
Base.iterate(S::UDT, ::Val{:done}) = nothing

Base.size(A::UDT, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.T, dim)
Base.size(A::UDT) = (size(A, 1), size(A, 2))

Base.similar(A::UDT) = UDT(similar(A.U), similar(A.D), similar(A.T))

# decomposition functions
"""
Compute the UDT decomposition of `A` and return an `UDT` object.

`U`, `D`, and `T`, can be obtained from the factorization `F`
with `F.U`, `F.D`, and `F.T` such that `A = U * Diagonal(D) * T`.

Iterating the decomposition produces the components `U`, `D`, and `V`.

Note that `T` is upper triangular only up to permutations of columns of `T`.
"""
function udt(A::AbstractMatrix{C}) where {C<:Number}
  F = qr(A, Val(true))
  _qr_to_udt(A, F)
end

"""
`udv!` is the same as `svd`, but saves space by overwriting the input `A`,
instead of creating a copy.
"""
function udt!(A::AbstractMatrix{C}) where {C<:Number}
  F = qr!(A, Val(true))
  _qr_to_udt(A, F)
end

@inline function _qr_to_udt(A::AbstractMatrix{C}, F::QRPivoted) where {C<:Number}
    n = size(A, 1)
    D = Vector{real(float(C))}(undef, n)
    R = F.R # F.R has regular matrix type
    @views F.p[F.p] = 1:n

    @inbounds for i in 1:n
    D[i] = abs(real(R[i,i]))
    end
    lmul!(Diagonal(1 ./ D), R)
    UDT(Matrix(F.Q), D, R[:, F.p])
end

function udt(x::Number)
    UDT(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
end
function udt(x::Integer)
    udt(float(x))
end

# operations
"""
    inv(F::UDT) -> AbstractMatrix

Computes the inverse matrix of the `UDT` decomposition of a matrix.
"""
function Base.inv(F::UDT)
    inv!(similar(F.U), F)
end

"""
    inv!(res, F::UDT) -> res

Same as `inv` but writes result into preallocated `res`.
"""
function inv!(res::M, F::UDT{E, Er, M}) where {E,Er,M}
    tmp = similar(F.U)
    ldiv!(tmp, lu(F.T), Diagonal(1 ./ F.D))
    mul!(res, tmp, F.U')
    return res
end

"""
    fact_mult(A::UDT, B::UDT) -> UDT

Stabilized multiplication of two `UDT` decompositions.
Returns a `UDT` factorization object.
"""
function fact_mult(A::UDT, B::UDT)
    mat = A.T * B.U
    lmul!(Diagonal(A.D), mat)
    rmul!(mat, Diagonal(B.D))
    F = udt!(mat)
    UDT(A.U * F.U, F.D, F.T * B.T)
end

##############################################################
#
#                   QR / UDT
#
##############################################################
"""
    udt_inv_one_plus(F::UDT) -> UDT

Stabilized calculation of [1 + UDT]^(-1). Returns and
`UDT` factorization object.

Optional preallocations via keyword arguments:

  * `u = similar(F.U)`
  * `t = similar(F.T)`
"""
function udt_inv_one_plus(F::UDT; u = similar(F.U), t = similar(F.T))
  U, D, T = F

  m = U' / T
  m[diagind(m)] .+= D

  utmp, d, ttmp = udt!(m)
  mul!(u, U, utmp)
  mul!(t, ttmp, T)

  UDT(inv(t), 1 ./ d, copy(u'))
end

"""
  inv_one_plus!(res, F::UDT) -> res

Same as `inv_one_plus` but stores the result in preallocated `res`.
"""
function inv_one_plus!(res, F::UDT;
                       u = similar(F.U),
                       d = similar(F.D),
                       t = similar(F.T))
  U, D, T = F

  m = U' / T
  m[diagind(m)] .+= D

  utmp, d, ttmp = udt!(m)
  mul!(u, U, utmp)
  mul!(t, ttmp, T)

  ldiv!(m, lu!(t), Diagonal(1 ./ d))

  mul!(res, m, u')
  res
end

"""
    inv_one_plus(F::UDT) -> AbstractMatrix

Stabilized calculation of `[1 + UDT]^(-1)`:

  * Use one intermediate UDT decomposition.

Faster but potentially less accurate than `inv_one_plus_loh`.

See `udt_inv_one_plus` for preallocation options.
"""
function inv_one_plus(F::UDT; kwargs...)
  res = similar(F.U)
  inv_one_plus!(res, F; kwargs...)
  return res
end

"""
    udt_inv_one_plus(A::UDT, Bdagger::UDT) -> UDT

Stabilized calculation of [1 + UlDlTl(UrDrTr)^†]^(-1). Returns and
`UDT` factorization object.

Optional preallocations via keyword arguments:

  * `tmp = similar(A.U)`
  * `tmp2 = similar(A.U)`
  * `tmp3 = similar(A.U)`
"""
function udt_inv_one_plus(A::UDT, Bdagger::UDT;
                          tmp = similar(A.U),
                          tmp2 = similar(A.U),
                          tmp3 = similar(A.U),
                          internaluse = false)
  Ul,Dl,Tl = A
  Ur,Dr,Tr = Bdagger

  mul!(tmp, Tl, adjoint(Tr))
  rmul!(tmp, Diagonal(Dr))
  lmul!(Diagonal(Dl), tmp)
  U1, D1, T1 = udt!(tmp)

  mul!(tmp3, Ul, U1)
  mul!(tmp2, T1, adjoint(Ur))
  mul!(tmp, adjoint(tmp3), inv(tmp2))

  tmp .+= Diagonal(D1)

  u, d, t = udt!(tmp)
  mul!(tmp, t, tmp2)
  mul!(tmp2, tmp3, u)

  if internaluse
    UDT(inv(tmp), 1 ./ d, tmp2')
  else
    UDT(inv(tmp), 1 ./ d, copy(tmp2'))
  end
end

"""
  inv_one_plus!(res, A::UDT, Bdagger::UDT) -> res

Stabilized calculation of [1 + UlDlTl(UrDrTr)^†]^(-1).
Writes the result into `res`.

See `udt_inv_one_plus` for preallocation options.
"""
function inv_one_plus!(res, A::UDT, Bdagger::UDT; kwargs...)
  F = udt_inv_one_plus(A, Bdagger; internaluse = true, kwargs...)
  rmul!(F.U, Diagonal(F.D))
  mul!(res, F.U, F.T)
  res
end

"""
  inv_one_plus(A::UDT, Bdagger::UDT) -> AbstractMatrix

Stabilized calculation of [1 + UlDlTl(UrDrTr)^†]^(-1).

See `udt_inv_one_plus` for preallocation options.
"""
function inv_one_plus(A::UDT, Bdagger::UDT; kwargs...)
  res = similar(A.U)
  inv_one_plus!(res, A, Bdagger; kwargs...)
  return res
end

"""
    udt_inv_sum(A::UDT, B::UDT) -> UDT

Stabilized calculation of [UaDaTa + UbDbTb]^(-1):

  * Use one intermediate UDT decompositions.

Optional preallocations via keyword arguments:

  * `m2 = similar(A.U)`
"""
function udt_inv_sum(A::UDT, B::UDT; m2 = similar(A.U))
  Ua,Da,Ta = A
  Ub,Db,Tb = B

  m1 = Ta / Tb
  lmul!(Diagonal(Da), m1)

  mul!(m2, Ua', Ub)
  rmul!(m2, Diagonal(Db))

  u,d,t = udt!(m1 + m2)

  mul!(m1, Ua, u)
  mul!(m2, t, Tb)

  UDT(inv(m2), 1 ./ d, m1')
end

"""
    inv_sum!(res, A::UDT, B::UDT) -> res

Stabilized calculation of [UaDaTa + UbDbTb]^(-1):

  * Use one intermediate UDT decompositions.

See `udt_inv_sum` for preallocation options.
"""
function inv_sum!(res, A::UDT, B::UDT; kwargs...)
  F = udt_inv_sum(A, B; kwargs...)
  rmul!(F.U, Diagonal(F.D))
  mul!(res, F.U, F.T)
  res
end

"""
    inv_sum(A::UDT, B::UDT) -> AbstractMatrix

Stabilized calculation of [UaDaTa + UbDbTb]^(-1):

  * Use one intermediate UDT decompositions.

See `udt_inv_sum` for preallocation options.
"""
function inv_sum(A::UDT, B::UDT; kwargs...)
  res = similar(A.U)
  inv_sum!(res, A, B; kwargs...)
  res
end

########################################
#
#         Loh et. al. schemes
#
########################################

"""
    udt_inv_one_plus_loh(F::UDT) -> UDT

Stabilized calculation of [1 + UDT]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Options for preallocation via keyword arguments:

  * `l = similar(F.U)`
  * `r = similar(F.U)`
  * `Dp = similar(F.D)`
  * `Dm = similar(F.D)`
"""
function udt_inv_one_plus_loh(F::UDT; r = similar(F.U),
                                      l = similar(F.U),
                                      Dp = similar(F.D),
                                      Dm = similar(F.D),
                                      internaluse = false)
  U,D,T = F

  Dp .= max.(D,1.)
  Dm .= min.(D,1.)

  Dp .\= 1
  Dpinv = Dp # renaming

  ldiv!(l, lu(T), Diagonal(Dpinv)) # Don't use lu!(T) because T is input

  mul!(r, U, Diagonal(Dm))
  r .+= l

  u, d, t = udt!(r)

  ldiv!(r, lu!(t), Diagonal(1 ./ d))
  mul!(l, r, u')

  lmul!(Diagonal(Dpinv), l)
  u, d, t = udt!(l)

  ldiv!(l, lu(T), u)

  if internaluse
    UDT(l, d, t)
  else
    UDT(copy(l), d, t)
  end
end

"""
  inv_one_plus_loh!(res, F::UDT) -> res

Stabilized calculation of [1 + UDT]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Writes the result into `res`.

See `udt_inv_one_plus_loh` for preallocation options.
"""
function inv_one_plus_loh!(res, F::UDT; kwargs...)
  X = udt_inv_one_plus_loh(F; internaluse = true, kwargs...)
  rmul!(X.U, Diagonal(X.D))
  mul!(res, X.U, X.T)
  res
end

"""
  inv_one_plus_loh(F::UDT) -> AbstractMatrix

Stabilized calculation of [1 + UDT]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

See `udt_inv_one_plus_loh` for preallocation options.
"""
function inv_one_plus_loh(F::UDT; kwargs...)
  res = similar(F.U)
  inv_one_plus_loh!(res, F; kwargs...)
  res
end

"""
    udt_inv_sum_loh(A::UDT, B::UDT) -> UDT

Stabilized calculation of [UaDaTa + UbDbTb]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Options for preallocations via keyword arguments:

  * `mat2 = similar(A.U)`
  * `Dap = similar(A.D)`
  * `Dam = similar(A.D)`
  * `Dbp = similar(A.D)`
  * `Dbm = similar(A.D)`

"""
function udt_inv_sum_loh(A::UDT, B::UDT; mat2 = similar(A.U),
                                         Dap = similar(A.D),
                                         Dam = similar(A.D),
                                         Dbp = similar(A.D),
                                         Dbm = similar(A.D),
                                         internaluse = false)
    Ua, Da, Ta = A
    Ub, Db, Tb = B

    d = length(Da)

    # separating scales larger and smaller than unity
    Dap .= max.(Da,1.)
    Dam .= min.(Da,1.)
    Dbp .= max.(Db,1.)
    Dbm .= min.(Db,1.)

    # mat1 = Dam * Vda * Vdb' / Dbp
    mat1 = Ta / Tb
    @inbounds for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] * Dam[j]/Dbp[k]
    end

    # mat2 = 1/(Dap) * Ua' * Ub * Dbm
    mul!(mat2, Ua', Ub)
    @inbounds for j in 1:d, k in 1:d
        mat2[j,k]=mat2[j,k] * Dbm[k]/Dap[j]
    end

    # mat1 = mat1 + mat2
    mat1 .+= mat2

    # decompose mat1: U, D, T
    U, D, T = udt!(mat1)

    # invert and combine inner part: mat1 = (U D T)^(-1)
    lmul!(Diagonal(D), T)
    ldiv!(mat1, lu!(T), U') # mat1 = T \ (U')

    # mat1 = 1/Dbp * mat1 /Dap
    @inbounds for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] / Dbp[j] / Dap[k]
    end

    #mat1 = U D T
    U, D, T = udt!(mat1)

    # U = Tb^(-1) * U , T = T * Ua'
    ldiv!(mat1, lu(Tb), U) # mat1 = Tb \ U
    mul!(mat2, T, Ua')

    if internaluse
      UDT(mat1, D, mat2)
    else
      UDT(mat1, D, copy(mat2))
    end
end

"""
    inv_sum_loh!(res, A::UDT, B::UDT) -> res

Stabilized calculation of [UaDaTa + UbDbTb]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Writes the result into `res`.

See `udt_inv_sum_loh` for preallocation options.
"""
function inv_sum_loh!(res, A::UDT, B::UDT; kwargs...)
  F = udt_inv_sum_loh(A, B; internaluse = true, kwargs...)
  rmul!(F.U, Diagonal(F.D))
  mul!(res, F.U, F.T)
  res
end

"""
    inv_sum_loh(A::UDT, B::UDT) -> AbstractMatrix

Stabilized calculation of [UaDaTa + UbDbTb]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

See `udt_inv_sum_loh` for preallocation options.
"""
function inv_sum_loh(A::UDT, B::UDT; kwargs...)
  res = similar(A.U)
  inv_sum_loh!(res, A, B; kwargs...)
  res
end


text:
category_trait(::Type{<:Ring}) = CommutativeRingTrait
category_trait(::Type{<:Number}) = IntegralDomainTrait

# promotions and conversions
function promote_rule(::Type{T}, ::Type{S}) where {T<:Ring,S<:Ring}
    dts = depth(T) - depth(S)
    if dts < 0
        promote_rule(S, T)
    elseif dts > 0
        B = basetype(T)
        if B == S
            T
        else
            promote_rule(B, S) == B ? T : _promote_rule(T, S)
        end
    else
        _promote_rule(T, S)
    end
end

promote_rule(::Type{R}, ::Type{S}) where {R<:Ring,S<:Rational}= _promote_rule(R, S)
promote_rule(::Type{S}, ::Type{R}) where {R<:Ring,S<:Rational} = _promote_rule(R, S)
_promote_rule(::Type{R}, ::Type{S}) where {R<:Ring,S<:Rational} = promote_rule(R, promote_type(basetype(R), S))

depth(::Type{<:Number}) = 0

for op in (:+, :-, :*, :/, :(==), :divrem, :div, :rem, :gcd, :gcdx, :pgcd, :pgcdx, :isless)
    @eval begin
        ($op)(a::Ring, b::Ring) = ($op)(promote(a, b)...)
        ($op)(a::Ring, b::Union{Integer,Rational}) = ($op)(promote(a, b)...)
        ($op)(a::Union{Integer,Rational}, b::Ring) = ($op)(promote(a, b)...)
    end
end

# generic operations
basetype(::T) where T<:Ring = basetype(T)
basetype(::Type{T}) where T = T

convert(::Type{T}, a) where T = T(a)
function convert(::Type{T}, a::S) where {T<:Ring,S<:Ring}
    if !(S <: basetype(T)) && depth(T) > depth(S)
            b = convert(basetype(T), a)
            convert(T, b)
    else
        T(a)
    end
end
(G::Type{<:Ring})(a) = G !== basetype(G) ? G(convert(basetype(G), a)) : throw(MethodError(G, a))
@generated function basetypes(a)
    _basetypes(::Type{a}) where a = begin b = basetype(a); a == b ? [a] : [a; _basetypes(b)] end
    bt = tuple(_basetypes(a.parameters[1])...)
    :($bt)
end

function /(a::T, b::T) where T<:Ring
    if b isa Union{FractionRing,QuotientRing}
        a * inv(b)
    else
        d, r = divrem(a, b)
        iszero(r) || throw(DomainError((a, b), "not dividable a/b."))
        T(d)
    end
end
\(a::Ring, b::Ring) = b / a
^(a::Ring, n::Integer) = Base.power_by_squaring(a, n)
zero(x::Ring) = zero(typeof(x))
one(x::Ring) = one(typeof(x))
inv(a::Ring) = isunit(a) ? 1 / a : throw(DomainError(a, "cannot divide by non-unit."))
# enable generic matrix factorization
abs(a::Ring) = isunit(a) ? 1 : 0
value(a::Ring) = a

import Base: literal_pow
@inline literal_pow(::typeof(^), x::Ring, ::Val{0}) = one(x)
@inline literal_pow(::typeof(^), x::Ring, ::Val{1}) = x
@inline literal_pow(::typeof(^), x::Ring, ::Val{2}) = x * x
@inline literal_pow(::typeof(^), x::Ring, ::Val{3}) = x * x * x

numerator(a::Ring) = a
denominator(::R) where R<:Ring = one(R)

"""
    order(::Type)

Returns the number of elements of (finite) ring `Z` or `0` if `|Z| == inf`.
"""
order(::Type) = 0
order(::Type{Z}) where Z<:ZZmod = uptype(modulus(Z))

"""
    dimension(::Type)

Returns the number of dimensions of vector-space like rings (like quotient rings over polynomials).
In al other cases return `1`
"""
dimension(R::Type) = 1

"""
    characteristic(::Type)

Returns characteristic `c` of ring `R`.
`c` is the smallest positive integer with `c * one(R) == 0`, or `0` if `c * one(R) != 0` for all positive integers `c`.
"""
characteristic(::Type) = 0

"""
    intpower(a, b)

Calculate `a ^ b` in appropriate result type.
"""
intpower(a::Integer, b::Integer) = uptype(a, mintype_for(a, b, false)) ^ b

"""
    uptype(a, [T::Type])

promote `a` to get at least type `promote_type(typeof(a), T)`.
"""
uptype(a::T, S::Type=Int) where T = promote_type(S, T)(a)

"""
    order(z::Ring)

Returns the order of the cyclic multiplicative group generated by z,
or `0` if `z^k != 1` for all integers `k`.
"""
function order(x::Ring)
    iszero(x) && return 0
    isone(x) && return 1
    n = mult_order(typeof(x))
    n <= 1 && return n
    r = trailing_zeros(n)
    f = factor(n >> r)
    for k in factors(f)
        y = x^k
        isone(y) && return k
        j = k
        for s = 1:r
            y *= y
            j += j
            isone(y) && return j
        end
    end
    0
end

"""
    mult_order(R::Type{<:Ring})

order of multiplicative subgroup of `R`.
"""
mult_order(R::Type) = 0
mult_order(R::Type{<:ZZmod}) = totient(modulus(R))
mult_order(R::Type{<:GaloisField}) = order(R) - 1
function mult_order(R::Type{<:Quotient{T,T}} where T<:UnivariatePolynomial)
    isirreducible(modulus(R)) ? order(R) - 1 : 0
end

"""
    characteristic(Z::Type{<:Ring})

Returns the characteristic `c` of ring `Z`. `c` is the smallest positive integer
with `c * one(Z) == 0`, or `0` if `c * one(Z) != 0` for all positive integers `c`.
"""
characteristic(::Type{Z}) where Z<:ZZmod = order(Z)
characteristic(::Type{T}) where {Z,T<:Quotient{Z}} = characteristic(Z)
characteristic(::Type{T}) where {Z,T<:Frac{Z}} = characteristic(Z)
characteristic(::Type{T}) where {Z,T<:Polynomial{Z}} = characteristic(Z)

"""
    isfield(::Type{<:Ring})

Return iff type is a field (all elements except `zero` are invertible).
"""
isfield(R::Type{<:Ring}) = category_trait(R) <: FieldTrait

"""
    deg(r::Union{Ring,Number})

Return the degree of ring element or number `r`.

For zero elements, that is `-1`, otherwise `0` for non-polynomials,
the ordinary degree for univariate polynomials and the
maximum sum of exponents for multivariate polynomials.
"""
deg(x::Union{Ring,Number}) = iszero(x) ? -1 : 0 # fallback

function divrem(a::T, b::T) where T<:Ring
    if isunit(b)
        (a * inv(b), zero(b))
    else
        throw(MethodError(divrem, (a, b)))
    end
end
div(a::T, b::T) where T<:Ring = divrem(a, b)[1]
rem(a::T, b::T) where T<:Ring = divrem(a, b)[2]
"""
    isdiv(a, b)

Return if ring element `a` is dividable by `b` without remainder.
"""
isdiv(a::T, b::T) where T <: Ring = iszero(rem(a, b))

"""
    iscoprime(a, b)

Return if there is a non-unit common divisor of `a` and `b`.
"""
iscoprime(a::T, b::T) where T<:Ring = isunit(a) || isunit(b)

divrem(a::T, b::T) where T<:QuotientRing = (a / b, zero(a))
"""
    LC(r::Ring)

Return the leading coefficient of `r`. For non-polynomials return `r` itself.
"""
LC(x::Ring) = x

"""
    lcunit(r::Ring)

Return `one(r)` if it's a unit, otherwise return a unit element `s` of the same ring or of
an object, which may be promoted to this ring, so `r / s` has a simplified form.
Example, for a polynomial over a field, the leading coefficient.
"""
lcunit(x::Ring) = one(x)

modulus(::T) where T<:Ring = modulus(T)
copy(p::QuotientRing) = typeof(p)(p.val)
# make Ring elements behave like scalars with broadcasting
Base.broadcastable(x::Ring) = Ref(x)

# apply homomorphism
(h::Hom{F,R,S})(a::R) where {F,R,S} = h.f(a)::S

divrem2(a::T, b::T) where T = divrem(a, b)
rem2(a::T, b::T) where T = divrem2(a, b)[2]

function divrem2(a::T, b::T) where T<:FractionRing
    c = div(a, b)
    d, r = divrem2(c.num, c.den)
    d, a - d * b
end

# generic Euclid's algorithm
function gcd(a::T, b::T) where T<:Ring
    iszero(b) && return a
    a, b = b, rem2(a, b)
    while !iszero(b)
        d, c = divrem2(a, b)
        if iszero(d)
            a = one(a)
            break
        end
        a, b = b, c
        issimpler(b, a) || throw(DomainError((a,b), "b is not simpler than a"))
    end
    u = lcunit(a)
    isone(a) ? a : a * inv(u)
end

# extension to array
function gcd(aa::Union{AbstractVector{T},NTuple{N,T}}) where {N,T<:Ring}
    n = length(aa)
    g = zero(T)
    n == 0 && return g
    for i = 1:n
        isone(g) && break
        g = gcd(aa[i], g)
    end
    g
end
gcd(a::T...) where T<:Ring = gcd(a)

pgcd(a::R, b::R) where R<:Ring = gcd(a, b)

# generic extended Euclid's algorithm
function gcdx(a::T, b::T) where T<:Ring
    s0, s1 = one(T), zero(T)
    t0, t1 = s1, s0
    # invariant: a * s0 + b * t0 == gcd(a, b)
    while !iszero(b)
        q, r = divrem2(a, b)
        a, b = b, r
        issimpler(b, a) || throw(DomainError((a,b), "b is not simpler than a"))
        s0, s1 = s1, s0 - q * s1
        t0, t1 = t1, t0 - q * t1
    end
    u = lcunit(a)
    if isone(a)
        a, s0, t0
    else
        a, s0, t0 = a / u, s0 / u, t0 / u
    end
end

function gcdx(a::Union{AbstractVector{T},NTuple{N,T}}) where {N,T<:Ring}
    n = length(a)
    n == 0 && return zero(T), zero(T), zero(T)
    u = similar(a)
    g = a[1]
    u[1] = one(T)
    for i = 2:n
        g, s, u[i] = gcdx(g, a[i])
        for j = 1:i-1
            u[j] *= s
        end
    end
    g, u
end
gcdx(a::T...) where T<:Ring = gcdx([a...])

# least common multiplier derived from gcd
function lcm(a::T, b::T) where T<:Ring
    div(a, gcd(a, b)) * b
end

function lcm(aa::AbstractVector{T}) where T<:Ring
    n = length(aa)
    n == 0 && return one(T)
    n == 1 && return aa[1]
    g = lcm(aa[1], aa[2])
    for i = 3:n
        g = lcm(aa[i], g)
    end
    g
end

function Base.powermod(x::R, p::Integer, m::R) where R<:Ring
    if p == 1
        return rem(x, m)
    elseif p == 0
        return one(x)
    elseif p == 2
        return rem(x * x, m)
    elseif p < 0
        if isunit(x)
            isone(x) && return copy(x)
            isone(-x) && return iseven(p) ? one(x) : copy(x)
            powermod(inv(x), -p, m)
        else
            throw(ArgumentError("negative powers not supported"))
        end
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x  = rem(x * x, m)
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x  = rem(x * x, m)
        end
        y = rem(y * x, m)
    end
    return y
end

const SUPERSCRIPTS = Char[0x2070, 0xb9, 0xb2, 0xb3, 0x2074, 0x2075, 0x2076, 0x2077, 0x2078, 0x2079, 0x207a, 0x207b]
const SUBSCRIPTS = Char[0x2080, 0x2081, 0x2082, 0x2083, 0x2084, 0x2085, 0x2086, 0x2087, 0x2088, 0x2089, 0x208a, 0x208b]

tosuper(a::Integer; sign::Bool=false) = _integer_to_script(a, SUPERSCRIPTS, sign)
tosub(a::Integer; sign::Bool=false) = _integer_to_script(a, SUBSCRIPTS, sign)
function _integer_to_script(a::Integer, chars::Vector{Char}, sign::Bool)
    io = IOBuffer()
    if a < 0
        print(io, chars[12])
    elseif sign
        print(io, chars[11])
    end
    for d in reverse(digits(a))
        print(io, chars[abs(d)+1])
    end
    String(take!(io))
end

function sort_unique!(A::AbstractVector; rev::Bool=false)
    n = length(A)
    n == 0 && return A
    a = sort!(A, rev=rev)
    j = 1
    aj = a[1]
    for i = 2:length(a)
        ai = a[i]
        if !isequal(aj, ai)
            j += 1
            aj = ai
        end
        if i != j
            a[j] = aj
        end
    end
    if j < n
        resize!(a, j)
    end
    a
end

"""
    testrules(io, g)

Test associativity and distributivity for all element combinations from `g`.
Print messages to `io` if errors found.
`g` can be a collection of `Ring` elements or a subtype of `Ring` (which is iterable).
"""
function testrules(io, gg)
    for a in gg
        if isunit(a)
            if !isone(inv(a) * a) || !isone(a * inv(a))
                println(io, "inv($a)")
            end
        end
        if a * a * a != a ^ 3
            println(io, "$(a) ^ 3")
        end
    end
    for (a,b,c) in Iterators.product(gg, gg, gg)
        if (a * b) * c != a * ( b * c)
            println(io, "assoc * $a $b $c")
        end
        if (a + b) + c != a + ( b + c)
            println(io, "assoc + $a $b $c")
        end
        if a * ( b + c) != a * b + a * c
            println(io, "distrib $a $b $c")
        end
    end
end


text:
"""
    BaseStats

Every stats common to any configuration.
"""
mutable struct BaseStats <: Statistics
    # federated learning stats
    num_comm_rounds::Int
    num_clients_per_round::Int
    num_total_clients::Int

    # machine learning stats
    num_weights::Int
    losses::Vector{Float32}
    accuracies::Vector{Float32}

    # changes inside data
    req_data::Any
    res_data::Vector{Any}

    BaseStats(
        num_comm_rounds::Int,
        fraction_clients::Float32,
        num_total_clients::Int,
        num_weights::Int
    ) = new(
        # fedetated learning
        num_comm_rounds,
        max(round(Int, fraction_clients * num_total_clients), 1),
        num_total_clients,
        # machine learning
        num_weights,
        Vector{Float32}(undef, 0),
        Vector{Float32}(undef, 0),
        # changes inside data
        Any,
        Vector{Any}(undef, 0)
    )
end

function update_stats!(
    stats::BaseStats,
    loss::Float32,
    accuracy::Float32,
)
    push!(stats.losses, loss)
    push!(stats.accuracies, accuracy)

    stats.req_data = Any # Vector(undef, 0)
    stats.res_data = Vector{Any}(undef, 0)
end


text:
module ControlHost

using Sockets
using StaticArrays

export @ip_str
export CHClient, CHTag, CHPrefix, CHMessage, subscribe

function num2byte(x::T) where {T<:Number}
    iob = IOBuffer()
    write(iob, x)
    seekstart(iob)
    return read(iob)
end

struct CHClient
    ip::IPv4
    port::UInt16
    socket::TCPSocket
    CHClient(ip, port) = begin
        socket = connect(ip, port)
        new(ip, port, socket)
    end
end

struct CHTag
    data::SVector{8, UInt8}
end

function CHTag(s::AbstractArray{T}) where {T<:UInt8}
    n = length(s)
    data = zeros(UInt8, 8)
    for i in 1:min(8, n)
        data[i] = UInt8(s[i])
    end
    CHTag(SVector{8, UInt8}(data))
end

CHTag(s::String) = CHTag(Vector{UInt8}(s))

CHTag(s::TCPSocket) = CHTag(read(s, 8))
String(t::CHTag) = String(Vector{Char}(t.data))

Base.show(io::IO, t::CHTag) = print(io, "CHTag: '" * String(t) * "'")

struct CHPrefix
    tag::CHTag
    length::UInt32
    data::SVector{16, UInt8}

    CHPrefix(tag, length::UInt32) = begin
        length_data = reverse(num2byte(length))
        data = vcat(tag.data, length_data, zeros(UInt8, 4))
        new(tag, length, data)
    end
end

function CHPrefix(s::TCPSocket)
    tag = CHTag(s)
    length = reinterpret(UInt32, reverse(read(s, 4)))[1]
    read(s, 4)  # dummy bytes
    CHPrefix(tag, length)
end

Base.show(io::IO, p::CHPrefix) = begin
    tag = String(p.tag)
    print(io, "CHPrefix with tag '$tag' and length $(p.length)")
end

struct CHMessage
    prefix::CHPrefix
    data::Vector{UInt8}
end

function CHMessage(s::TCPSocket)
    prefix = CHPrefix(s)
    length = prefix.length
    data = read(s, length)
    return CHMessage(prefix, data)
end

CHMessage(c::CHClient) = CHMessage(c.socket)

Base.show(io::IO, m::CHMessage) = begin
    tag = String(m.prefix.tag)
    print(io, "CHMessage with tag '$tag' and length $(m.prefix.length)")
end

function subscribe(c::CHClient, tag::AbstractString; mode::AbstractString="wait")
    chtag = CHTag("_Subscri")
    prefix = CHPrefix(chtag, UInt32(length(tag)+3))
    message = CHMessage(prefix, Vector{UInt8}(" $(mode[1]) $tag"))
    data = vcat(prefix.data, message.data)
    write(c.socket, data)
    write(c.socket, CHPrefix(CHTag("_Always"), UInt32(0x00)).data)
end

Base.close(c::CHClient) = close(c.socket)

end # module


text:
import Base:append!, getindex

mutable struct Container
    obj
    deps::AbstractArray
    function Container(obj, deps)
        deps_cont = filter(x->typeof(x) == Container, deps)
        deps_notcont = filter(x->typeof(x) != Container, deps)
        if !isempty(deps_cont) # add all previous container objects and dependencies to the new container
            return new(obj, [deps_notcont..., map((x->getfield(x, :obj)), deps_cont)..., (map((x->getfield(x, :deps)), deps_cont)...)...])
        else
            return new(obj, deps)
        end
    end
end

function append!(original::Container, new::Container...)
    append!(original.deps, (map((x->getfield(x, :deps)), new)...))
end

getindex(cont::Container) = cont.obj

function Base.show(io::IO, cont::Container)
    println("Container:")
    println("    Object: $(cont[])")
    println("    Dependencies:")
    for dep in cont.deps
        println("        $dep")
    end
end

"""
    troubleshoot(container)

Recursively loads all pointers.
"""
function troubleshoot(obj)
    if typeof(obj) <: Ptr
        o = unsafe_load(obj)
        typeof(o) <: Ptr ? troubleshoot(o) : @warn o
    elseif typeof(obj) <: Ref
        troubleshoot(obj[])
    elseif typeof(obj) <: Union{Tuple,AbstractArray}
        for o in obj
            troubleshoot(o)
        end
    end
end

# this function has to be generated so that there is no change in scope, which may invalide pointers
@generated function unsafe_pointer(obj)
    quote
        if typeof(obj) == Container
            unsafe_pointer(obj[])
        elseif typeof(obj) <: Union{Base.RefValue,Base.AbstractArray}
            Base.unsafe_convert(Ptr{typeof(obj[])}, obj)
        end
    end
end

macro safe_contain(expr)
    str_expr = repr(expr)
    unsafe_ptr_vars_regex = r"unsafe_pointer\((.*?)\)"
    vars = Any[Symbol(x.captures[1]) for x in eachmatch(unsafe_ptr_vars_regex, str_expr)]
    quote
        main_obj_cont = $(esc(expr))
        loc = Base.@locals()
        deps = Any[loc[symbol] for symbol in $vars]
        container = Container(main_obj_cont, deps)
    end
end


text:
# plot first 20,000 time steps (before stabilization)
include("../Model_scripts/Functions.jl");
# plot the number of infected & uninfected hosts
plot_population_numbers(
    hcat(labels...),
    "c=0.5\n \nNumber of infected and uninfected hosts",
    100,
    20000,
    "$figure_directory/supplemental/supp_3a.png"
)
# plot the average weighted R0
plot_reproductive_rate(5, 20000, "$figure_directory/supplemental/supp_3c.png")
# plot the average V0
plot_vertical_reproductive_ratio(1, 20000, "$figure_directory/supplemental/supp_3e.png")
# plot the average horizontal transmission
plot_horizontal_and_virulence(1, 20000, "$figure_directory/supplemental/supp_3g.png")
# plot the evenness through time
plot_evenness(1, 20000, "$figure_directory/supplemental/supp_3i.png")