text:
using LAP_julia, TimerOutputs, DataFrames, TableView, JLD2, FileIO
##
sp_lap_point_count = DataFrame(
    index = Int[],
    diag_size = Int[],
    whs = Int[],
    timer = TimerOutput[],
    median_results = Dict{String,Float64}[],
    point_count = Int[],
    spacing = Int[],
    points_found = Float16[]
    );

window_half_sizes = cat(collect(15:15:51), dims=1)
diag_sizes = [400, 600, 900, 1200];
point_counts = vcat(collect(50:150:700), collect(800:300:1550))
spacings = vcat(5, 15, collect(35:15:40), 50, collect(60:20:155))

##
using LAP_julia: fun_on_dict_values, resize_to_diag_size
using BenchmarkTools

img_count = 5;
errors = Array{Any}(undef, 0)
imgs = Array{Image}(undef, 5)
# prepare images:
for k in 1:img_count
    imgs[k], _ = gen_anhir(mutate=false)
end

imgshow(imgs[2])
imgshow(resize_to_diag_size(imgs[2], 450))
size(resize_to_diag_size(imgs[2], 500))

total_iterations = sum(broadcast(*, length.([window_half_sizes, diag_sizes, point_counts, spacings])...)) * img_count

# erase old_results
df = sp_lap_point_count;

let index = 0
    for diag_size in diag_sizes
        for whs in window_half_sizes
            for point_count in point_counts, spacing in spacings
                timer = TimerOutput("reg alg: sp lap")
                results_dicts = Array{Dict}(undef, img_count)
                inds_array = Array{Array{Any}}(undef, img_count)

                for k in 1:img_count
                    img = resize_to_diag_size(imgs[k], diag_size)
                    whs_limited = round(Int, min(whs, min(size(img)...)/4))
                    rand_flow = gen_quad_flow(size(img), whs_limited)
                    imgw = warp_img(img, real(rand_flow), imag(rand_flow))
                    try
                        # run registration function
                        _, _, timer, results_dicts[k], (_, inds_array[k]) = test_registration_alg(sparse_lap,
                            img, imgw, rand_flow, [whs_limited], Dict(:timer => timer, :point_count =>
                                point_count, :spacing => spacing), timer=timer, display=false)
                    catch e
                        println(e)
                        push!(errors, [e, img, imgw, rand_flow, whs_limited])
                        continue
                    end
                end
                # if no data gathered:
                if any(x-> 0 .== x, length.([results_dicts, inds_array]))
                    println("SKIPPED: diag_size:", diag_size, " whs: ", whs)
                    continue
                end

                median_results = fun_on_dict_values(filter_defined(results_dicts), median)
                avg_points_found = LAP_julia.mean(map(x -> length(x), filter_defined(inds_array)))

                index = index + 1
                println("at index: ", index, " (of $total_iterations", " diag_size: ", diag_size, " whs: ", whs)
                push!(df, Dict(:index => index,
                               :diag_size => diag_size,
                               :whs => whs,
                               :timer => timer,
                               :median_results => median_results,
                               :point_count => point_count,
                               :spacing => spacing,
                               :points_found => avg_points_found
                               ))
            end # point_count in point_counts, spacing in spacings
        end
    end
end

defined_perms(arraylike) = filter(k -> isassigned(arraylike, k), 1:length(arraylike))
filter_defined(arraylike) = arraylike[defined_perms(arraylike)]

results_dicts[defined_perms(results_dicts)]
filter_defined(results_dicts)

errors

results_dicts = Array{Dict}(undef, img_count)

results_dicts[2] = Dict("sd" => 12)

filter(k -> isassigned(results_dicts, k), 1:length(results_dicts))

isassigned(results_dicts, 2)

filter(, results_dicts)

# df

# win_sum1_df = df[(df.reg_fun .== :sparse_lap_win_sum1), :];
# win_sum3_df = df[(df.reg_fun .== :sparse_lap), :]; # [1, :timer]
#
# TimerOutputs.ncalls(df[(df.reg_fun .== :sparse_lap), :][1, :timer][sections[1]][sections[2]][sections[3]])
#
# get_timer(win_sum3_df[1, :timer], sections)

import LAP_julia: sparse_lap

win_sum3_df[1, :timer]

sections = ["reg alg: sp lap", "sparse lap", "prepare A and b"]
sections = ["reg alg: sp lap", "sparse lap", "filtering"]

fig, axs = subplots(6, sharex=true, figsize=(10,10))
for (diag_size, ax) in zip(diag_sizes[1:end-1], axs)
    ax.set_title("Image Size: " * string(diag_size))

    n = count(k -> k < (diag_size/4), window_half_sizes)
    timers3 = df[(df.reg_fun .== :sparse_lap) .& (df.diag_size .== diag_size), :][:, :timer][1:n]

    y3 = map(x -> TimerOutputs.time(get_timer(x, sections))/
                  (TimerOutputs.ncalls(get_timer(x, sections)) * 10e8), timers3)
    x = filter(x -> x < (diag_size/4), window_half_sizes)[1:n]

    println(length(y3), "  ", length(x))

    xscale("log")
    yscale("log")
    subplots_adjust(top=1.5)
    ax.plot(x, y3)
    xlabel("whs [pixels]")
    ax.set_ylabel("time [s]")
end
gcf()

showtable(df)

TimerOutputs.time(get_timer(win_sum3_df[:, :timer][1], ["reg alg: sp lap", "sparse lap", "prepare A and b"]))
TimerOutputs.tottime(get_timer(win_sum3_df[:, :timer][1], ["reg alg: sp lap", "sparse lap", "prepare A and b"]))
TimerOutputs.tottime(get_timer(win_sum3_df[:, :timer][1], ["reg alg: sp lap", "interpolate flow"]))
TimerOutputs.time(win_sum3_df[:, :timer][1]["reg alg: sp lap"]["interpolate flow"])

print_timer(win_sum3_df[:, :timer][1])

function get_timer(timer, sections)
    private_sections = copy(sections)
    try
        timer[private_sections[1]]
    catch e
        println(e)
        return TimerOutput()
    end
    out = timer
    while private_sections != []
        out = out[popfirst!(private_sections)]
    end
    return out
end

full_flow_estim, source_reg, flow_estim_at_inds, inds = sparse_lap_filt_one_by_one(img, imgw, 12)

inds

using ImageFiltering: KernelFactors.gaussian
using LAP_julia: inds_to_points, filt_onebyone!, window_sum3_at_inds!, multi_mat_div_qr
function sparse_lap_filt_one_by_one(target,
                                    source,
                                    fhs,
                                    window_size=(2fhs+1, 2fhs+1);
                                    spacing::Integer=15,
                                    point_count::Integer=100,
                                    timer::TimerOutput=TimerOutput("sparse_lap"),
                                    display::Bool=false,
                                    flow_interpolation_method::Symbol=:quad,
                                    base_method_kwargs=Dict(:timer => timer, :display => display))

    mask = parent(padarray(trues(size(target).-(2*fhs, 2*fhs)), Fill(false, (fhs, fhs), (fhs, fhs))))
    @timeit_debug timer "find edge points" begin
        inds = find_edge_points(target, spacing=spacing, number=point_count, mask=mask)
    end
    @timeit_debug timer "sparse lap" begin
        flow_estim_at_inds, inds = single_lap_at_points_filt_one_by_one(target, source, fhs, window_size, inds; base_method_kwargs...)
    end
    if all(isnan, flow_estim_at_inds)
        @timeit_debug timer "interpolate flow" begin
            full_flow_estim = zeros(size(target)) .* im .+ zeros(size(target))
        end
    else
        @timeit_debug timer "interpolate flow" begin
            full_flow_estim = interpolate_flow(flow_estim_at_inds, inds, size(target), method=flow_interpolation_method)
        end
    end
    @timeit_debug timer "generate source_reg" begin
        source_reg = warp_img(source, real(full_flow_estim), imag(full_flow_estim))
    end
    return full_flow_estim, source_reg, flow_estim_at_inds, inds
end

function single_lap_at_points_filt_one_by_one(target::Image,
                                              source::Image,
                                              filter_half_size::Integer,
                                              window,
                                              inds::Array{CartesianIndex,1};
                                              filter_count::Integer=3,
                                              timer::TimerOutput=TimerOutput("Sparse LAP"),
                                              display::Bool=false)

    image_size = size(target)
    pixel_count = length(target)
    ind_count = length(inds)
    filter_size = 2 * filter_half_size + 1

    # Prepare filters:
    # Calculate separable filters from basis:
    sigma = (filter_half_size + 2) / 4
    centered_inds = centered(-filter_half_size:filter_half_size)
    gaus = gaussian(sigma, filter_size)
    gaus_der = gaus .* centered_inds .* (-1)/sigma^2
    gaus_der_flip = reverse(gaus_der, dims=1)

    # prepare 2D filter basis for later use
    basis = similar(target, (filter_size, filter_size, filter_count))

    # dif = (points of forward image - points of backward image) for each filter
    dif = similar(target, (image_size..., filter_count))
    points = inds_to_points(inds)

    @timeit_debug timer "filtering" begin
        if filter_count == 3
            # temporary place to store filtered images
            # tmp_filtered_1 = fill(NaN, image_size)
            # tmp_filtered_2 = fill(NaN, image_size)
            tmp_filtered_1 = similar(target)
            tmp_filtered_2 = similar(target)

            # basis 1 - (gaus, gaus) both forward and backward
            kernf_1 = kernelfactors((gaus, gaus))
            basis[:, :, 1] = broadcast(*, kernf_1...)
            filt_onebyone!(tmp_filtered_1, target, kernf_1, filter_half_size, points)
            filt_onebyone!(tmp_filtered_2, source, kernf_1, filter_half_size, points)
            # imfilter!(tmp_filtered_1, target, kernf_1, "symmetric")
            # imfilter!(tmp_filtered_2, source, kernf_1, "symmetric")
            dif[:, :, 1] = tmp_filtered_2 - tmp_filtered_1

            # basis 2 - (gaus, gaus_der_flip) as forward, (gaus, gaus_der) as backward
            kernf_2f = kernelfactors((gaus, gaus_der_flip))
            basis[:, :, 2] = broadcast(*, kernf_2f...)
            filt_onebyone!(tmp_filtered_1, target, kernf_2f, filter_half_size, points)
            # imfilter!(tmp_filtered_1, target, kernf_2f, "symmetric")
            kernf_2b = kernelfactors((gaus, gaus_der))
            filt_onebyone!(tmp_filtered_2, source, kernf_2b, filter_half_size, points)
            # imfilter!(tmp_filtered_2, source, kernf_2b, "symmetric")
            dif[:, :, 2] = tmp_filtered_2 - tmp_filtered_1

            # basis 3 - (gaus_der_flip, gaus) as forward, (gaus_der, gaus) as backward
            kernf_3f = kernelfactors((gaus_der_flip, gaus))
            basis[:, :, 3] = broadcast(*, kernf_3f...)
            filt_onebyone!(tmp_filtered_1, target, kernf_3f, filter_half_size, points)
            # imfilter!(tmp_filtered_1, target, kernf_3f, "symmetric")
            kernf_3b = kernelfactors((gaus_der, gaus))
            filt_onebyone!(tmp_filtered_2, source, kernf_3b, filter_half_size, points)
            # imfilter!(tmp_filtered_2, source, kernf_3b, "symmetric")
            dif[:, :, 3] = tmp_filtered_2 - tmp_filtered_1
        end
    end # "filtering"

    dif = reshape(dif, (:, filter_count))

    # prepare matrices for linear system of equations
    A = similar(target, (filter_count-1, filter_count-1, ind_count))
    b = similar(target, (filter_count-1, ind_count))

    @timeit_debug timer "prepare A and b" begin
        # TODO: check whether window_sum_around_points calculates only around points and not on whole window.
        for k in 1:filter_count-1
            for l in k:filter_count-1
                @timeit_debug timer "window sum part 1" begin
                    # @views window_sum_around_points!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window, points)
                    @views window_sum3_at_inds!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window, inds)
                    A[l, k, :] = A[k, l, :]
                end
            end
            @timeit_debug timer "window sum part 2" begin
                # @views window_sum_around_points!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window, points)
                @views window_sum3_at_inds!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window, inds)
            end
        end
    end # "prepare A and b"

    # Perform Gauss elimination on all pixels in parallel:
    # coeffs will be of shape: pixel_count, filter_count-1
    @timeit_debug timer "multi mat div" begin
        @views coeffs = multi_mat_div_qr(A, b)
    end

    # adding ones so that all base filters have their coefficients even the first one
    all_coeffs = [ones(ind_count) coeffs]

    k = (-filter_half_size:filter_half_size)

    # Get the displacement vector field from the filters
    @timeit_debug timer "calculate flow" begin
        u_est_at_inds = similar(target, ComplexF64, ind_count)
        @views u_est_at_inds[:] = 2 .* ((im .* (.-1 .* all_coeffs[:, 3]) ./ all_coeffs[:, 1]) .+ ((.-1 .* all_coeffs[:, 2]) ./ all_coeffs[:, 1]));
    end # "calculate flow"

    # dont use estimations whose displacement is larger than the filter_half_size
    displacement_mask = (real(u_est_at_inds).^2 .+ imag(u_est_at_inds).^2) .<= filter_half_size^2

    return u_est_at_inds[displacement_mask], inds[displacement_mask]

end


text:
# Rational number reconstruction implementation borrowed from CLUE
# and modified a bit to suit the 'Modern Computer Algebra' definitions
# Returns a rational r // h of QQ field in a canonical form such that
#   r // h ≡ a (mod m)
#
# let n = max( λ(a), λ(m) ) , where λ(x) is a number of bits for x
# O(n^2)
function rational_reconstruction(a::I, m::I) where {I<:Union{Int, BigInt}}
    bnd = sqrt(float(m) / 2)

    steps = 0

    # @warn "" bnd
    U = (I(1), I(0), m)
    V = (I(0), I(1), a)
    while abs(V[3]) >= bnd
        q = div(U[3], V[3])
        T = U .- q .* V
        U = V
        V = T

        steps += 1
    end

    t = abs(V[2])
    r = V[3] * sign(V[2])

    # @warn "" r t

    # changed from `<= bnd` to `<= m / bnd`
    # we can speed up this !
    if t <= bnd && gcd(r, t) == 1
        # return (r, t)
        return (true, steps)
    end

    return (false, steps)
end

function maximal_quotient_rational_reconstruction(a::I, m::I) where {I<:Union{Int, BigInt}}

    steps = 0

    T = 2^20 * log2(m)

    n, d = 0, 0

    t0, r0 = 0, m
    t1, r1 = 1, a

    U = (I(1), I(0), m)
    V = (I(0), I(1), a)
    while r1 != 0 && r0 > T
        q = div(r0, r1)

        if q > T
            n, d, T = r1, t1, q
        end

        steps += 1

        r0, r1 = r1, r0 - q*r1
        t0, t1 = t1, t0 - q*t1
    end

    if d == 0 || gcd(n, d) != 1
        return (false, steps)
    end

    if d < 0
        n, d = -n, -d
    end

    return (true, steps)
end

ext: jl


text:
include("calculator.yy.jl")

src = "2 x (5 + 3) + 2 * 2"

println("---")
result = parse(src)
println("Result : $result")


text:
"""
`DirImageSource(path; nargs...)` → `DirImageSource`

Description
============

Concrete implementation of `ImageSource` for (nested) directories
of a local filesystem. The purpose of `DirImageSource` is to
seamlessly load images from a directory and its subdirectories.
Note that subdirectories will be treated as categories and their
relative paths will be stored as metadata within the loaded image
object.

Every subtype of `ImageSource` implements most read-only functions
of the array interface. Additionally, it fully implements the
iterator interface and can be used as for-loop to iterate over
all the stored images.

Usage
======

    DirImageSource(path = "."; hidden=false, expand=false, recursive=true, formats=[".png", ".jpg", ".jpeg", ".bmp"])

Arguments
==========

- **`path`** : The directory that contains the desired images.
This can also be just a container for subdirectories, in which
the subdirectories contain the actual images. However, for this
to work `recursive` has to be set to `true`

- **`hidden`** : If `true`, all the hidden files within `path`
will be processed as well.

- **`expand`** : If `true`, all the paths will be extended to
absolute paths instead of being relative to the root directory
specified by `path`.
It is generally recommended to set `expand = false` (default).

- **`recursive`** : If `true`, then all the subdirectories of
`path` will be processed as well. That implies that if any
subdirectory, or their subdirectories, contain any images of
a format specified by `formats`, then those images will be part
of the `DirImageSource`. Defaults to `true`

- **`formats`** : Array of strings. Specifies which file endings
should be considered an image. Any file of such ending will be
available as part of the `DirImageSource`

Methods
========

- **`rand`** : Loads and returns a random image, or an array of
images, from the specified directory (depending on the parameters).
If `recursive` was set to `true`, then all the images from the
subdirectories are considered as well.

- **`length`** : Returns the total number of images that are
available in the specified directory. If `recursive` was set
`true`, then all the images from the subdirectories are counted
as well

- **`getindex`** : Loads and returns the image of the given index.
The numerical index is more or less arbitrary and depends on the
order the images were parsed.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    # create the image source by indexing the directory's content.
    # Note: this command won't actually load any images.
    src = DirImageSource("mydatadir/")

    # number of available images
    n = length(src)

    # images can be indexed like an array
    img1 = src[1]

    # loop through all images in the source
    for img in src
        println(img)
    end

see also
=========

`rand`, `length`

`ImageSource`
"""
immutable DirImageSource <: ImageSource
    path::UTF8String
    files::Vector{UTF8String}
end

function _name(filename)
    ending = match(r"\.[a-zA-Z]+$", filename)
    if ending == nothing || ending.offset == 1
        return filename
    else
        filename[1:(ending.offset-1)]
    end
end

function DirImageSource(path = "."; expand = false, nargs...)
    filepaths = listfiles(path; expand = expand, nargs...)
    filenames = map(_ -> relpath(_, path), filepaths)
    DirImageSource(expand ? abspath(path) : path, filenames)
end

function Base.getindex(s::DirImageSource, index::Integer)
    filename = s.files[index]
    img_dir = dirname(filename)
    img_dir = length(img_dir) == 0 ? "." : img_dir
    img_name = _name(relpath(filename, img_dir))
    img = load(joinpath(s.path, filename))::Image
    img["name"] = img_name
    img["category"] = img_dir
    img
end

Base.length(s::DirImageSource) = length(s.files)
Base.start(::DirImageSource) = 1
Base.next(s::DirImageSource, state) = (s[state], state+1)
Base.done(s::DirImageSource, state) = state > length(s)


text:
using LightGraphs.Experimental.Traversals

"""
    struct BFS <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use the [Breadth-First Search algorithm](https://en.m.wikipedia.org/wiki/Breadth-first_search).

An optional sorting algorithm may be specified (default = no sorting).
Sorting helps maintain cache locality and will improve performance on
very large graphs; for normal use, sorting will incur a performance
penalty.

`BFS` is the default algorithm used when a source is specified
but no distance matrix is specified.

### Implementation Notes
`BFS` supports the following shortest-path functionality:
- (optional) multiple sources
- all destinations
"""
struct BFS{T} <: ShortestPathAlgorithm
    traversal::T
end

BFS() = BFS(Traversals.BFS())
BFS(a::Base.Sort.Algorithm) = BFS(Traversals.BFS(a))

mutable struct BFSSPState{U} <: Traversals.AbstractTraversalState
    parents::Vector{U}
    dists::Vector{U}
    n_level::U
end

@inline function initfn!(s::BFSSPState, u)
    s.dists[u] = 0
    return true
end
@inline function newvisitfn!(s::BFSSPState, u, v)
    s.dists[v] = s.n_level
    s.parents[v] = u
    return true
end
@inline function postlevelfn!(s::BFSSPState{U}) where U
    s.n_level += one(U)
    return true
end

struct BFSResult{U<:Integer} <: ShortestPathResult
    parents::Vector{U}
    dists::Vector{U}
end

function shortest_paths(
    g::AbstractGraph{U},
    ss::AbstractVector{U},
    alg::BFS
   ) where U <: Integer

    n = nv(g)
    dists = fill(typemax(U), n)
    parents = zeros(U, n)
    state = BFSSPState(parents, dists, one(U))
    Traversals.traverse_graph!(g, ss, alg.traversal, state)
    return BFSResult(state.parents, state.dists)
end

shortest_paths(g::AbstractGraph{U}, ss::AbstractVector{<:Integer}, alg::BFS) where {U<:Integer} = shortest_paths(g, U.(ss), alg)
shortest_paths(g::AbstractGraph{U}, s::Integer, alg::BFS) where {U<:Integer} = shortest_paths(g, Vector{U}([s]), alg)


text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1]])

dec = 64
modulus = 4800
base_frequency = 125000000
periods_per_step = 5
samples_per_period = div(modulus, dec)
periods_per_frame = 50 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency
slow_dac_periods_per_frame = div(50, periods_per_step)

stopADC(rp)
masterTrigger(rp, false)
clearSequence(rp)

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(master(rp), true)

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
signalTypeDAC(rp, 1 , "SINE")
amplitudeDAC(rp, 1, 1, 0.1)
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "EXTERNAL")

# Sequence
numSlowDACChan(master(rp), 1)
slowDACStepsPerFrame(rp, slow_dac_periods_per_frame)

fig = figure(1)
clf()

amplitudeDAC(rp, 1, 1, 0.1)
#First sequence
lut1 = [0.2]
seq1 = ConstantSequence(lut1, nothing, slow_dac_periods_per_frame, 2, computeRamping(master(rp), frame_period * 1, 0.5))

#Second sequence
lut2 = collect(range(0,0.7,length=slow_dac_periods_per_frame))
seq2 = ArbitrarySequence(lut2, nothing, slow_dac_periods_per_frame, 2, computeRamping(master(rp), 0.0, 0.0))

# Prepare both
appendSequence(master(rp), seq1)
appendSequence(master(rp), seq2)
success = prepareSequence(master(rp))
startADC(rp)
masterTrigger(rp, true)
sleep(0.1)

uCurrentFrame = readFrames(rp, 0, 7)
stopADC(rp)
masterTrigger(rp, false)
clearSequence(rp)

plot(vec(uCurrentFrame[:, 1, :, :]))
title("Two sequences")
fig

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
using HologramX
using Test

@testset "HologramX.jl" begin
    # Write your tests here.
end


text:
export play_sound, play_music

# NOTE this function is all we need for pong
function play_sound(sound_file, loops=0, ticks=-1)
  sample=Mix_LoadWAV_RW(SDL_RWFromFile(sound_file, "rb"), 1);
  if sample == C_NULL
      @warn "Could not load sound file: $sound_file\n$(getSDLError())"
      return
  end
  r = Mix_PlayChannelTimed(Int32(-1), sample, loops, ticks)
  if r == -1
      @warn "Unable to play sound $sound_file\n$(getSDLError())"
  end
end


text:
"Rejection Sampling"
struct RejectionSampleAlg <: SamplingAlgorithm end

"Rejection Sampling"
const RejectionSample = RejectionSampleAlg()

isapproximate(::RejectionSampleAlg) = false

"`n` samples from `x` with rejection sampling"
function Base.rand(rng,
                   ΩT::Type{OT},
                   pred::RandVar,
                   n::Integer,
                   alg::RejectionSampleAlg) where {OT <: Ω}
  ωsamples = OT[]
  accepted = 0
  i = 1
  while accepted < n
    ω = ΩT() # FIXME, use rng to select random points
    issat = pred(Omega.Space.tagrng(ω, rng))
    if issat
      push!(ωsamples, ω)
      accepted += 1
    end
    # cb((ω = ω, accepted = accepted, p = float(issat), i = i), IterEnd)
    lens(Loop, (ω = ω, accepted = accepted, p = float(issat), i = i))
    i += 1
  end
  ωsamples
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer,
                   alg::RejectionSampleAlg;
                   ΩT::Type{OT} = defΩ(alg),
                   memoize = true) where {OT <: Ω}
  pred = Omega.mem(Omega.indomain(x))
  ωsamples = rand(rng, ΩT, pred, n, alg)
  map(Omega.mem(x), ωsamples)
end


text:
# Routines for interaction with pymatgen, e.g. converting to
# its structures from the DFTK equivalents

function pymatgen_lattice(lattice::AbstractArray)
    # Notice: Pymatgen uses rows as lattice vectors, so we unpeel
    # our lattice column by column. The default unit in pymatgen is Ǎngström
    mg = pyimport("pymatgen")
    bohr_to_A = 1 / units.Ǎ
    mg.Lattice([Array(bohr_to_A .* lattice[:, i]) for i in 1:3])
end
pymatgen_lattice(model::Model) = pymatgen_lattice(model.lattice)

function pymatgen_structure(model_or_lattice, atoms)
    mg = pyimport("pymatgen")
    pylattice = pymatgen_lattice(model_or_lattice)

    n_species = sum(length(pos) for (spec, pos) in atoms)
    pyspecies = Vector{Int}(undef, n_species)
    pypositions = Array{Vector{Float64}}(undef, n_species)
    ispec = 1
    for (spec, pos) in atoms
        for coord in pos
            pyspecies[ispec] = spec.Z
            pypositions[ispec] = Vector{Float64}(coord)
            ispec = ispec + 1
        end
    end
    @assert ispec == n_species + 1

    mg.Structure(pylattice, pyspecies, pypositions)
end
pymatgen_structure(model::Model) = pymatgen_structure(model, model.atoms)

function pymatgen_bandstructure(basis, λ, εF, klabels)
    Spin = pyimport("pymatgen.electronic_structure.core").Spin
    bandstructure = pyimport("pymatgen.electronic_structure.bandstructure")

    # This assumes no spin polarization
    @assert basis.model.spin_polarization in (:none, :spinless)

    eigenvals_spin_up = Matrix{eltype(λ[1])}(undef, length(λ[1]), length(basis.kpoints))
    for (ik, λs) in enumerate(λ)
        eigenvals_spin_up[:, ik] = λs
    end
    eigenvals = Dict(Spin.up => eigenvals_spin_up)

    kcoords = [kpt.coordinate for kpt in basis.kpoints]
    pylattice = pymatgen_lattice(basis.model.lattice)
    bandstructure.BandStructureSymmLine(
        kcoords, eigenvals, pylattice.reciprocal_lattice, εF,
        labels_dict=klabels, coords_are_cartesian=true
    )
end

function load_lattice_pymatgen(T, pyobj::PyObject)
    mg = pyimport("pymatgen")

    if pyisinstance(pyobj, mg.Structure)
        load_lattice_pymatgen(T, pyobj.lattice)
    elseif pyisinstance(pyobj, mg.Lattice)
        lattice = Matrix{T}(undef, 3, 3)
        for i in 1:3, j in 1:3
            lattice[i, j] = units.Ǎ * get(get(pyobj.matrix, j-1), i-1)
        end
        Mat3{T}(lattice)
    else
        error("load_lattice_pymatgen not implemented for python type $pyobj")
    end
end

"""
Load a DFTK-compatible atoms representation from a supported pymatgen object.
All atoms are using a Coulomb model.
"""
function load_atoms_pymatgen(T, pyobj::PyObject)
    @assert pyisinstance(pyobj, pyimport("pymatgen").Structure)
    map(unique(pyobj.species)) do spec
        coords = [s.frac_coords for s in pyobj.sites if s.specie == spec]
        ElementCoulomb(spec.number) => coords
    end
end


text:
for T in (:(AlignedSequence),
    :(AnnotatedAlignedSequence),
    :(MultipleSequenceAlignment),
    :(AnnotatedMultipleSequenceAlignment))
    @eval Base.IndexStyle(::Type{$(T)}) = Base.IndexLinear()
end

@inline Base.getindex(x::AbstractAlignedObject,
              args...) = getindex(namedmatrix(x), args...)

@inline function Base.setindex!(x::AbstractAlignedObject, value, args...)
    setindex!(namedmatrix(x), value, args...)
end

# Special getindex/setindex! for sequences to avoid `seq["seqname","colname"]`

@inline function Base.getindex(x::AbstractAlignedSequence,
                       i::Union{Int,AbstractString})
    getindex(namedmatrix(x), 1, i)
end

@inline function Base.setindex!(x::AbstractAlignedSequence, value,
                        i::Union{Int,AbstractString})
    setindex!(namedmatrix(x), value, 1, i)
end

# ## Special getindex to conserve annotations and types

"""
Helper function to create a boolean mask to select sequences.
"""
function _get_selected_sequences(msa, selector)
    type = eltype(selector)
    if type == Bool
        return selector
    else
        to_select = Set(selector)
        if type  <: Number
            Bool[i in to_select for i in 1:nsequences(msa)]
        elseif type <: AbstractString
            Bool[i in to_select for i in sequencenames(msa)]
        end
    end
end

_get_names_to_index(msa::AbstractAlignedObject) = NamedArrays.names(msa.matrix, 2)
_get_names_to_index(msa::NamedResidueMatrix{T}) where T = NamedArrays.names(msa, 2)
function _get_names_to_index(msa::Matrix{Residue})
    throw(ErrorException("`Matrix{Residue}` cannot be indexed with `AbstractString`s."))
end

function _column_indices(msa, selector)
    if eltype(selector) <: AbstractString
        colnames = _get_names_to_index(msa)
        return Int[findfirst(==(i), colnames) for i in selector]
    end
    selector
end

function Base.getindex(msa::AnnotatedMultipleSequenceAlignment,
               seqs::AbstractArray, cols::AbstractArray)
    annot = copy(annotations(msa))
    seq_selector = _get_selected_sequences(msa, seqs)
    filtersequences!(annot, sequencenames(msa), seq_selector)
    _annotate_seq_modification!(annot, seq_selector)
    col_selector = _column_indices(msa, cols)
    filtercolumns!(annot, col_selector)
    _annotate_col_modification!(annot, col_selector)
    AnnotatedMultipleSequenceAlignment(msa.matrix[seqs, cols], annot)
end

function Base.getindex(msa::AnnotatedMultipleSequenceAlignment,
               seqs::AbstractArray, cols::Colon)
    annot = copy(annotations(msa))
    seq_selector = _get_selected_sequences(msa, seqs)
    filtersequences!(annot, sequencenames(msa), seq_selector)
    _annotate_seq_modification!(annot, seq_selector)
    AnnotatedMultipleSequenceAlignment(msa.matrix[seqs, cols], annot)
end

function Base.getindex(msa::AnnotatedMultipleSequenceAlignment,
               seqs::Colon, cols::AbstractArray)
    annot = copy(annotations(msa))
    col_selector = _column_indices(msa, cols)
    filtercolumns!(annot, col_selector)
    _annotate_col_modification!(annot, col_selector)
    AnnotatedMultipleSequenceAlignment(msa.matrix[seqs, cols], annot)
end

Base.getindex(msa::AnnotatedMultipleSequenceAlignment, seqs::Colon, cols::Colon) = copy(msa)

function Base.getindex(msa::MultipleSequenceAlignment,
seqs::Union{AbstractArray,Colon},
cols::Union{AbstractArray,Colon})
    MultipleSequenceAlignment(msa.matrix[seqs, cols])
end

function Base.getindex(seq::AnnotatedAlignedSequence, cols::AbstractArray)
    seq_copy = copy(seq)
    col_selector = _column_indices(msa, cols)
    filtercolumns!(seq_copy, col_selector)
    _annotate_col_modification!(seq_copy, col_selector)
    seq_copy
end

function Base.getindex(seq::AlignedSequence, cols::Union{AbstractArray,Colon})
    AlignedSequence(seq.matrix[cols])
end


text:
# 0.18 deprecations. Remove before releasing 0.19
@deprecate path(t::Trie, str::AbstractString) partial_path(t::Trie, str::AbstractString)
@deprecate find_root find_root!


text:
using Test
using ScottishTaxBenefitModel
using .ModelHousehold: Household, Person, People_Dict, default_bu_allocation

using .ExampleHouseholdGetter
using .Definitions
using .STBIncomes
using .NationalInsuranceCalculations:
    calculate_national_insurance!,
    calc_class1_secondary

using .STBParameters: NationalInsuranceSys,weeklyise!
using .GeneralTaxComponents: WEEKS_PER_YEAR
using .FRSHouseholdGetter: get_household
using .Results: IndividualResult, map_incomes
using .ExampleHelpers

@testset "Run on actual Data" begin
    nhhs,npeople = init_data()
    nisys = NationalInsuranceSys{Float64}()
    weeklyise!( nisys )
    hh = get_household(27)
    person =  hh.people[120150022701]
    println( person )
    pres = IndividualResult{Float64}()
    calculate_national_insurance!( pres, person, nisys )
    @test pres.ni.class_1_primary >= 0.0
    @test pres.ni.class_1_secondary >= 0.0

end

@testset "Run on actual Data" begin
    nhhs,npeople = init_data()
    nisys = NationalInsuranceSys{Float64}()
    weeklyise!( nisys )
    for hno in 1:nhhs
        hh = get_household(hno)
        for (pid,person) in hh.people
            println( "on hh $hno; pid=$pid")
            pres = IndividualResult{Float64}()
            pres.income = map_incomes( person )
            calculate_national_insurance!( pres, person, nisys )
            @test pres.ni.class_1_primary >= 0.0
            @test pres.ni.class_1_secondary >= 0.0
        end # people loop
    end # hhld loop
end #

@testset "Melville 2019 ch16 examples 1; Class 1 NI" begin
    # BASIC IT Calcaulation on
    nisys = NationalInsuranceSys{Float64}()
    weeklyise!( nisys )
    @time names = ExampleHouseholdGetter.initialise()
    income = [110.0,145.0,325,755.0,1_000.0]
    nidue = [(0.0,false),(0.0,true),(19.08,true),(70.68,true),(96.28,true)]
    niclass1sec = [0.0,0.0,21.94,81.28,115.09]
    ntests = size(income)[1]
    @test ntests == size( nidue )[1]
    hh = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = hh.people[RUK_PERSON]
    pers.age = 50
    for i in 1:ntests
        # pers.income[wages] = income[i]
        println( "case $i income = $(income[i])")
        pres = IndividualResult{Float64}()
        pres.income[WAGES] = income[i]
        calculate_national_insurance!( pres, pers, nisys )
        class1sec = calc_class1_secondary( income[i], pers, nisys )
        @test pres.ni.class_1_primary ≈ nidue[i][1]
        @test pres.ni.above_lower_earnings_limit == nidue[i][2]
        @test round(class1sec,digits=2) ≈ niclass1sec[i]
        print( pres.ni )
    end
end

@testset "Melville 2019 ch16 examples 6,7; Class 2,4 NI" begin
    # BASIC IT Calcaulation on
    nisys = NationalInsuranceSys{Float64}()
    weeklyise!( nisys )
    # self employment testing
    hh = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = hh.people[RUK_PERSON]
    pers.age = 50
    pers.employment_status = Full_time_Self_Employed
    seinc = [6_280.0, 7_200.0]
    class2 = [0.0,3.0]
    nisys.class_2_threshold *= WEEKS_PER_YEAR
    nisys.class_4_bands *= WEEKS_PER_YEAR
    for i in 1:size(seinc)[1]
        pres = IndividualResult{Float64}()
        pres.income[SELF_EMPLOYMENT_INCOME] = seinc[i]
        println( "case $i seinc = $(seinc[i])")
        calculate_national_insurance!( pres, pers, nisys )
        @test pres.ni.class_2 ≈ class2[i]
    end
    seinc = [15_140.0, 55_000.0,7_500.0]
    class4 = [585.72, 3_823.12, 0.0 ]
    for i in 1:size(seinc)[1]
        pres = IndividualResult{Float64}()
        pres.income[SELF_EMPLOYMENT_INCOME] = seinc[i]
        println( "case $i seinc = $(seinc[i])")
        calculate_national_insurance!( pres, pers, nisys )
        @test pres.ni.class_4 ≈ class4[i]
    end

end # example 6,7


text:
module TestConfiguration

using Test
using TreeParzen

# threshold
@test_throws ArgumentError Config(-1.0, 25, 24, 20, 1.0)
@test_throws ArgumentError Config(2.0, 25, 24, 20, 1.0)

# linear_forgetting
@test_throws ArgumentError Config(0.25, -1, 24, 20, 1.0)

# draws
@test_throws ArgumentError Config(0.25, 25, -1, 20, 1.0)

# random_trials
@test_throws ArgumentError Config(0.25, 25, 24, -1, 1.0)

# prior_weight
@test_throws ArgumentError Config(0.25, 25, 24, -1, -1.0)
@test_throws ArgumentError Config(0.25, 25, 24, -1, 2.0)

end
true


text:
# The cube, 41063625 (345^3), can be permuted to produce two other cubes:
# 56623104 (384^33) and 66430125 (405^33). In fact, 41063625 is the smallest
# cube which has exactly three permutations of its digits which are also cube.
#
# Find the smallest cube for which exactly five permutations of its digits are
# cube.

using ProjectEulerSolutions

# For each fixed number of digits, search through the right ranges of numbers
# to build up sets of cubes.  Store them as sorted tuples, counting the number
# that match to the same key.  When we hit 5, return the smallest.
function p062solution(n::Integer=5)::Integer

    ndig_max = ndigits(2^63-1) - 1  # Bigger than 18 and we will need to use BigInts
    for ndig in 9:ndig_max

        dict = Dict{Tuple, Array{Integer,1}}()

        num_low = Integer(ceil(10^((ndig-1)/3)))
        num_high = Integer(floor(10^(ndig/3)))
        max_count = 1
        max_key = 1

        for i in num_low:num_high
            num = i^3
            key = Tuple(sort(digits(num)))
            if haskey(dict, key)
                append!(dict[key], num)
            else
                dict[key] = [num]
            end
            c = length(dict[key])
            if max_count < c
               max_count = c
               if max_count == n
                   return minimum(dict[key])
               end
            end
        end
    end
    return -1
end

p062 = Problems.Problem(p062solution)

Problems.benchmark(p062, 5)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
import Base

const NOISE_SIGMA = 15.0
const INITIAL_VARIANCE = 4NOISE_SIGMA^2
# const MIN_VARIANCE = NOISE_SIGMA^2
const MIN_VARIANCE = 5.0^2
const BACKGROUND_RATIO = 0.7
const INITIAL_WEIGHT = .05
const INITIAL_SORTKEY = INITIAL_WEIGHT / sqrt(INITIAL_VARIANCE)
const VARIANCE_THRESHOLD = 2.5^2
const INITIALIZATION_WINDOW = 200
const MIN_LEARNING_RATE = 1.0 / INITIALIZATION_WINDOW

struct Gaussian
    μ :: Float64
    σ² :: Float64
    ω :: Float64
    sortkey :: Float64
end
Gaussian(μ, σ², ω) = Gaussian(μ, σ², ω, ω / sqrt(σ²))
Gaussian(μ, σ²) = Gaussian(μ, σ², INITIAL_WEIGHT)
Gaussian(μ) = Gaussian(μ, INITIAL_VARIANCE, INITIAL_WEIGHT, INITIAL_SORTKEY)
Base.isless(K::Gaussian, L::Gaussian) = K.sortkey < L.sortkey
neglog_likelihood(K::Gaussian, x) = (x - K.μ)^2 / K.σ²

Base.:*(K::Gaussian, x) = Gaussian(K.μ, K.σ², K.ω * x, K.sortkey * x)
Base.:*(x, K::Gaussian) = K * x
Base.:/(K::Gaussian, x) = Gaussian(K.μ, K.σ², K.ω / x, K.sortkey / x)

struct MixtureParams
    threshold :: Float64
    background_threshold :: Float64
    learning_rate :: Float64
    number_frames :: Int
end

MixtureParams() = MixtureParams(VARIANCE_THRESHOLD, BACKGROUND_RATIO, 1.0, 0)

function increment_frames(p::MixtureParams)
    if p.number_frames > INITIALIZATION_WINDOW
        return MixtureParams(p.threshold, p.background_threshold,
                             MIN_LEARNING_RATE, p.number_frames+1)
    else
        return MixtureParams(p.threshold, p.background_threshold,
                             1.0 / (p.number_frames + 2), p.number_frames+1)
    end
end

mutable struct MixtureModel
    kernels::Array{Gaussian, 3}
    n::Int
    params::MixtureParams
end

function MixtureModel(number_gaussians::Int, rows::Int, cols::Int, μ=0.0, σ²=INITIAL_VARIANCE)
    kernels = Array{Gaussian}(number_gaussians, rows, cols)
    weight = 1.0 / number_gaussians
    kernels .= Gaussian(μ, σ², weight)
    MixtureModel(kernels, number_gaussians, MixtureParams())
end
MixtureModel(n::Int, m::Int) = MixtureModel(5, n, m)

#Base.size(m :: MixtureModel) = Base.size(m.kernels)
#Base.IndexStyle(::Type{MixtureModel}) = Base.IndexLinear()
#Base.getindex(M::MixtureModel, i::Int) = Base.getindex(M.kernels, i)
#Base.setindex!(M::MixtureModel, kernel::Gaussian, i::Int) = Base.setindex!(M.kernels, kernel, i)

function normalize!(M, idx)
    total_weight = 0.0
    n = size(M.kernels)[1]
    for i in 1:n
        total_weight += M.kernels[i, idx].ω
    end
    for i in 1:n
        M.kernels[i, idx] = M.kernels[i, idx] / total_weight
    end
    return M
end

function apply!(M::MixtureModel, input, output)
    for idx in CartesianRange(size(input))
        x = input[idx]
        output[idx] = _apply!(M, idx, x)
    end
    M.params = increment_frames(M.params)
    return output
end

function apply!(M::MixtureModel, input)
    output = similar(input)
    return apply!(M, input, output)
end

function _apply!(M::MixtureModel, idx, x)
    k = 0
    n = size(M.kernels)[1]
    matched = false
    τ = M.params.threshold
    λ = M.params.learning_rate
    for k in 1:n
        kernel = M.kernels[k, idx]
        if kernel.ω < eps() break end
        if neglog_likelihood(kernel, x) < τ
            k = update!(M, k, idx, x, λ)
            matched = true
            break
        end
    end
    if !matched
        kernel = Gaussian(x)
        k = replace_kernel_at!(k, idx, M, kernel)
    end
    weight_sum = 0.0
    average = 0.0
    for j in 1:k-1
        weight_sum += M.kernels[j, idx].ω
        average += M.kernels[j, idx].ω * M.kernels[j, idx].μ
    end
    if weight_sum < M.params.background_threshold
        return zero(x)
    elseif x < average
        return zero(x)
    else
        return one(x)
    end
end

function bg_energy!(M::MixtureModel, input)
    output = Array{Float64}(size(input))
    for idx in CartesianRange(size(input))
        x = input[idx]
        output[idx] = _bg_energy!(M, idx, x)
    end
    M.params = increment_frames(M.params)
    return output
end

function _bg_energy!(M, idx, x)
    energy = Inf
    weight_sum = 0.0
    average = 0.0
    _apply!(M, idx, x)
    n = size(M.kernels)[1]
    for i in 1:n
        kernel = M.kernels[i, idx]
        weight_sum += kernel.ω
        average += kernel.ω * kernel.μ
        energy = min(neglog_likelihood(kernel, x), energy)
        if weight_sum >= M.params.background_threshold
            break
        end
    end
    if x < average
        return min(energy, 0.9 * M.params.threshold)
    else
        return energy
    end
end

function update!(M, matched_pos, idx, x, λ)
    num_kernels = size(M.kernels)[1]
    for n in 1:num_kernels
        K = M.kernels[n, idx]
        if n == matched_pos
            μ = K.μ + λ*(x - K.μ)
            σ² = max((1-λ)*K.σ² + λ*(x - K.μ)^2, MIN_VARIANCE)
            ω = (1-λ) * K.ω + λ
            M.kernels[n, idx] =  Gaussian(μ, σ², ω)
        else
            M.kernels[n, idx] = (1-λ) * K
        end
    end
    return move_kernel_into_order!(M, matched_pos, idx)
end

function move_kernel_into_order!(M, k, idx)
    n = k
    for n in k-1:-1:1
        if M.kernels[n, idx].ω < M.kernels[n+1, idx].ω
            M.kernels[n, idx], M.kernels[n+1, idx] = M.kernels[n+1, idx], M.kernels[n, idx]
        else
            n = n+1
            break
        end
    end
    normalize!(M, idx)
    return n
end

function replace_kernel_at!(k, idx, M, kernel::Gaussian)
    M.kernels[k, idx] = kernel
    return move_kernel_into_order!(M, k, idx)
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Raibert Policy
include("policy/2D.jl")

# ## Simulation
s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
model = s.model
env = s.env

# ## Reference Trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_forward.jld2"),
    load_type = :joint_traj))
H = ref_traj.H
h = ref_traj.h
N_sample = 5
h_sim = h / N_sample
H_sim = 1000# 100*H*N_sample

# ## Policy
v0 = 0.2
Tstance = 0.13 # measure using hop-in-place gait
Tflight = 0.62 # measure using hop-in-place gait
p = raibert_policy(s.model, v0=v0, Tstance=Tstance, Tflight=Tflight, h=h);

# ## Initial conditions
off0 = ContactImplicitMPC.SVector{model.nq}([0.0, 0.5, 0.0, 0.0])
off1 = ContactImplicitMPC.SVector{model.nq}([0.0 * v0 * h_sim, 0.5, 0.0, 0.0])
q_ref = ContactImplicitMPC.SVector{model.nq}([0.0, 0.5, 0.0, 0.5])
q0_sim = copy(q_ref) + off0
q1_sim = copy(q_ref) + off1

# ## Simulator
sim = simulator(s, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = InteriorPointOptions(
      γ_reg = 0.0,
      undercut = Inf,
      r_tol = 1.0e-8,
      κ_tol = 1.0e-8,),
    sim_opts = SimulatorOptions(warmstart = true));

# ## Simulation
@time status = simulate!(sim, verbose = true)

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
anim = visualize_robot!(vis, model, sim.traj, sample=5);


text:
@testset "read-only maps" begin
    mktempdir() do tmpdir_ro
        open(joinpath(tmpdir_ro, "myroinputfile"), "w") do io
            println(io, "Hello RO world!")
        end
        x1 = read(joinpath(tmpdir_ro, "myroinputfile"), String)
        @test x1 == "Hello RO world!\n"

        config = DockerConfig(;
            base_image = "julia:latest",
            read_only_maps = Dict(
                "/home/myuser/workdir_ro" => tmpdir_ro,
            ),
        )

        Utils.make_world_readable_recursively(tmpdir_ro)

        with_container() do container
            code = """
            y1 = read("/home/myuser/workdir_ro/myroinputfile", String)
            if y1 != "Hello RO world!\n"
                throw(ErrorException("File contents did not match"))
            end
            mktempdir() do tmpdir
                open(joinpath(tmpdir, "sometemporaryfile"), "w") do io
                    println(io, "I wrote this file into a temporary directory.")
                end
            end
            """
            @test success(container, config, `julia -e $(code)`)
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkWatcherPropertiesFormat <: SwaggerModel
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function NetworkWatcherPropertiesFormat(;provisioningState=nothing)
        o = new()
        validate_property(NetworkWatcherPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type NetworkWatcherPropertiesFormat

const _property_map_NetworkWatcherPropertiesFormat = Dict{Symbol,Symbol}(Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_NetworkWatcherPropertiesFormat = Dict{Symbol,String}(Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ NetworkWatcherPropertiesFormat }) = collect(keys(_property_map_NetworkWatcherPropertiesFormat))
Swagger.property_type(::Type{ NetworkWatcherPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_NetworkWatcherPropertiesFormat[name]))}
Swagger.field_name(::Type{ NetworkWatcherPropertiesFormat }, property_name::Symbol) =  _property_map_NetworkWatcherPropertiesFormat[property_name]

const _allowed_NetworkWatcherPropertiesFormat_provisioningState = ["Succeeded", "Updating", "Deleting", "Failed"]

function check_required(o::NetworkWatcherPropertiesFormat)
    true
end

function validate_property(::Type{ NetworkWatcherPropertiesFormat }, name::Symbol, val)
    if name === Symbol("provisioningState")
        Swagger.validate_param(name, "NetworkWatcherPropertiesFormat", :enum, val, _allowed_NetworkWatcherPropertiesFormat_provisioningState)
    end
end


text:
@scalar_rule(SpecialFunctions.lgamma(x), SpecialFunctions.digamma(x))
@scalar_rule(SpecialFunctions.erf(x), (2 / sqrt(π)) * exp(-x * x))
@scalar_rule(SpecialFunctions.erfc(x), -(2 / sqrt(π)) * exp(-x * x))
@scalar_rule(SpecialFunctions.erfi(x), (2 / sqrt(π)) * exp(x * x))
@scalar_rule(SpecialFunctions.digamma(x), SpecialFunctions.trigamma(x))
@scalar_rule(SpecialFunctions.trigamma(x), SpecialFunctions.polygamma(2, x))
@scalar_rule(SpecialFunctions.airyai(x), SpecialFunctions.airyaiprime(x))
@scalar_rule(SpecialFunctions.airyaiprime(x), x * SpecialFunctions.airyai(x))
@scalar_rule(SpecialFunctions.airybi(x), SpecialFunctions.airybiprime(x))
@scalar_rule(SpecialFunctions.airybiprime(x), x * SpecialFunctions.airybi(x))
@scalar_rule(SpecialFunctions.besselj0(x), -SpecialFunctions.besselj1(x))
@scalar_rule(SpecialFunctions.bessely0(x), -SpecialFunctions.bessely1(x))
@scalar_rule(SpecialFunctions.invdigamma(x), inv(SpecialFunctions.trigamma(SpecialFunctions.invdigamma(x))))
@scalar_rule(SpecialFunctions.besselj1(x), (SpecialFunctions.besselj0(x) - SpecialFunctions.besselj(2, x)) / 2)
@scalar_rule(SpecialFunctions.bessely1(x), (SpecialFunctions.bessely0(x) - SpecialFunctions.bessely(2, x)) / 2)
@scalar_rule(SpecialFunctions.gamma(x), Ω * SpecialFunctions.digamma(x))
@scalar_rule(SpecialFunctions.erfinv(x), (sqrt(π) / 2) * exp(Ω^2))
@scalar_rule(SpecialFunctions.erfcinv(x), -(sqrt(π) / 2) * exp(Ω^2))
@scalar_rule(SpecialFunctions.erfcx(x), (2 * x * Ω) - (2 / sqrt(π)))
@scalar_rule(SpecialFunctions.dawson(x), 1 - (2 * x * Ω))


text:
@testset "aggregate" begin
for i in 1:10
    testimage = rand(1:10.0, 10, 10)
    testmask = rand(0:1, 10, 10)
    @test aggregate(testimage, testmask) >= 0
end
end

@testset "aggregate_timeseries" begin
for i in 1:10
    testcube = rand(1:10.0, 10, 10, 10)
    testmask = rand(0:1, 10, 10)
    @test length(aggregate_timeseries(testcube, testmask)) == size(testcube)[3]
end
end


text:
# Copyright © 2021 Alexander L. Hayes

# Run from the base of the repository:
#   julia --project=. docs/make.jl
#
# This is partially to make sure the project README.md is copied and turned
# into the "Home" page.

push!(LOAD_PATH,"../src/")

cp("README.md", "docs/src/index.md", force=true)

using Documenter, RelationalDatasets

makedocs(
    sitename="RelationalDatasets.jl",
    pages = [
        "Home" => "index.md",
        "API Docs" => "api.md",
        "Types" => "types.md",
    ]
)

deploydocs(
    repo = "github.com/srlearn/RelationalDatasets.jl.git",
    devbranch = "main",
)


text:
### type definition ###

"""
A container with information about a series of patched-conic trajectories for a ballistic transfer from a starting orbit
to an ending orbit, where the two orbits are in the same solar system.
"""
struct Transfer
    startorbit::Orbit
    endorbit::Orbit
    departure_time::Float64
    arrival_time::Float64
    patch_positions::Vector{SVector{3,Float64}}
    transfer_orbits::Vector{Orbit}
    ejection_orbits::Vector{Orbit}
    insertion_orbits::Vector{Orbit}
    burns::Vector{Burn}
    Δv::Float64
    planechange::Bool
    match_sorb_Mo::Bool
    match_eorb_Mo::Bool
end

### alternate constructors ###

function Transfer(startorb::Orbit, endorb::Orbit, starttime, endtime,
                  patch_posns::Union{Nothing,AbstractVector{<:AbstractVector{<:Real}}} = nothing;
                  transferpath = path_to_body(startorb.primary, endorb.primary),
                  commonparent = closest_common_parent(startorb, endorb),
                  torb_fun::Union{typeof(p_lambert), typeof(plane_change_p_lambert)} = p_lambert,
                  dorb_fun::Union{typeof(departure_orbit), typeof(fast_departure_orbit)} = departure_orbit,
                  aorb_fun::Union{typeof(arrival_orbit), typeof(fast_arrival_orbit)} = arrival_orbit,
                  match_sorb_Mo=true, match_eorb_Mo=true,)

    tidx = findfirst(x->x==commonparent, transferpath)
    n_ejc = tidx - 1
    n_ins = length(transferpath) - tidx
    if isnothing(patch_posns)
        patch_posns = fill(@SVector(zeros(3)), length(transferpath)-1)
    end

    # transfer
    sorb = tidx == 1 ? startorb : transferpath[tidx-1].orbit
    eorb = tidx == length(transferpath) ? endorb : transferpath[tidx+1].orbit
    start_pos = n_ejc > 0 ? patch_posns[tidx-1] : @SVector(zeros(3))
    end_pos = n_ins > 0 ? patch_posns[tidx] : @SVector(zeros(3))
    torbs = torb_fun(sorb, eorb, starttime, endtime, start_pos, end_pos)
    if typeof(torbs)<:Orbit
        torbs = Orbit[torbs]
    else
        torbs = Orbit[torbs...]
    end
    v̄dep = time_orbital_velocity(starttime,torbs[1])   - time_orbital_velocity(starttime,sorb)
    v̄arr = time_orbital_velocity(endtime,  torbs[end]) - time_orbital_velocity(endtime,  eorb)

    # backward to start
    dorbs = Orbit[]
    stime = starttime
    didx = tidx-1
    while didx > 0
        sorb = didx == 1 ? startorb : transferpath[didx-1].orbit
        if didx==1 && !match_sorb_Mo
            dorb = dorb_fun(startorb, v̄dep, stime; match_pkorb_Mo=false)[1]
            r̄ₒ = time_orbital_position(dorb.epoch, dorb)
            θₚₖ = angle_in_plane(r̄ₒ, startorb)
            startorb = Orbit(startorb.a, startorb.e, startorb.i, startorb.Ω, startorb.ω,
                             true_to_mean(θₚₖ, startorb),
                             dorb.epoch, startorb.primary,
            )
        else
            dorb = dorb_fun(sorb, v̄dep, stime)[1]
        end
        if didx>1
            if norm(patch_posns[didx-1])>0
                depoch = dorb.epoch
                dorb = departarrive_true_anomaly(time_orbital_position(dorb.epoch, dorb) + patch_posns[didx-1],
                                                 dorb.primary,
                                                 v̄dep,
                                                 stime,
                                                 1)[1]
                dorb = Orbit(dorb.a, dorb.e, dorb.i, dorb.Ω, dorb.ω, dorb.Mo, depoch, dorb.primary)
            end
        end
        push!(dorbs, dorb)
        stime = dorb.epoch
        v̄dep = time_orbital_velocity(stime,dorb) - time_orbital_velocity(stime,sorb)
        didx = didx - 1
    end

    # forward to end
    aorbs = Orbit[]
    etime = endtime
    aidx = tidx+1
    while aidx < length(transferpath)+1
        eorb = aidx == length(transferpath) ? endorb : transferpath[aidx+1].orbit
        if aidx==length(transferpath) && !match_eorb_Mo
            aorb = aorb_fun(endorb, v̄arr, etime; match_pkorb_Mo=false)[1]
            r̄ₒ = time_orbital_position(aorb.epoch, aorb)
            θₚₖ = angle_in_plane(r̄ₒ, endorb)
            endorb = Orbit(endorb.a, endorb.e, endorb.i, endorb.Ω, endorb.ω,
                           true_to_mean(θₚₖ, endorb),
                           aorb.epoch, endorb.primary,
            )
        else
            aorb = aorb_fun(eorb, v̄arr, etime)[1]
        end
        if aidx<=length(patch_posns)
            if norm(patch_posns[aidx])>0
                aepoch = aorb.epoch
                aorb = departarrive_true_anomaly(time_orbital_position(aorb.epoch, aorb) + patch_posns[aidx],
                                                 aorb.primary,
                                                 v̄arr,
                                                 etime,
                                                 -1)[1]
                aorb = Orbit(aorb.a, aorb.e, aorb.i, aorb.Ω, aorb.ω, aorb.Mo, aepoch, aorb.primary)
            end
        end
        push!(aorbs, aorb)
        etime = aorb.epoch
        v̄arr = time_orbital_velocity(etime,aorb) - time_orbital_velocity(etime,eorb)
        aidx = aidx + 1
    end

    # get burns from beggining and end
    dep_orb, dep_time = !isempty(dorbs) ? (dorbs[end], dorbs[end].epoch) : (torbs[1], starttime)
    dΔv̄ = time_orbital_velocity(dep_orb.epoch, dep_orb) - time_orbital_velocity(dep_orb.epoch, startorb)
    departure_burn = Burn(startorb, dΔv̄, dep_orb.epoch)

    arr_orb, arr_time = !isempty(aorbs) ? (aorbs[end], aorbs[end].epoch) : (torbs[end], endtime)
    aΔv̄ = time_orbital_velocity(arr_time, endorb) - time_orbital_velocity(arr_time, arr_orb)
    arrival_burn = Burn(endorb, aΔv̄, arr_time)

    burns = [departure_burn]
    for i=1:length(torbs)-1
        Δv̄ = time_orbital_velocity(torbs[i+1].epoch, torbs[i+1]) - time_orbital_velocity(torbs[i+1].epoch, torbs[i])
        push!(burns, Burn(torbs[i], Δv̄, torbs[i+1].epoch))
    end
    push!(burns, arrival_burn)

    Δv = sum([norm(b.Δv̄) for b in burns])

    # construct Transfer
    return Transfer(startorb,
                    endorb,
                    starttime,
                    endtime,
                    patch_posns,
                    torbs,
                    reverse(dorbs),
                    aorbs,
                    burns,
                    Δv,
                    torb_fun==plane_change_p_lambert,
                    match_sorb_Mo,
                    match_eorb_Mo)
end

Transfer(tfer::Transfer, patch_posns::Union{Nothing,AbstractVector{<:AbstractVector{<:Real}}}=nothing; kwargs...
    ) = Transfer(tfer.startorbit, tfer.endorbit, tfer.departure_time, tfer.arrival_time, patch_posns;
                 torb_fun = tfer.planechange ? plane_change_p_lambert : p_lambert,
                 match_sorb_Mo=tfer.match_sorb_Mo, match_eorb_Mo=tfer.match_eorb_Mo,
                 # match_dorb_plane=tfer.match_dorb_plane, match_aorb_plane=tfer.match_aorb_plane,
                 kwargs...)

fastTransfer(startorb::Orbit, endorb::Orbit, starttime, endtime, patch_posns::Union{Nothing,AbstractVector{<:AbstractVector{<:Real}}} = nothing; kwargs...
    ) = Transfer(startorb, endorb, starttime, endtime, patch_posns; dorb_fun=fast_departure_orbit, aorb_fun=fast_arrival_orbit, kwargs...)

### helper methods ###

# ensure continuous time across SoI patches

function patch_time_errors(tfer::Transfer)
    up_time_mismatches   = Float64[]
    down_time_mismatches = Float64[]
    for (i,dorb) in enumerate(tfer.ejection_orbits)
        if i==length(tfer.ejection_orbits)
            t1 = tfer.departure_time
        else
            t1 = tfer.ejection_orbits[i+1].epoch
        end
        t2 = ejection_time(dorb)
        push!(up_time_mismatches, t2-t1)
    end
    for (j,aorb) in enumerate(tfer.insertion_orbits)
        if j==1
            t1 = tfer.arrival_time
        else
            t1 = tfer.insertion_orbits[j-1].epoch
        end
        t2 = insertion_time(aorb)
        push!(down_time_mismatches, t2-t1)
    end
    return up_time_mismatches, down_time_mismatches
end

function transfer_patch_time_errors(tfer::Transfer)
    up_time_mismatches, down_time_mismatches = patch_time_errors(tfer)
    start_err = isempty(up_time_mismatches)   ? 0. : up_time_mismatches[end]
    end_err =   isempty(down_time_mismatches) ? 0. : down_time_mismatches[1]
    return up_time_mismatches[end], down_time_mismatches[1]
end

"""
    tfer = match_transfer_patch_times(startorb, endorb, stime, etime, patch_posns; atol=1e-6, maxit=100)
    tfer = match_transfer_patch_times(tfer; kwargs...)

Attempts to optimize transfer start and end times to ensure that patch times are continuous
on either side of a patch, and returns the optimized Transfer.
"""
function match_transfer_patch_times(startorb, endorb, stime, etime, patch_posns; atol=1e-3, maxit=100, kwargs...)
    start_err, end_err = transfer_patch_time_errors(Transfer(startorb, endorb, stime, etime, patch_posns; kwargs...))
    err = sum(abs, [start_err, end_err])
    it=0
    # naive iteration
    while err>atol && it<maxit
        it+=1
        stime = stime + start_err
        etime = etime + end_err

        start_err, end_err = transfer_patch_time_errors(Transfer(startorb, endorb, stime, etime, patch_posns; kwargs...))
        err = sum(abs, [start_err, end_err])
    end
    # # attempt at optimization
    # function err_fun(se, startorb, endorb, patch_posns; kwargs...)
    #     start_err, end_err = transfer_patch_time_errors(Transfer(startorb, endorb, se[1], se[2], patch_posns; kwargs...))
    #     return sum(abs, [start_err, end_err])
    # end
    # obj_fun(se) = err_fun(se, startorb, endorb, patch_posns; kwargs...)
    # res = optimize(obj_fun, [stime, etime], SimulatedAnnealing())
    # @show res.minimum
    # stime, etime = res.minimizer
    return Transfer(startorb, endorb, stime, etime, patch_posns; kwargs...)
end
match_transfer_patch_times(tfer::Transfer; kwargs...
    ) = match_transfer_patch_times(tfer.startorbit, tfer.endorbit, tfer.departure_time, tfer.arrival_time, tfer.patch_positions;
                                   torb_fun = tfer.planechange ? plane_change_p_lambert : p_lambert,
                                   match_sorb_Mo=tfer.match_sorb_Mo, match_eorb_Mo=tfer.match_eorb_Mo,
                                   # match_dorb_plane=tfer.match_dorb_plane, match_aorb_plane=tfer.match_aorb_plane,
                                   kwargs...)

# ensure continuous position across SoI patches

function new_patch_positions(tfer::Transfer, coef=1)
    if tfer.startorbit.primary == tfer.endorbit.primary
       return nothing
    end
    if coef == 1
        return vcat(SVector{3,Float64}[ejection_position(dorb) for dorb in tfer.ejection_orbits],
                    SVector{3,Float64}[insertion_position(aorb) for aorb in tfer.insertion_orbits])
    else
        return vcat(SVector{3,Float64}[normalize(coef*ejection_position(dorb)+(1-coef)*tfer.patch_positions[i])*norm(ejection_position(dorb)) for (i,dorb) in enumerate(tfer.ejection_orbits)],
                    SVector{3,Float64}[normalize(coef*insertion_position(aorb)+(1-coef)*tfer.patch_positions[length(tfer.ejection_orbits)+j])*norm(insertion_position(aorb)) for (j,aorb) in enumerate(tfer.insertion_orbits)])
    end
end

function patch_position_errors(tfer::Transfer)
    if tfer.startorbit.primary == tfer.endorbit.primary
        return 0.
     end
    return tfer.patch_positions - new_patch_positions(tfer)
end

"""
    tfer = match_transfer_patch_positions(tfer; atol=1., maxit=100)

Attempts to optimize transfer start and end positions to ensure that patch times and positions are continuous
on either side of a patch, and returns the optimized Transfer.
"""
function match_patch_positions(tfer::Transfer; atol=1., maxit=100, kwargs...)
    err = sum(norm, patch_position_errors(tfer))
    best_err = err
    best_tfer = tfer
    coef = 1
    it = 0
    while best_err > atol && it<maxit
        it += 1
        tfer = match_transfer_patch_times(tfer; kwargs...)
        tfer = Transfer(tfer, new_patch_positions(tfer, coef); kwargs...)
        err = sum(norm, patch_position_errors(tfer))
        if err < best_err
            best_err = err
            best_tfer = tfer
        end
    end
    return best_tfer
end

# align start orbit with departure orbit plane
function match_departure_inclination(tfer::Transfer)
    # align start orbit with departure orbit plane
    if isempty(tfer.ejection_orbits)
        return tfer
    end
    rot = align_vectors(tfer.startorbit.basis[:,3], tfer.ejection_orbits[1].basis[:,3])
    r̄ₒ, v̄ₒ = time_state_vector(tfer.ejection_orbits[1].epoch, tfer.startorbit)
    v̄ₒ = rot * v̄ₒ
    startorb = Orbit(tfer.ejection_orbits[1].epoch, r̄ₒ, v̄ₒ, tfer.startorbit.primary, tfer.startorbit.epoch)

    dΔv̄ = time_orbital_velocity(tfer.ejection_orbits[1].epoch, tfer.ejection_orbits[1]) - time_orbital_velocity(tfer.ejection_orbits[1].epoch, startorb)
    departure_burn = Burn(startorb, dΔv̄, tfer.ejection_orbits[1].epoch)
    burns = [departure_burn, tfer.burns[2:end]...]

    Δv = sum([norm(b.Δv̄) for b in burns])

    return Transfer(startorb, tfer.endorbit, tfer.departure_time, tfer.arrival_time,
                    tfer.patch_positions, tfer.transfer_orbits, tfer.ejection_orbits, tfer.insertion_orbits,
                    burns, Δv,
                    tfer.planechange, tfer.match_sorb_Mo, tfer.match_eorb_Mo)
end

# align end orbit with arrival orbit plane
function match_arrival_inclination(tfer::Transfer)
    # align start orbit with departure orbit plane
    if isempty(tfer.insertion_orbits)
        return tfer
    end
    rot = align_vectors(tfer.endorbit.basis[:,3], tfer.insertion_orbits[end].basis[:,3])
    r̄ₒ, v̄ₒ = time_state_vector(tfer.insertion_orbits[end].epoch, tfer.endorbit)
    v̄ₒ = rot * v̄ₒ
    endorb = Orbit(tfer.insertion_orbits[end].epoch, r̄ₒ, v̄ₒ, tfer.endorbit.primary, tfer.endorbit.epoch)

    aΔv̄ = time_orbital_velocity(tfer.insertion_orbits[end].epoch, endorb) - time_orbital_velocity(tfer.insertion_orbits[end].epoch, tfer.insertion_orbits[end])
    arrival_burn = Burn(tfer.insertion_orbits[end], aΔv̄, tfer.insertion_orbits[end].epoch)
    burns = [tfer.burns[1:end-1]..., arrival_burn]

    Δv = sum([norm(b.Δv̄) for b in burns])

    return Transfer(tfer.startorbit, endorb, tfer.departure_time, tfer.arrival_time,
                    tfer.patch_positions, tfer.transfer_orbits, tfer.ejection_orbits, tfer.insertion_orbits,
                    burns, Δv,
                    tfer.planechange, tfer.match_sorb_Mo, tfer.match_eorb_Mo)
end

## Get additional inforation about transfers that won't be stored as part of the type ##

function phase_angle(tfer::Transfer)
    sorb = isempty(tfer.ejection_orbits)  ? tfer.startorbit : tfer.ejection_orbits[end].primary.orbit
    eorb = isempty(tfer.insertion_orbits) ? tfer.endorbit   : tfer.insertion_orbits[1].primary.orbit
    r̄ₛ = time_orbital_position(tfer.departure_time, sorb)
    r̄ₑ = time_orbital_position(tfer.departure_time, eorb)
    return wrap_angle(angle_in_plane(r̄ₑ, r̄ₛ, sorb.basis), -π)
end

function ejection_angle(tfer::Transfer)
    if isempty(tfer.ejection_orbits)
        return nothing
    end
    v̂ₚ = normalize(time_orbital_velocity(tfer.ejection_orbits[1].epoch, tfer.startorbit.primary.orbit))
    r̄ₒ = time_orbital_position(tfer.ejection_orbits[1].epoch, tfer.ejection_orbits[1])
    wrap_angle(angle_in_plane(r̄ₒ, v̂ₚ, tfer.startorbit.basis), -π)
end

function transfer_angle(tfer::Transfer)
    r̄ₛ = time_orbital_position(tfer.departure_time, tfer.transfer_orbits[1])
    r̄ₑ = time_orbital_position(tfer.arrival_time, tfer.transfer_orbits[end])
    return wrap_angle(angle_in_plane(r̄ₑ, r̄ₛ, basis_rotation(r̄ₛ, r̄ₑ)))
end

### descriptive display ###

Base.show(io::IO, tfer::Transfer) = println(io,
    "Transfer from ",
    "orbit around $(tfer.startorbit.primary.name) to ",
    "orbit around $(tfer.endorbit.primary.name)\n",
    "  departure time:\t$(tfer.departure_time) s\n",
    "  arrival time:\t\t$(norm(tfer.arrival_time)) s\n",
    "  Δv:\t\t\t$(tfer.Δv) m/s",
)


text:
# Converts an FFT into a gammatonegram. Based on the gammatonegram/fft2gammatonemx MATLAB
# functions by Dan Ellis.
# D. P. W. Ellis (2009). "Gammatone-like spectrograms", web resource.
# http://www.ee.columbia.edu/~dpwe/resources/matlab/gammatonegram/

type Gammatonegram{T, F<:Real} <: DSP.Periodograms.TFR{T}
	amplitude::Matrix{T}
	frequencies::Vector{F}
	time::FloatRange{Float64}
end

function gammatonegram(x,sr::Integer,twin::Real,thop::Real,N::Integer,fmin,fmax,width)
	nfft = int(2^(ceil(log(2*twin*sr)/log(2))))
    nhop = iround(thop*sr)
    nwin = iround(twin*sr)
    (W,F) = fft2gammatonemx(nfft, sr, N, width, fmin, fmax)
    # perform FFT and weighting in amplitude domain
	S = stft(x, nwin, nwin-nhop; nfft=nfft, fs=sr, window=hanning)
	Y = 1/nfft*W*abs(S)
	Gammatonegram(Y, vec(F), ((0:size(Y,2)-1)*(nhop)+nwin/2)/sr)
end

function fft2gammatonemx(nfft::Integer, sr::Integer, N::Integer, width, fmin, fmax)
	# Allocating matrix to store the filterbank weights
	W = zeros(N, int(nfft/2+1));
	# Gammatone filterbank constants
	EarQ = 9.26449
	minBW = 24.7
	order = 1

	cfreqs = -(EarQ*minBW) + exp((1:N)'*(-log(fmax + EarQ*minBW) + log(fmin + EarQ*minBW))/N) * (fmax + EarQ*minBW);
	cfreqs = fliplr(cfreqs)
	GTord = 4
	ucirc = exp(im*2*pi*[0:(nfft/2)]/nfft)
	for k=1:N
	    cf = cfreqs[k];
	    ERB = width*((cf/EarQ).^order + minBW^order).^(1/order);
	    B = 1.019*2*pi*ERB;
	    r = exp(-B/sr);
	    theta = 2*pi*cf/sr;
	    pole = r*exp(im*theta);
	    T = 1/sr;

		A11 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3+2^1.5)*T*sin(2*
	                                                      cf*pi*T)./exp(B*T))/2;
	    A12 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3+2^1.5)*T*sin(2*
	                                                      cf*pi*T)./exp(B*T))/2;
	    A13 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3-2^1.5)*T*sin(2*
	                                                      cf*pi*T)./exp(B*T))/2;
	    A14 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3-2^1.5)*T*sin(2*
	                                                      cf*pi*T)./exp(B*T))/2;
	    zros = -[A11 A12 A13 A14]/T;

	    gain =  abs((-2*exp(4*im*cf*pi*T)*T +
	                2*exp(-(B*T) + 2*im*cf*pi*T).*T.*
	                (cos(2*cf*pi*T) - sqrt(3 - 2^(3/2))*
	                 sin(2*cf*pi*T))) .*
	               (-2*exp(4*im*cf*pi*T)*T +
	                2*exp(-(B*T) + 2*im*cf*pi*T).*T.*
	                (cos(2*cf*pi*T) + sqrt(3 - 2^(3/2)) *
	                 sin(2*cf*pi*T))).*
	               (-2*exp(4*im*cf*pi*T)*T +
	                2*exp(-(B*T) + 2*im*cf*pi*T).*T.*
	                (cos(2*cf*pi*T) -
	                 sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) .*
	               (-2*exp(4*im*cf*pi*T)*T + 2*exp(-(B*T) + 2*im*cf*pi*T).*T.*
	                (cos(2*cf*pi*T) + sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) ./
	               (-2 ./ exp(2*B*T) - 2*exp(4*im*cf*pi*T) +
	                2*(1 + exp(4*im*cf*pi*T))./exp(B*T)).^4);
	    W[k,:] = ((T^4)/gain) *
	        abs(ucirc-zros[1]).*abs(ucirc-zros[2]) .*
	        abs(ucirc-zros[3]).*abs(ucirc-zros[4]) .*
	        (abs((pole-ucirc).*(pole'-ucirc)).^-GTord);
	end
	W, cfreqs
end


text:
# ---
# title: 218. The Skyline Problem
# id: problem218
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Divide and Conquer, Heap, Binary Indexed Tree, Segment Tree, Line Sweep
# link: <https://leetcode.com/problems/the-skyline-problem/description/>
# hidden: true
# ---
#
# A city's **skyline** is the outer contour of the silhouette formed by all the
# buildings in that city when viewed from a distance. Given the locations and
# heights of all the buildings, return _the **skyline** formed by these
# buildings collectively_.
#
# The geometric information of each building is given in the array `buildings`
# where `buildings[i] = [lefti, righti, heighti]`:
#
#   * `lefti` is the x coordinate of the left edge of the `ith` building.
#   * `righti` is the x coordinate of the right edge of the `ith` building.
#   * `heighti` is the height of the `ith` building.
#
# You may assume all buildings are perfect rectangles grounded on an absolutely
# flat surface at height `0`.
#
# The **skyline** should be represented as a list of "key points" **sorted by
# their x-coordinate** in the form `[[x1,y1],[x2,y2],...]`. Each key point is
# the left endpoint of some horizontal segment in the skyline except the last
# point in the list, which always has a y-coordinate  `0` and is used to mark
# the skyline's termination where the rightmost building ends. Any ground
# between the leftmost and rightmost buildings should be part of the skyline's
# contour.
#
# **Note:** There must be no consecutive horizontal lines of equal height in the
# output skyline. For instance, `[...,[2 3],[4 5],[7 5],[11 5],[12 7],...]` is
# not acceptable; the three lines of height 5 should be merged into one in the
# final output as such: `[...,[2 3],[4 5],[12 7],...]`
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/12/01/merged.jpg)
#
#
#
#     Input: buildings = [[2,9,10],[3,7,15],[5,12,12],[15,20,10],[19,24,8]]
#     Output: [[2,10],[3,15],[7,12],[12,0],[15,10],[20,8],[24,0]]
#     Explanation:
#     Figure A shows the buildings of the input.
#     Figure B shows the skyline formed by those buildings. The red points in figure B represent the key points in the output list.
#
#
# **Example 2:**
#
#
#
#     Input: buildings = [[0,2,3],[2,5,3]]
#     Output: [[0,3],[5,0]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= buildings.length <= 104`
#   * `0 <= lefti < righti <= 231 - 1`
#   * `1 <= heighti <= 231 - 1`
#   * `buildings` is sorted by `lefti` in non-decreasing order.
#
#
## @lc code=start
using LeetCode

function get_skyline(buildings::Vector{Vector{Int}})
    all = Tuple{Int, Int}[]
    res = Tuple{Int, Int}[]
    @inbounds for b in buildings
        push!(all, (b[1], -b[3]))
        push!(all, (b[2], b[3]))
    end
    sort!(all)

    heights = SortedMultiDict{Int, Nothing}()
    push!(heights, 0 => nothing)
    lst = (0, 0)
    for p in all
        p[2] < 0 ? push!(heights, -p[2] => nothing) : delete!((heights, searchsortedfirst(heights, p[2])))
        mh = last(heights).first
        if lst[2] != mh
            lst = (p[1], mh)
            push!(res, lst)
        end
    end
    return res
end
## @lc code=end


text:
using Lynx
using Luxor

Lynx.init("Events", 800, 600)

mouse = Point(0, 0)

function update(dt)
    background("black")
    sethue("#1fe")
    circle(mouse, 50, :fill)
end

onkeypress(@window) do event
    key = keyname(event)
    @info "Key $key pressed!" code=event.keyval
end

onevent(:destroy, @window) do win
    @info "Window has been closed!"
end

onmousemove(@canvas) do event
    global mouse = Point(event.x, event.y)
    @info "Mouse moved to $(mouse)"
end

onmousedown(@canvas) do event
    x, y = event.x, event.y
    btn = event.button
    @info "Mouse pressed at ($x, $y) with button $btn"
end

onmouseup(@canvas) do event
    x, y = event.x, event.y
    btn = event.button
    @info "Mouse released at ($x, $y) with button $btn"
end

onmousedrag(@canvas) do event
    x, y = event.x, event.y
    btn = event.button
    @info "Mouse drag at ($x, $y) with button $btn"
end

run!(update)


text:
"""
Unsupported argument types
"""
const RESERVED_ARG_TYPES = ["va_list"]

"""
    wrap!(ctx::AbstractContext, cursor::CLFunctionDecl)
Subroutine for handling function declarations. Note that VarArg functions are not supported.
"""
function wrap!(ctx::AbstractContext, cursor::CLFunctionDecl)
    func_type = type(cursor)
    if kind(func_type) == CXType_FunctionNoProto
        @warn "No Prototype for $cursor - assuming no arguments"
    elseif isvariadic(func_type)
        @warn "Skipping VarArg Function $cursor"
        return ctx
    end

    func_name = isempty(ctx.force_name) ? Symbol(spelling(cursor)) : ctx.force_name
    ret_type = clang2julia(return_type(cursor))
    args = function_args(cursor)
    arg_types = [argtype(func_type, i) for i in 0:length(args)-1]
    arg_reps = clang2julia.(arg_types)
    for (i, arg) in enumerate(arg_reps)
        # constant array argument should be converted to Ptr
        # e.g. double f[3] => Ptr{Cdouble} instead of NTuple{3, Cdouble}
        if Meta.isexpr(arg, :curly) && first(arg.args) == :NTuple
            arg_reps[i] = Expr(:curly, :Ptr, last(arg.args))
        end
    end

    # check whether any argument types are blocked
    for t in arg_types
        if spelling(t) in RESERVED_ARG_TYPES
            @warn "Skipping $(name(cursor)) due to unsupported argument: $(spelling(t))"
            return ctx
        end
    end

    # handle unnamed args and convert names to symbols
    arg_count = 0
    arg_names = map(args) do x
                    n = name_safe(name(x))
                    s = !isempty(n) ? n : "arg"*string(arg_count+=1)
                    Symbol(s)
                end

    isstrict = get(ctx.options, "is_function_strictly_typed", true)
    signature = isstrict ? efunsig(func_name, arg_names, arg_reps) : Expr(:call, func_name, arg_names...)

    ctx.libname == "libxxx" && @warn "default libname: \":libxxx\" are being used, did you forget to specify `context.libname`?"
    body = eccall(func_name, Symbol(ctx.libname), ret_type, arg_names, arg_reps)

    push!(ctx.api_buffer, Expr(:function, signature, Expr(:block, body)))

    return ctx
end

function is_ptr_type_expr(@nospecialize t)
    (t === :Cstring || t === :Cwstring) && return true
    isa(t, Expr) || return false
    t = t::Expr
    t.head === :curly && t.args[1] === :Ptr
end

function efunsig(name::Symbol, args::Vector{Symbol}, types)
    x = [is_ptr_type_expr(t) ? a : Expr(:(::), a, t) for (a,t) in zip(args,types)]
    Expr(:call, name, x...)
end

function eccall(func_name::Symbol, libname::Symbol, rtype, args, types)
  :(ccall(($(QuoteNode(func_name)), $libname),
            $rtype,
            $(Expr(:tuple, types...)),
            $(args...))
    )
end

"""
    wrap!(ctx::AbstractContext, cursor::CLEnumDecl)
Subroutine for handling enum declarations.
"""
function wrap!(ctx::AbstractContext, cursor::CLEnumDecl)
    cursor_name = name(cursor)
    # handle typedef anonymous enum
    idx = ctx.children_index
    if 0 < idx < length(ctx.children)
        next_cursor = ctx.children[idx+1]
        if is_typedef_anon(cursor, next_cursor)
            cursor_name = name(next_cursor)
        end
    end
    !isempty(ctx.force_name) && (cursor_name = ctx.force_name;)
    cursor_name == "" && (@warn("Skipping unnamed EnumDecl: $cursor"); return ctx)

    enum_sym = symbol_safe(cursor_name)
    enum_type = INT_CONVERSION[clang2julia(cursor)]
    name2value = Tuple{Symbol,Int}[]
    # extract values and names
    for item_cursor in children(cursor)
        kind(item_cursor) == CXCursor_PackedAttr && (@warn("this is a `__attribute__((packed))` enum, the underlying alignment of generated structure may not be compatible with the original one in C!"); continue)
        item_name = spelling(item_cursor)
        isempty(item_name) && continue
        item_sym = symbol_safe(item_name)
        push!(name2value, (item_sym, value(item_cursor)))
    end

    expr = Expr(:macrocall, Symbol("@cenum"), nothing, Expr(:(::), enum_sym, enum_type))
    enum_pairs = Expr(:block)
    ctx.common_buffer[enum_sym] = ExprUnit(expr)
    for (name,value) in name2value
        ctx.common_buffer[name] = ctx.common_buffer[enum_sym]  ##???
        push!(enum_pairs.args, :($name = $value))
    end
    push!(expr.args, enum_pairs)

    return ctx
end

"""
    wrap!(ctx::AbstractContext, cursor::CLStructDecl)
Subroutine for handling struct declarations.
"""
function wrap!(ctx::AbstractContext, cursor::CLStructDecl)
    # make sure a empty struct is indeed an opaque struct typedef/typealias
    # cursor = canonical(cursor)  # this won't work
    cursor = type(cursor) |> canonical |> typedecl
    cursor_name = name(cursor)
    # handle typedef anonymous struct
    idx = ctx.children_index
    if 0 < idx < length(ctx.children)
        next_cursor = ctx.children[idx+1]
        if is_typedef_anon(cursor, next_cursor)
            cursor_name = name(next_cursor)
        end
    end
    !isempty(ctx.force_name) && (cursor_name = ctx.force_name;)
    cursor_name == "" && (@warn("Skipping unnamed StructDecl: $cursor"); return ctx)

    struct_sym = symbol_safe(cursor_name)
    ismutable = get(ctx.options, "is_struct_mutable", false)
    buffer = ctx.common_buffer

    # generate struct declaration
    block = Expr(:block)
    expr = Expr(:struct, ismutable, struct_sym, block)
    deps = OrderedSet{Symbol}()
    struct_fields = children(cursor)
    for (field_idx, field_cursor) in enumerate(struct_fields)
        field_name = name(field_cursor)
        field_kind = kind(field_cursor)
        if field_kind == CXCursor_StructDecl || field_kind == CXCursor_UnionDecl
            continue
        elseif field_kind == CXCursor_FirstAttr
            continue
        elseif field_kind != CXCursor_FieldDecl || field_kind == CXCursor_TypeRef
            buffer[struct_sym] = ExprUnit(Poisoned())
            @warn "Skipping struct: \"$cursor\" due to unsupported field: $field_cursor"
            return ctx
        elseif isempty(field_name)
            error("Unnamed struct member in: $cursor ... cursor: $field_cursor")
        end

        if occursin("anonymous", string(clang2julia(field_cursor)))
            idx = field_idx-1
            anonymous_record = struct_fields[idx]
            while idx != 0 && kind(anonymous_record) == CXCursor_FieldDecl
                idx -= 1
                anonymous_record = struct_fields[idx]
            end
            if idx == field_idx-1
                ctx.anonymous_counter += 1
                anon_name = "ANONYMOUS$(ctx.anonymous_counter)_"*spelling(field_cursor)
                ctx.force_name = anon_name
                wrap!(ctx, anonymous_record)
                ctx.force_name = ""
                repr = symbol_safe(anon_name)
            else
                anon_name = "ANONYMOUS$(ctx.anonymous_counter)_"*spelling(struct_fields[idx+1])
                repr = symbol_safe(anon_name)
            end
        else
            repr = clang2julia(field_cursor)
        end
        push!(block.args, Expr(:(::), symbol_safe(field_name), repr))
        push!(deps, target_type(repr))
    end

    # check for a previous forward ordering
    if !(struct_sym in keys(buffer)) || buffer[struct_sym].state == :empty
        if !isempty(struct_fields)
            buffer[struct_sym] = ExprUnit(expr, deps)
        else
            # opaque struct typedef/typealias
            buffer[struct_sym] = ExprUnit(:(const $struct_sym = Cvoid), deps, state=:empty)
        end
    end

    return ctx
end

"""
    wrap!(ctx::AbstractContext, cursor::CLUnionDecl)
Subroutine for handling union declarations.
"""
function wrap!(ctx::AbstractContext, cursor::CLUnionDecl)
    # make sure a empty union is indeed an opaque union typedef/typealias
    cursor = type(cursor) |> canonical |> typedecl
    cursor_name = name(cursor)
    # handle typedef anonymous union
    idx = ctx.children_index
    if 0 < idx < length(ctx.children)
        next_cursor = ctx.children[idx+1]
        if is_typedef_anon(cursor, next_cursor)
            cursor_name = name(next_cursor)
        end
    end
    !isempty(ctx.force_name) && (cursor_name = ctx.force_name;)
    cursor_name == "" && (@warn("Skipping unnamed UnionDecl: $cursor"); return ctx)

    union_sym = symbol_safe(cursor_name)
    ismutable = get(ctx.options, "is_struct_mutable", false)
    buffer = ctx.common_buffer

    block = Expr(:block)
    expr = Expr(:struct, ismutable, union_sym, block)
    deps = OrderedSet{Symbol}()
    # find the largest union field and declare a block of bytes to match.
    union_fields = children(cursor)
    if !isempty(union_fields)
        max_size = 0
        largest_field_idx = 0
        for i = 1:length(union_fields)
            field_cursor = union_fields[i]
            field_kind = kind(field_cursor)
            (field_kind == CXCursor_StructDecl || field_kind == CXCursor_UnionDecl) && continue
            field_kind == CXCursor_FirstAttr && continue
            field_size = typesize(type(field_cursor))
            if field_size > max_size
                max_size = field_size
                largest_field_idx = i
            end
        end
        largest_field = union_fields[largest_field_idx]
        if occursin("anonymous", string(clang2julia(largest_field)))
            idx = largest_field_idx-1
            anonymous_record = union_fields[idx]
            while idx != 0 && kind(anonymous_record) == CXCursor_FieldDecl
                idx -= 1
                anonymous_record = union_fields[idx]
            end
            ctx.anonymous_counter += 1
            anon_name = "ANONYMOUS$(ctx.anonymous_counter)_"*spelling(largest_field)
            ctx.force_name = anon_name
            wrap!(ctx, anonymous_record)
            ctx.force_name = ""
            repr = symbol_safe(anon_name)
        else
            repr = clang2julia(largest_field)
        end
        largest_field_sym = symbol_safe(spelling(largest_field))
        push!(block.args, Expr(:(::), largest_field_sym, repr))
        push!(deps, target_type(repr))
        buffer[union_sym] = ExprUnit(expr, deps)
    elseif !(union_sym in keys(buffer)) || buffer[union_sym].state == :empty
        buffer[union_sym] = ExprUnit(:(const $union_sym = Cvoid), deps, state=:empty)
    else
        @warn "Skipping union: \"$cursor\" due to unknown cases."
    end

    return ctx
end

"""
    wrap!(ctx::AbstractContext, cursor::CLTypedefDecl)
Subroutine for handling typedef declarations.
"""
function wrap!(ctx::AbstractContext, cursor::CLTypedefDecl)
    td_type = underlying_type(cursor)
    td_sym = isempty(ctx.force_name) ? Symbol(spelling(cursor)) : ctx.force_name
    buffer = ctx.common_buffer
    if kind(td_type) == CXType_Unexposed
        # TODO: which corner case will trigger this pass?
        @error "Skipping Typedef: CXType_Unexposed, $cursor, please report this on Github."
    end

    if kind(td_type) == CXType_FunctionProto
        # TODO: need to find a test case too
        if !haskey(buffer, td_sym)
            buffer[td_sym] = ExprUnit(string("# Skipping Typedef: CXType_FunctionProto ", spelling(cursor)))
        end
        return ctx
    end

    td_target = clang2julia(td_type)
    if !haskey(buffer, td_sym)
        buffer[td_sym] = ExprUnit(:(const $td_sym = $td_target), [td_target])
    end
    return ctx
end

"""
    handle_macro_exprn(tokens::TokenList, pos::Int)
For handling of #define'd constants, allows basic expressions but bails out quickly.
"""
function handle_macro_exprn(tokens::TokenList, pos::Int)
    function trans(tok)
        ops = ["+" "-" "*" "~" ">>" "<<" "/" "\\" "%" "|" "||" "^" "&" "&&" "=="]
        token_kind = kind(tok)
        (token_kind == CXToken_Literal || token_kind == CXToken_Identifier) && return 0
        token_kind == CXToken_Punctuation && tok.text ∈ ops && return 1
        return -1
    end

    # normalize literal with a size suffix
    function literally(tok)
        # note: put multi-character first, or it will break out too soon for those!
        literalsuffixes = ["ULL", "Ull", "uLL", "ull", "LLU", "LLu", "llU", "llu",
                           "LL", "ll", "UL", "Ul", "uL", "ul", "LU", "Lu", "lU", "lu",
                           "U", "u", "L", "l", "F", "f"]

        function literal_totype(literal, txt)
          literal = lowercase(literal)

          # Floats following http://en.cppreference.com/w/cpp/language/floating_literal
          float64 = occursin(".", txt) && occursin("l", literal)
          float32 = occursin("f", literal)

          if float64 || float32
            float64 && return "Float64"
            float32 && return "Float32"
          end

          # Integers following http://en.cppreference.com/w/cpp/language/integer_literal
          unsigned = occursin("u", literal)
          nbits = count(x -> x == 'l', literal) == 2 ? 64 : 32
          return "$(unsigned ? "U" : "")Int$nbits"
        end

        token_kind = kind(tok)
        txt = tok.text |> strip
        if token_kind == CXToken_Identifier || token_kind == CXToken_Punctuation
            # pass
        elseif token_kind == CXToken_Literal
            for sfx in literalsuffixes
                if endswith(txt, sfx)
                    type = literal_totype(sfx, txt)
                    txt = txt[1:end-length(sfx)]
                    txt = "$(type)($txt)"
                    break
                end
            end
        end
        return txt
    end

    # check whether identifiers and literals alternate
    # with punctuation
    exprn = ""
    pos > length(tokens) && return exprn, pos

    prev = 1 >> trans(tokens[pos])
    for lpos = pos:length(tokens)
        pos = lpos
        tok = tokens[lpos]
        state = trans(tok)
        if xor(state, prev) == 1
            prev = state
        else
            break
        end
        exprn = exprn * literally(tok)
    end
    return exprn, pos
end

# TODO: This really returns many more symbols than we want,
# Functionally, it shouldn't matter, but eventually, we
# might want something more sophisticated.
# (Check: Does this functionality already exist elsewhere?)
get_symbols(s) = Any[]
get_symbols(s::Symbol) = Any[s]
get_symbols(e::Expr) = vcat(get_symbols(e.head), get_symbols(e.args))
get_symbols(xs::Array) = reduce(vcat, [get_symbols(x) for x in xs])

"""
    wrap!(ctx::AbstractContext, cursor::CLMacroDefinition)
Subroutine for handling macro declarations.
"""
function wrap!(ctx::AbstractContext, cursor::CLMacroDefinition)
    tokens = tokenize(cursor)
    # Skip any empty definitions
    tokens.size < 2 && return ctx
    startswith(name(cursor), "_") && return ctx

    buffer = ctx.common_buffer
    pos = 1; exprn = ""
    if tokens[2].text == "("
        exprn, pos = handle_macro_exprn(tokens, 3)
        if pos != lastindex(tokens) || tokens[pos].text != ")" || exprn == ""
            mdef_str = join([c.text for c in tokens], " ")
            buffer[Symbol(mdef_str)] = ExprUnit(string("# Skipping MacroDefinition: ", replace(mdef_str, "\n"=>"\n#")))
            return ctx
        end
        exprn = "(" * exprn * ")"
    else
        exprn, pos = handle_macro_exprn(tokens, 2)
        if pos != lastindex(tokens)
            mdef_str = join([c.text for c in tokens], " ")
            buffer[Symbol(mdef_str)] = ExprUnit(string("# Skipping MacroDefinition: ", replace(mdef_str, "\n"=>"#\n")))
            return ctx
        end
    end

    # Occasionally, skipped definitions slip through
    (exprn == "" || exprn == "()") && return buffer

    use_sym = symbol_safe(tokens[1].text)

    try
        target = Meta.parse(exprn)
        e = Expr(:const, Expr(:(=), use_sym, target))
        deps = get_symbols(target)
        buffer[use_sym] = ExprUnit(e, deps)
    catch err
        # this assumes all parsing failures are due to string-parsing
        ## TODO: find a elegant way to solve this
        e = :(const $use_sym = $(exprn[2:end-1]))
        buffer[use_sym] = ExprUnit(e,[])
    end

    return ctx
end

"""
    wrap!(ctx::AbstractContext, cursor::CLTypeRef)
For now, we just skip CXCursor_TypeRef cursors.
"""
function wrap!(ctx::AbstractContext, cursor::CLTypeRef)
    @warn "Skipping CXCursor_TypeRef cursor: $cursor"
    return ctx
end

function wrap!(ctx::AbstractContext, cursor::CLCursor)
    @warn "not wrapping $(cursor)"
    return ctx
end

function wrap!(ctx::AbstractContext, cursor::Union{CLLastPreprocessing,CLMacroInstantiation})
    @debug "not wrapping $(cursor)"
    return ctx
end


text:
lines = readlines("input.txt")

edges = [Tuple(split(line, '-')) for line = lines]
nodes = unique(vcat(collect.(edges)...))
N = length(nodes)
adjacency = falses((N,N))
node2idx = n -> findfirst(nodes .== n)
idx2node = i -> nodes[i]

function isupper(s)
    return isuppercase.(collect(s)) |> all
end

function DFS(from, to, cond, state)
    state = [state..., from]
    if from == to
        return [state]
    else
        C_n = findall(adjacency[from, :] .== true)
        paths = [DFS(n, to, cond, state) for n ∈ C_n if cond(from, n, state)]
        return length(paths) > 0 ? reduce(vcat, paths) : []
    end
end

for (from, to) = edges
    ix, iy = node2idx(from), node2idx(to)
    adjacency[ix,iy] = adjacency[iy,ix] = true
end

cond = (from, to, state) -> from != node2idx("end") && (isupper(idx2node(to)) || !(to ∈ state))
paths = DFS(node2idx("start"), node2idx("end"), cond, [])
println(length(paths))


text:
######################################################################
# General purpose patterns realized through abstract types & structs
# -----
# Licensed under MIT License

export Ident
struct Ident{S} end

export Mutable
mutable struct Mutable{T}
    value::T
end

Base.getindex(mutable::Mutable) = mutable.value
Base.setindex!(mutable::Mutable{T}, value::T) where T = mutable.value = value
Base.:(==)(lhs::Mutable{T}, rhs::Mutable{T}) where T = lhs.value == rhs.value

export Unknown, unknown
struct Unknown end
const unknown = Unknown()

export Optional
mutable struct Optional{S, T}
    value::Union{T, Unknown}
    Optional{S, T}(value) where {S, T} = new(convert(T, value))
    Optional{S, T}(::Unknown) where {S, T} = new(unknown)
    Optional{S, T}(value::T) where {S, T} = new(value)
end
Optional{S, T}() where {S, T} = Optional{S, T}(unknown)
Optional{S}(value::T) where {S, T} = Optional{S, T}(value)
Optional{S}(::Unknown) where S = Optional{S, Any}(unknown)
Optional{S}() where S = Optional{S}(unknown)
Optional(S::Symbol, value::T) where T = Optional{S, T}(value)
Optional(value::T) where T = Optional{:generic, T}(value)
Optional(::Unknown) = Optional{:generic, Any}(unknown)
Optional() = Optional(unknown)

function Base.getindex(opt::Optional)
    if isunknown(opt)
        load(opt)
    end
    opt.value
end
Base.setindex!(opt::Optional, value) = opt.value = value

Base.convert(::Type{Optional{S}}, value) where S = Optional{S}(value)
Base.convert(::Type{Optional{S}}, opt::Optional) where S = (println("convert $(opt.value)"); Optional{S}(opt.value))
Base.convert(::Type{Optional{S, T}}, value) where {S, T} = Optional{S, T}(value)
Base.convert(::Type{Optional{S, T}}, opt::Optional) where {S, T} = Optional{S, T}(opt.value)

Base.convert(::Type{Optional{S} where S}, value) = Optional{:generic}(value)
Base.convert(::Type{Optional{S, T} where S}, value) where T = Optional{:generic, T}(value)
Base.convert(::Type{Optional{S} where S}, opt::Optional) = opt
Base.convert(::Type{Optional{S1, T} where S1}, opt::Optional{S2, T} where S2) where T = opt
Base.convert(::Type{Optional{S1, T} where S1}, opt::Optional{S2}) where {S2, T} = Optional{S2, T}(opt.value)

export isunknown
isunknown(opt::Optional) = opt.value === unknown

export load
"""`load(::Optional)` loads the `Optional`'s associated value from remote storage and caches it."""
load(opt::Optional) = opt.value = unknown

export Dirty, markdirty!, isdirty
mutable struct Dirty{T}
    value::T
    dirty::Bool
end
Dirty(value) = Dirty(value, false)

Base.getindex(dirty::Dirty) = dirty.value
Base.setindex!(dirty::Dirty, value) = (markdirty!(dirty); dirty.value = value)
markdirty!(dirty::Dirty) = dirty.dirty = true
isdirty(dirty::Dirty) = dirty.dirty
clear(dirty::Dirty) = dirty.dirty = false


text:
##################################################
## Qt utilities

## connect a callback to clicked event
qconnect(callback::Function, object::PyObject, signal::Symbol) = object[signal][:connect](callback)
qconnect(object::PyObject, signal::Symbol, callback::Function) = object[signal][:connect](callback)
qconnect(object::PyObject, signal::Symbol, slot::PyObject) = object[signal][:connect](slot)

qemit(object::PyObject, signal::ASCIIString) = object[:emit](QtCore[:SIGNAL](signal))
qemit(object::PyObject, signal::Symbol) = qemit(object, "$signal()")

## o[:symbol](xs...) -> qcall(o, :symbol, xs...)
qcall(object::PyObject, member::Symbol, args...) = object[member](args...)

## get Qt.Qt.enum as pyobject (http://qt-project.org/doc/qt-4.8/qt.html)
## TODO: generelize to others by passing in second argument

qt_enum(attr::Symbol) = qt_enum(string(attr))
qt_enum(attr::ASCIIString) = QtNamespace[attr] # want pyobject here

qt_enum(attr::Vector{ASCIIString}; how="|") = PyCall.pyeval(join(map(u -> "QtCore.Qt.$u", attr), " $how "))

## some converstion
## This seems unnecessary with PySide
## convert text to string ## this is python2 http://pyqt.sourceforge.net/Docs/PyQt4/gotchas.html#python-strings-qt-strings-and-unicode
### from_QString(o::PyObject)  = pyeval("str(x)", Dict([(:x, o)]))

qexec(object::PyObject) = convert(Function, object[:exec])()

##################################################
## Methods for extending classes. Hacky.
newclass_tpl = Mustache.template_from_file(Pkg.dir("PySide", "tpl", "newclass.tpl"))

## Create new class from parent
## qclass("Example", "QtGui.QWidget") ## not Qt.QWidget...
## can pass in methods as strings.
## e.g. qnew_class("Julia", "QtCore.QObject", meths=[{:meth_dfn => "evaluate = QtCore.Signal(str)"}])
## then can connect to evalute signal, qconnect(julia, :evaluate, (cmd) -> print(cmd))
function qnew_class(name::ASCIIString, parent::ASCIIString; meths=nothing)
    tmp = tempname() * ".py"
    io = open(tmp, "w")
    d = Dict{Symbol, Any}()
    d[:NewClass] = name
    d[:OldClass] = parent
    if !isa(meths, Nothing) d[:methods] = meths end
    out = Mustache.render(io, newclass_tpl, d)
    println(out)
    close(io)
    if pyversion < v"3"
        PyCall.pyeval("execfile('$tmp', globals())")
    else
        PyCall.pyeval("exec(open('$tmp', 'rb'), globals())")
    end
end

## make an instance of the new class
qnew_class_instance(name::ASCIIString) = PyCall.pyeval("$name()")
qnew_class_instance(name::ASCIIString, parent) = PyCall.pyeval("$name(x)", x=parent)

## add a method to the class.
## qnew_class("Example", "QtGui.QWidget")
## o = qnew_class_instance("Example")
## qset_method(o, :enterEvent, e -> println("enter"))
## or
## qset_method(o, :enterEvent) do e
##    println("enter")
## end
qset_method(object::PyObject, meth::Symbol, callback::Function) = object[meth] = callback
qset_method(callback::Function, object::PyObject, meth::Symbol) = qset_method(object, meth, callback)


text:
# This file was generated, do not modify it. # hide
using MLJ, Statistics, PrettyPrinting
MLJ.color_off() # hide
X, y = @load_iris;

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
@discord_object struct AuditLogEntry
    target_id::String
    changes::Vector{AuditLogChange}
    user_id::Snowflake
    id::Snowflake
    action_type::AuditLogEvent.AuditLogEventEnum
    options::AuditLogInfo
    reason::String
end


text:
@testset "wrap" begin
    @testset "basic" begin
        @system SWrap(Controller) begin
            a => 1 ~ preserve
            b(a) => 2a ~ track
            c(wrap(a)) => 2a' ~ track
        end
        s = instance(SWrap)
        @test s.b' == s.c' == 2
    end
end