text:
import Base: merge, vcat, map

###### merge ####################

function merge(ta1::TimeArray{T, N, D}, ta2::TimeArray{T, M, D}, method::Symbol=:inner;
               colnames::Vector=[], meta::Any=Void) where {T, N, M, D}

    if ta1.meta == ta2.meta && meta == Void
        meta = ta1.meta
    elseif typeof(ta1.meta) <: AbstractString && typeof(ta2.meta) <: AbstractString && meta == Void
        meta = string(ta1.meta, "_", ta2.meta)
    else
        meta = meta
    end

    if method == :inner

        idx1, idx2 = overlaps(ta1.timestamp, ta2.timestamp)
        vals = [ta1[idx1].values ta2[idx2].values]
        ta = TimeArray(ta1[idx1].timestamp, vals, [ta1.colnames; ta2.colnames], meta)

    elseif method == :left

        new_idx2, old_idx2 = overlaps(ta1.timestamp, ta2.timestamp)
        right_vals = NaN * zeros(length(ta1), length(ta2.colnames))
        right_vals[new_idx2, :]  = ta2.values[old_idx2, :]
        ta = TimeArray(ta1.timestamp, [ta1.values right_vals], [ta1.colnames; ta2.colnames], meta)

    elseif method == :right

        ta = merge(ta2, ta1, :left)
        ncol2 = length(ta2.colnames)
        vals = [ta.values[:, (ncol2+1):end] ta.values[:, 1:ncol2]]
        ta = TimeArray(ta.timestamp, vals, [ta1.colnames; ta2.colnames], meta)

    elseif method == :outer

        timestamps = sorted_unique_merge(ta1.timestamp, ta2.timestamp)
        ta = TimeArray(timestamps, zeros(length(timestamps), 0), String[], Void)
        ta = merge(ta, ta1, :left)
        ta = merge(ta, ta2, :left, meta=meta)

    else
        throw(ArgumentError(
            "merge method must be one of :inner, :left, :right, :outer"))
    end

    return setcolnames!(ta, colnames)

end

# collapse ######################

function collapse(ta::TimeArray{T, N, D}, period::Function, timestamp::Function,
                  value::Function=timestamp) where {T, N, D}

    length(ta) == 0 && return ta

    ncols = length(ta.colnames)
    collapsed_tstamps = D[]
    collapsed_values = ta.values[1:0, :]

    tstamp = ta.timestamp[1]
    mapped_tstamp = period(tstamp)
    cluster_startrow = 1

    for i in 1:length(ta)-1

        next_tstamp = ta.timestamp[i+1]
        next_mapped_tstamp = period(next_tstamp)

        if mapped_tstamp != next_mapped_tstamp
          push!(collapsed_tstamps, timestamp(ta.timestamp[cluster_startrow:i]))
          collapsed_values = [collapsed_values; T[value(ta.values[cluster_startrow:i, j]) for j in 1:ncols]']
          cluster_startrow = i+1
        end #if

        tstamp = next_tstamp
        mapped_tstamp = next_mapped_tstamp

    end #for

    push!(collapsed_tstamps, timestamp(ta.timestamp[cluster_startrow:end]))
    collapsed_values = [collapsed_values; T[value(ta.values[cluster_startrow:end, j]) for j in 1:ncols]']

    N == 1 && (collapsed_values = vec(collapsed_values))
    return TimeArray(collapsed_tstamps, collapsed_values, ta.colnames, ta.meta)

end

# vcat ######################

function vcat(TA::TimeArray...)
    # Check all meta fields are identical.
    prev_meta = TA[1].meta
    for ta in TA
        if ta.meta != prev_meta
            throw(ArgumentError("metadata doesn't match"))
        end
    end

    # Check column names are identical.
    prev_colnames = TA[1].colnames
    for ta in TA
        if ta.colnames != prev_colnames
            throw(ArgumentError("column names don't match"))
        end
    end

    # Concatenate the contents.
    timestamps = vcat([ta.timestamp for ta in TA]...)
    values = vcat([ta.values for ta in TA]...)

    order = sortperm(timestamps)
    if length(TA[1].colnames) == 1 # Check for 1D to ensure values remains a 1D vector.
        return TimeArray(timestamps[order], values[order], TA[1].colnames, TA[1].meta)
    else
        return TimeArray(timestamps[order], values[order, :], TA[1].colnames, TA[1].meta)
    end
end

# map ######################

function map(f::Function, ta::TimeArray)
    timestamps = similar(ta.timestamp)
    values = similar(ta.values)

    for i in 1:length(ta)
        timestamps[i], values[i, :] = f(ta.timestamp[i], vec(ta.values[i, :]))
    end

    order = sortperm(timestamps)
    if length(ta.colnames) == 1 # Check for 1D to ensure values remains a 1D vector.
        return TimeArray(timestamps[order], values[order], ta.colnames, ta.meta)
    else
        return TimeArray(timestamps[order], values[order, :], ta.colnames, ta.meta)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

@userplot ParetoPlot

@recipe function f(p::ParetoPlot)
  # get user input
  data = p.args[1]

  x = sort(data, rev=true)
  n = length(x)
  logp = log.([i/(n+1) for i in 1:n])
  logx = log.(x)

  seriestype --> :scatter
  xguide --> "-log(i/(n+1))"
  yguide --> "log(x*)"

  -logp, logx
end

ext: jl


text:
module BenchIndexSet

using BenchmarkTools
using ITensors

suite = BenchmarkGroup()

i = Index(2)
is = IndexSet(n -> i^(n - 1), Order(10))
is_tuple = ntuple(n -> i^(n - 1), Val(10))

suite["constructor"] = BenchmarkGroup()

suite["constructor"]["function"] = @benchmarkable IndexSet($(n -> i^(n - 1)), $(Order(6)))
suite["constructor"]["function, tuple"] = @benchmarkable ntuple(
  $(n -> i^(n - 1)), $(Val(6))
)

suite["filter"] = BenchmarkGroup()

suite["filter"]["kwargs"] = @benchmarkable filter($is; plev=3)
suite["filter"]["function"] = @benchmarkable filter($(i -> plev(i) < 2), $is)
suite["filter"]["function, tuple"] = @benchmarkable filter($(i -> plev(i) < 2), $is_tuple)

i, j, k, l = Index.(2, ("i", "j", "k", "l"))

Iij = (i, j)
Ijl = (j, l)
Ikl = (k, l)
Iijk = (i, j, k)

suite["set_functions"] = BenchmarkGroup()

suite["set_functions"]["uniqueinds"] = BenchmarkGroup()
suite["set_functions"]["uniqueinds"]["nofilter2"] = @benchmarkable uniqueinds($Iijk, $Ikl)
suite["set_functions"]["uniqueinds"]["nofilter0"] = @benchmarkable uniqueinds($Iij, $Iijk)
suite["set_functions"]["uniqueinds"]["filter_tags"] = @benchmarkable uniqueinds(
  $Iijk, $Ikl; tags=$(ts"i")
)
suite["set_functions"]["uniqueinds"]["filter_not_tags"] = @benchmarkable uniqueinds(
  $Iijk, $Ikl; tags=$(not("i"))
)
suite["set_functions"]["uniqueinds"]["3_inputs"] = @benchmarkable uniqueinds(
  $Iijk, $Ijl, $Ikl
)

A = randomITensor(i'', i', i)

suite["set_functions"]["prime, ITensor"] = @benchmarkable prime($A)
suite["set_functions"]["uniqueinds, ITensor"] = @benchmarkable uniqueinds($(A'), $A)
suite["set_functions"]["commoninds, ITensor"] = @benchmarkable commoninds($(A'), $A)
suite["set_functions"]["unioninds, ITensor"] = @benchmarkable unioninds($(A'), $A)
suite["set_functions"]["noncommoninds, ITensor"] = @benchmarkable noncommoninds($(A'), $A)

end

BenchIndexSet.suite


text:
using Documenter, PAPI

makedocs(;
    modules = [PAPI],
    sitename = "PAPI.jl",
    authors = "Tom Haber <tom.haber@uhasselt.be>",
    pages = [
        "Home" => "index.md",
        "Events" => "events.md",
        "Event Sets" => "eventsets.md",
        "Counting" => "counters.md",
        "Library" => [
            "Public" => "public.md",
            "Internals" => "internals.md",
        ]
    ]
)

deploydocs(;
    repo="github.com/tomhaber/PAPI.jl.git"
)


text:
# ---
# title: 202. Happy Number
# id: problem202
# author: zhwang
# date: 2022-03-22
# difficulty: Easy
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/happy-number/description/>
# hidden: true
# ---
#
# Write an algorithm to determine if a number `n` is "happy".
#
# A happy number is a number defined by the following process: Starting with any
# positive integer, replace the number by the sum of the squares of its digits,
# and repeat the process until the number equals 1 (where it will stay), or it
# **loops endlessly in a cycle** which does not include 1. Those numbers for
# which this process **ends in 1** are happy numbers.
#
# Return True if `n` is a happy number, and False if not.
#
# **Example:  **
#
#
#
#     Input: 19
#     Output: true
#     Explanation: 12 + 92 = 82
#     82 + 22 = 68
#     62 + 82 = 100
#     12 + 02 + 02 = 1
#
#
#
## @lc code=start
using LeetCode
function ishappy(n::Int)::Bool
    while n != 1
        n == 4 && return false ## all loops must contain 4
        n = sum(digits(n) .^ 2)
    end
    true
end
## @lc code=end

## @lc test=start
@testset "202.happy-number.jl" begin
    @test all(ishappy(i) for i in (1, 7, 10, 13, 19, 23, 28, 31, 123456))
    @test !any(ishappy(i) for i in 2:6)
    @test !any(ishappy(i) for i in (8, 9, 11, 12, 12345))
end
## @lc test=end

## @lc info=start
# link: [solution to 202]()
## @lc info=end


text:
import Base: show, push!

import .config.get_qubit_params
import .config.get_marker_params
import .config.get_edge_params
import .config.get_instrument_params

export compile_to_hardware

immutable Event
	label::String
end

SequenceEntry = Union{Pulse, PulseBlock, ZPulse, ControlFlow, Event}

function flatten_seqs{T}(seqs::Vector{Vector{T}})
	flat_seq = Vector{SequenceEntry}()
	for seq = seqs
		push!(flat_seq, wait())
		for e in seq
			push!(flat_seq, e)
		end
	end
	return flat_seq
end

function compile_to_hardware{T}(seqs::Vector{Vector{T}}, base_filename)
	compile_to_hardware(flatten_seqs(seqs), base_filename)
end

function compile_to_hardware{T}(seq::Vector{T}, base_filename; suffix="")

	# TODO: save input code to file as metadata
	#save_code(seq)

	# add slave trigger to every WAIT
	add_slave_trigger!(seq, Marker("slave_trig"))

	# propagate frame changes to edges
	propagate_frame_change!(seq)

	# TODO: add gating/blanking pulses
	#add_gate_pulses!(seq)

	# TODO: gating constraints
	# TODO: move to device drivers
	#apply_gating_constraints!(seq)

	seqs, pulses, chans = compile(seq)

	# normalize and inject the channel delays
	qubit_params = get_qubit_params()
	instr_params = get_instrument_params()
	chan_delays = Dict{Channel, Float64}()
	for chan in chans
		# delays are shifted by AWG delay too
		#chan_awg = channel_params[chan.awg_channel]["instrument"]
		ch_type = typeof(chan) == Marker ? "markers" : "tx_channels"
		chan_delays[chan] = instr_params[split(chan.awg_channel)[1]][ch_type][split(chan.awg_channel)[2]]["delay"]
		#chan_delays[chan] = channel_params[chan.awg_channel]["delay"] + awg_delay $TODO: is this obsolete?
	end
	normalize_channel_delays!(chan_delays)
	inject_channel_delays!(seqs, pulses, chan_delays)

	# map the labeled channels to physical channels and bundle per APS/AWG
	AWGs = Dict{String, Dict}()
	chan_str_map = Dict("12"=>:ch12, "12m1"=>:m1, "12m2"=>:m2, "12m3"=>:m3, "12m4"=>:m4)
	for chan in chans
		awg = split(chan.awg_channel)[1]
		# get channel string from AWG-chstr convention
		chan_str = split(chan.awg_channel)[2]
		# TODO: map is currently only for APS2 - should be looked up from somewhere
		# there can be multiple logical channels mapped to the same physical channel
		if haskey(AWGs, awg) && haskey(AWGs[awg], chan_str_map[chan_str])
			push!(AWGs[awg][chan_str_map[chan_str]], chan)
		else
			get!(AWGs, awg, Dict{Symbol, Array{QGL.Channel}}())[chan_str_map[chan_str]] = [chan]
	  end
	end

	translator_map = Dict("APS2" => APS2)
	for (awg, ch_map) in AWGs
		# use first channel to lookup translator
		translator = translator_map[instr_params[awg]["type"]]
		translator.write_sequence_file(base_filename*"-$awg.h5", seqs, pulses, ch_map)
	end
	return seqs
end

function add_slave_trigger!(seq, slave_trig_chan)
	wait_entry = wait()
	slave_trig = Pulse("TRIG", slave_trig_chan, slave_trig_chan.shape_params[:length], 1.0)
	for (ct,e) in enumerate(seq)
		if e == wait_entry
			# try to add to next entry with non-zero length
			ct2 = 1
			while true
				if length(seq[ct+ct2]) > 0
					seq[ct+ct2] = slave_trig ⊗ seq[ct+ct2]
					break
				elseif (typeof(seq[ct+ct2]) == ControlFlow)
					# if we've hit another ControlFlow we'll just have to inject
					insert!(seq, ct+ct2, slave_trig)
					break
				end
				ct2 += 1
				if ct+ct2 > length(seq)
					error("Could not find a place to add slave trigger")
				end
			end
		end
	end
end

function propagate_frame_change!(seq)
	#get a dictionary mapping qubits to edges in the sequence which hold them as target
	chans = channels(seq)
	seq_edges = filter(x -> typeof(x) == Edge, chans)
	qs = qubits(seq)
	edges = Dict(q => Set{Edge}(filter(e -> e.target == q, seq_edges)) for q in qs)

	# if there are no edges then we're finished here
	all(isempty(s) for s in values(edges)) && return

	# keep a cache of propagated Z pulses
	Z_edges = Dict{Tuple{Qubit, Float64}, PulseBlock}()

	#for any ZPulse at the target qubit, add ZPulse on its edges
	for (ct,entry) in enumerate(seq)
		if typeof(entry) == PulseBlock
			for (ch, e) in entry.pulses
				for pulse in e
					if typeof(pulse) == ZPulse && typeof(ch) == Qubit && ~isempty(edges[ch])
						if !((ch, pulse.angle) in keys(Z_edges))
							Z_edges[(ch, pulse.angle)] = reduce(⊗, [Z(x, pulse.angle) for x in edges[ch]])
						end
						insert!(seq, ct+1, Z_edges[(ch, pulse.angle)])
					end
				end
			end
		elseif typeof(entry) == ZPulse && typeof(entry.channel) == Qubit && ~isempty(edges[entry.channel])
			if !((entry.channel, entry.angle) in keys(Z_edges))
				Z_edges[(entry.channel, entry.angle)] = reduce(⊗, [Z(x, entry.angle) for x in edges[entry.channel]])
			end
			insert!(seq, ct+1, Z_edges[(entry.channel, entry.angle)])
		end
	end
end

function compile(seq)
	# find what channels we're dealing with here
	chans = channels(seq)
	seqs = Vector{SequenceEntry}()
	pulses = Dict(chan => Set{Pulse}() for chan in chans)
	paddings = Dict(chan => 0.0 for chan in chans)

	# step through sequence and schedule
	for entry in seq
		if typeof(entry) == ControlFlow
			# control flow resets clock
			if length(seqs) > 0 && typeof(seqs[end]) != ControlFlow
				for chan in chans
					apply_padding!(chan, seqs[end], paddings, pulses)
				end
			end
			push!(seqs, entry)
		else
			if (length(seqs) == 0) || (typeof(seqs[end]) == ControlFlow)
				push!(seqs, PulseBlock(chans))
			end
			push!(seqs[end], entry, pulses, paddings)
		end
	end

	# assume we want to loop to the beginning
	for chan in chans
		apply_padding!(chan, seqs[end], paddings, pulses)
	end
	push!(seqs, goto(0))

	return seqs, pulses, chans
end

function channels(seq)
	chans = Set{Channel}()
	for e in seq
		if typeof(e) == Pulse
			push!(chans, e.channel)
		elseif typeof(e) == PulseBlock
			for chan in keys(e.pulses)
				push!(chans, chan)
			end
		end
	end
	return chans
end

function qubits(seq)
	chans = channels(seq)
	for chan in chans
		if typeof(chan) != Qubit
			delete!(chans, chan)
		end
	end
	return chans
end

"""
Shifts global channel delay such that all delays are >= 0
"""
function normalize_channel_delays!(chan_delays)
	min_delay = minimum(values(chan_delays))
	for (c,d) in chan_delays
		chan_delays[c] -= min_delay
	end
end

function inject_channel_delays!(seqs, pulses, chan_delays)

	all(map(x -> x == 0, values(chan_delays))) && return

	delay_block = PulseBlock(collect(keys(chan_delays)))

	for (c,d) in chan_delays
		if d > 0
			p = QGL.Id(c, d)
			push!(delay_block.pulses[c], p)
			push!(pulses[c], p)
		end
	end

	for ct = length(seqs):-1:1
		if typeof(seqs[ct]) == QGL.ControlFlow && (seqs[ct].op == QGL.WAIT || seqs[ct].op == QGL.SYNC)
			insert!(seqs, ct+1, delay_block)
		end
	end
end

""" Append a Pulse to a specific channel in a PulseBlock """
function push_channel!(pb::PulseBlock, p::Pulse, pulses, paddings)
	length(p) == 0.0 && return 	# elide zero-length pulses
	apply_padding!(p.channel, pb, paddings, pulses)
	push!(pb.pulses[p.channel], p)
	push!(pulses[p.channel], p)
end

""" Append a ZPulse to a specific channel in a PulseBlock """
function push_channel!(pb::PulseBlock, zp::ZPulse, pulses, paddings)
	zp.angle == 0.0 && return 	# elide noop Z pulses
	apply_padding!(zp.channel, pb, paddings, pulses)
	push!(pb.pulses[zp.channel], zp)
end

""" Concatenate two PulseBlocks """
function push!(pb_cur::PulseBlock, pb_new::PulseBlock, pulses, paddings)
	pb_new_length = length(pb_new)
	for chan in channels(pb_cur)
		if chan in channels(pb_new)
			for p in pb_new.pulses[chan]
				push_channel!(pb_cur, p, pulses, paddings)
			end
			paddings[chan] += pb_new_length - sum(length(p) for p in pb_new.pulses[chan])
		else
			paddings[chan] += pb_new_length
		end
	end
end

""" Append a Pulse to a PulseBlock """
function push!(pb_cur::PulseBlock, p::Pulse, pulses, paddings)
	length(p) == 0.0 && return 	# elide zero-length pulses
	# push pulse onto appropriate channel and pad the other channels
	for chan in channels(pb_cur)
		if chan == p.channel
			push_channel!(pb_cur, p, pulses, paddings)
		else
			paddings[chan] += length(p)
		end
	end
end

""" Append a ZPulse to a PulseBlock """
push!(pb_cur::PulseBlock, zp::ZPulse, pulses, paddings) =	push_channel!(pb_cur, zp, pulses, paddings)

function apply_padding!(chan, pb, paddings, pulses)
	if paddings[chan] > 1e-16 #arbitrarily eps for Float64 relative to 1.0
		pad_pulse = Id(chan, paddings[chan])
		push!(pb.pulses[chan], pad_pulse)
		push!(pulses[chan], pad_pulse)
		paddings[chan] = 0.0
	end
end


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = [1.0, 1.0, 1.0, 1.0];
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    gam = mu[1];
    gam1 = gam - 1.0;
    r = u[1];
    ru = u[2];
    rv = u[3];
    rE = u[4];
    r1 = 1/r;
    uv = ru*r1;
    vv = rv*r1;
    E = rE*r1;
    p = gam1*(rE-r*0.5*(uv*uv+vv*vv));
    h = E+p*r1;
    f = [ru, ru*uv+p, rv*uv, ru*h, rv, ru*vv, rv*vv+p, rv*h];
    f = reshape(f,(4,2));
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = [0.0, 0.0, 0.0, 0.0];
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = [0.0, 0.0, 0.0, 0.0];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[:,1]*n[1] + f[:,2]*n[2] + tau[1]*(u-uhat);
    return fb;
end
function initu(x, mu, eta)

    t = 0.0;
    phi = 5.0;
    x1 = x[1];
    x2 = x[2];
    gam = mu[1];
    M_ref = mu[2];

    r = sqrt((x1-t)^2 + x2^2);
    u01 = (1 - ((gam-1)/(16*pi^2)) * phi^2 * exp(2*(1-r^2)) )^(1/(gam-1));
    u02 = u01 * (1 - M_ref^(-1) * phi * exp(1-r^2) * x2/(2*pi));
    u03 = u01 * M_ref^(-1) * phi * exp(1-r^2) * x1/(2*pi);
    p = u01^gam/(gam*M_ref^2);
    u04 = p/(gam-1) + 0.5 * (u02*u02 + u03*u03) / u01;

    u0 = [u01, u02, u03, u04];
    return u0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using FountainCodes, Test

@test (M = QMatrix{Int}(64, 10); M[1,1] == 0)
@test (M = QMatrix{Int}(64, 10); M[1,1] = 0; M[1,1] == 0)
@test (M = QMatrix{Int}(64, 10); M[1,1] = 10; M[1,1] == 10)
@test (M = QMatrix{Int}(64, 10); M[1,1] = 1; M[1,1] == 1)
@test (M = QMatrix{Int}(64, 10); M[1,1] = 1; M[2,2] = 1; subtract!(M,1,2); M[2,1] == 1)
@test (M = QMatrix{Int}(64, 10); M[1,1] = 1; M[2,2] = 1; subtract!(M,5,1,2); M[2,1] == -5)
@test (M = QMatrix{Int}(64, 10); M[64,1] = 1; resize!(M, 128, 10); M[64,1] == 1 && M[65,1] == 0)
@test (M = QMatrix{Int}(64, 10); M[64,10] = 1; resize!(M, 64, 20); M[64,10] == 1 && M[64,11] == 0)
@test (M = QMatrix{Int}(64, 10); M[64,10] = 1; resize!(M, 128, 10); M[64,10] == 1)
@test (M = QMatrix{Int}(64, 10); M[64,10] = 1; resize!(M, 256, 20); M[64,10] == 1)
@test (M = QMatrix{Int}(64, 10); M[64,10] = 2; resize!(M, 256, 20); M[64,10] == 2 && M[65,11] == 0)
@test (M = QMatrix{Int}(64, 10); M[64,10] = 1; M[:,10] = 0; M[64,10] == 0)
@test (M = QMatrix{Int}(64, 10); M[64,10] = 1; M[:,10] = 2; M[1,10] == 2)
@test (M = QMatrix{Int}(64, 10); M[:,1] = 1; M[:,10] = 1; subtract!(M, 1, 1, 10); iszero(getcolumn(M, 1)))

# tests for subtracting one column from another

"""test both cols binary"""
function test_subtract_1()
    m, n = 64, 10
    M = QMatrix{Float64}(m, n)
    M[:, 2] .= true
    subtract!(M, 1, 2)
    if M[:, 1] != trues(m)
        error("expected $(trues(m)), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_1()

"""test target col binary and source col non-binary"""
function test_subtract_2()
    m, n = 64, 10
    v = 1.2
    M = QMatrix{Float64}(m, n)
    M[:, 2] .= v
    subtract!(M, 1, 2)
    if !all(M[:, 1] .== -v)
        error("expected all values to be $(-v), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_2()

"""test source col binary and target col non-binary"""
function test_subtract_3()
    m, n = 64, 10
    v = 1.2
    M = QMatrix{Float64}(m, n)
    M[:, 1] .= v
    M[:, 2] .= true
    subtract!(M, 1, 2)
    if !all(M[:, 1] .== (v-1))
        error("expected all values to be $(v-1), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_3()

"""test both target and source col non-binary"""
function test_subtract_4()
    m, n = 64, 10
    v1 = 1.2
    v2 = 2.0
    M = QMatrix{Float64}(m, n)
    M[:, 1] .= v1
    M[:, 2] .= v2
    subtract!(M, 1, 2)
    if !all(M[:, 1] .== (v1-v2))
        error("expected all values to be $(v1-v2), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_4()

# tests for subtracting a multiple of one column from another

"""test both cols binary and coef = 1"""
function test_subtract_coef_1()
    m, n = 64, 10
    M = QMatrix{Float64}(m, n)
    M[:, 2] .= true
    subtract!(M, 1.0, 1, 2)
    if !FountainCodes.isbinary(M, 1)
        error("expected column 1 to be binary, but got $(M[:, 1])")
    end
    if M[:, 1] != trues(m)
        error("expected $(trues(m)), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_coef_1()

"""test both cols binary and coef != 1"""
function test_subtract_coef_2()
    m, n = 64, 10
    M = QMatrix{Float64}(m, n)
    M[:, 2] .= true
    subtract!(M, 2.0, 1, 2)
    if FountainCodes.isbinary(M, 1)
        error("expected column 1 to be non-binary, but got $(M[:, 1])")
    end
    if !all(M[:, 1] .== (-2.0))
        error("expected all values to be -2, but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_coef_2()

"""test target col binary and source col non-binary"""
function test_subtract_coef_3()
    m, n = 64, 10
    v = 1.2
    M = QMatrix{Float64}(m, n)
    M[:, 2] .= v
    subtract!(M, 2.0, 1, 2)
    if !all(M[:, 1] .== -2.4)
        error("expected all values to be -2.4, but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_coef_3()

"""test source col binary and target col non-binary"""
function test_subtract_coef_4()
    m, n = 64, 10
    v = 1.2
    M = QMatrix{Float64}(m, n)
    M[:, 1] .= v
    M[:, 2] .= true
    subtract!(M, 2.0, 1, 2)
    if !all(M[:, 1] .== (1.2-2.0))
        error("expected all values to be $(1.2-2.0), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_coef_4()

"""test both target and source col non-binary"""
function test_subtract_coef_5()
    m, n = 64, 10
    v1 = 1.2
    v2 = 2.1
    M = QMatrix{Float64}(m, n)
    M[:, 1] .= v1
    M[:, 2] .= v2
    subtract!(M, 10.0, 1, 2)
    if !all(M[:, 1] .== (v1-10*v2))
        error("expected all values to be $(v1-10*v2), but got $(M[:, 1])")
    end
    return true
end
@test test_subtract_coef_5()

# resize tests

"""test that subtracted values persist after resize"""
function test_resize_1()
    M = QMatrix{UInt8}(64, 10)
    M[:,1] = UInt8(1)
    for i in 2:10
        subtract!(M, UInt8(1), i, 1)
    end
    for i in 2:10
        for (j, v) in enumerate(getcolumn(M, i))
            if v != one(v)
                error("M[$j, $i] should be one but is $v")
            end
        end
    end
    resize!(M, 128, 10)
    for i in 2:10
        subtract!(M, UInt8(1), i, 1)
    end
    for i in 2:10
        for (j, v) in enumerate(getcolumn(M, i))
            if !iszero(v)
                error("M[$j, $i] should be zero but is $v")
            end
        end
    end
    return true
end
@test test_resize_1()

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
############################################################################################
include("float.jl")
include("float_vector.jl")
include("float_array.jl")
include("integer.jl")

############################################################################################
# Export


text:
immutable TrajdataFrame
    lo::Int # index in states of first vehicle in scene
    hi::Int # index in states of last vehicle in scene
    t::Float64 # time
end
immutable TrajdataState
    id::Int # vehicle ID
    state::VehicleState
end
Base.length(frame::TrajdataFrame) = frame.hi - frame.lo + 1 # number of cars in the frame

type Trajdata
    roadway::Roadway
    vehdefs::Dict{Int, VehicleDef} # vehicle id -> vehdef
    states::Vector{TrajdataState} # list of vehicle states (for each scene)
    frames::Vector{TrajdataFrame} # list of frames
end
Trajdata(roadway::Roadway) = Trajdata(roadway, Dict{Int, VehicleDef}(), TrajdataState[], TrajdataFrame[])
function Trajdata(trajdata::Trajdata, frame_lo::Int, frame_hi::Int)
    roadway = trajdata.roadway
    vehdefs = Dict{Int, VehicleDef}()

    nstates = 0
    for frame in trajdata.frames[frame_lo : frame_hi]
        nstates += length(frame)
    end

    states = Array(TrajdataState, nstates)
    copy!(states, 1, trajdata.states, trajdata.frames[frame_lo].lo, nstates)
    for s in states
        if !haskey(vehdefs, s.id)
            vehdefs[s.id] = trajdata.vehdefs[s.id]
        end
    end

    frames = Array(TrajdataFrame, frame_hi - frame_lo + 1)
    Δstate_index = trajdata.frames[frame_lo].lo - 1
    for (i,frame) in enumerate(trajdata.frames[frame_lo : frame_hi])
        frames[i] = TrajdataFrame(frame.lo - Δstate_index, frame.hi - Δstate_index, frame.t)
    end

    Trajdata(roadway, vehdefs, states, frames)
end

function Base.write(io::IO, trajdata::Trajdata)
    # writes to a text file
    # - does not write the roadway

    println(io, "TRAJDATA")

    # vehdefs
    println(io, length(trajdata.vehdefs)) # number of vehdefs
    for def in values(trajdata.vehdefs)
        @printf(io, "%d %d %.3f %.3f\n", def.id, def.class, def.length, def.width)
    end

    # states
    println(io, length(trajdata.states)) # number of states
    for trajdatastate in trajdata.states
        state = trajdatastate.state
        @printf(io, "%d (%.4f %.4f %.4e) (%d %.4f %d %d) (%.4f %.4f %.4e) %.4f\n",
                     trajdatastate.id,
                     state.posG.x, state.posG.y, state.posG.θ,
                     state.posF.roadind.ind.i, state.posF.roadind.ind.t,
                     state.posF.roadind.tag.segment, state.posF.roadind.tag.lane,
                     state.posF.s, state.posF.t, state.posF.ϕ,
                     state.v
                )
    end

    # frames
    println(io, length(trajdata.frames)) # number of frames
    for tf in trajdata.frames
        @printf(io, "%d %d %.4f\n", tf.lo, tf.hi, tf.t)
    end
end
function Base.read(io::IO, ::Type{Trajdata}, roadway::Roadway=Roadway())
    lines = readlines(io)
    line_index = 1
    if contains(lines[line_index], "TRAJDATA")
        line_index += 1
    end

    function advance!()
        line = strip(lines[line_index])
        line_index += 1
        line
    end

    vehdefs = Dict{Int, VehicleDef}()
    N = parse(Int, advance!())
    for i in 1 : N
        tokens = split(advance!(), ' ')
        id = parse(Int, tokens[1])
        class = parse(Int, tokens[2])
        length = parse(Float64, tokens[3])
        width = parse(Float64, tokens[4])
        vehdefs[id] = VehicleDef(id, class, length, width)
    end

    N = parse(Int, advance!())
    states = Array(TrajdataState, N)
    for i in 1 : N
        line = advance!()
        cleanedline = replace(line, r"(\(|\))", "")
        tokens = split(cleanedline, ' ')
        id = parse(Int, tokens[1])
        x = parse(Float64, tokens[2])
        y = parse(Float64, tokens[3])
        θ = parse(Float64, tokens[4])
        ind_i = parse(Int, tokens[5])
        ind_t = parse(Float64, tokens[6])
        seg = parse(Int, tokens[7])
        lane = parse(Int, tokens[8])
        s = parse(Float64, tokens[9])
        t = parse(Float64, tokens[10])
        ϕ = parse(Float64, tokens[11])
        v = parse(Float64, tokens[12])
        states[i] = TrajdataState(id, VehicleState(VecSE2(x,y,θ), Frenet(RoadIndex(CurveIndex(ind_i, ind_t), LaneTag(seg, lane)), s, t, ϕ), v))
    end

    N = parse(Int, advance!())
    frames = Array(TrajdataFrame, N)
    for i in 1:N
        tokens = split(advance!(), ' ')
        lo = parse(Int, tokens[1])
        hi = parse(Int, tokens[2])
        t = parse(Float64, tokens[3])
        frames[i] = TrajdataFrame(lo, hi, t)
    end

    Trajdata(roadway, vehdefs, states, frames)
end

get_ids(trajdata::Trajdata) = collect(keys(trajdata.vehdefs))
nframes(trajdata::Trajdata) = length(trajdata.frames)
frame_inbounds(trajdata::Trajdata, frame::Int) = 1 ≤ frame ≤ nframes(trajdata)
carsinframe(trajdata::Trajdata, frame::Int) = length(trajdata.frames[frame])
nth_carid(trajdata::Trajdata, frame::Int, n::Int=1) = trajdata.states[trajdata.frames[frame].lo + n-1].id
get_time(trajdata::Trajdata, frame::Int) = trajdata.frames[frame].t
get_elapsed_time(trajdata::Trajdata, frame_lo::Int, frame_hi::Int) = trajdata.frames[frame_hi].t - trajdata.frames[frame_lo].t
get_mean_timestep(trajdata::Trajdata) = (trajdata.frames[end].t - trajdata.frames[1].t) / (nframes(trajdata)-1)

function get_first_frame_with_id(trajdata::Trajdata, id::Int)
    for (frame_index, frame) in enumerate(trajdata.frames)
        for i in frame.lo : frame.hi
            if trajdata.states[i].id == id
                return frame_index
            end
        end
    end

    -1
end
function get_last_frame_with_id(trajdata::Trajdata, id::Int)
    for frame_index in reverse(1:length(trajdata.frames))
        frame = trajdata.frames[frame_index]
        for i in frame.lo : frame.hi
            if trajdata.states[i].id == id
                return frame_index
            end
        end
    end

    -1
end

function iscarinframe(trajdata::Trajdata, id::Int, frame::Int)
    frame = trajdata.frames[frame]
    for i in frame.lo : frame.hi
        s = trajdata.states[i]
        if s.id == id
            return true
        end
    end
    false
end

function get_vehiclestate(trajdata::Trajdata, id::Int, frame::Int)
    frame = trajdata.frames[frame]
    for i in frame.lo : frame.hi
        s = trajdata.states[i]
        if s.id == id
            return s.state
        end
    end
    error("vehicle not found for id $id and frame $(frame)!")
end
get_vehicledef(trajdata::Trajdata, id::Int) = trajdata.vehdefs[id]

function get_vehicle!(veh::Vehicle, trajdata::Trajdata, id::Int, frame::Int)
    veh.state = get_vehiclestate(trajdata, id, frame)
    veh.def = get_vehicledef(trajdata, id)
    veh
end
get_vehicle(trajdata::Trajdata, id::Int, frame::Int) = get_vehicle!(Vehicle(), trajdata, id, frame)

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

immutable TrajdataVehicleIterator
    trajdata::Trajdata
    id::Int
end
Base.length(iter::TrajdataVehicleIterator) = sum(frame->iscarinframe(iter.trajdata, iter.id, frame), 1 : nframes(iter.trajdata))
function Base.start(iter::TrajdataVehicleIterator)
    frame = 1
    while frame < nframes(iter.trajdata) &&
          !iscarinframe(iter.trajdata, iter.id, frame)

        frame += 1
    end
    frame
end
Base.done(iter::TrajdataVehicleIterator, frame::Int) = frame > nframes(iter.trajdata)
function Base.next(iter::TrajdataVehicleIterator, frame::Int)
    item = (frame, get_vehicle(iter.trajdata, iter.id, frame))
    frame += 1
    while frame < nframes(iter.trajdata) &&
          !iscarinframe(iter.trajdata, iter.id, frame)

        frame += 1
    end
    (item, frame)
end


text:
using Knet.KnetArrays: KnetArray
using CUDA: CuArray
using AutoGrad: Value
const GPUVal{T,N} = Union{KnetArray{T,N},CuArray{T,N},Value{KnetArray{T,N}},Value{CuArray{T,N}}}

# cut-and-paste method: see Knet/prof/zeropad for alternatives
# cpu=147ms (x=(56,56,128,128), p=((0,1),(0,1)))

function zeropad1(x, p; y=nothing, fillzero=true)
    p = padtuple(p, size(x))
    n = ndims(x)
    d = ntuple(i->(i < n-1 ? size(x,i)+sum(p[i]) : size(x,i)), n)
    c = ntuple(i->(i < n-1 ? (p[i][1]+1:p[i][1]+size(x,i)) : Colon()), n)
    y === nothing ? y = similar(x,d) : @assert typeof(y)===typeof(x) && size(y)==d
    fillzero && fill!(y, 0)
    y[c...] .= x
    return y
end

# concatenate: 185ms (x=(56,56,128,128), p=((0,1),(0,1)))
function zeropad2(x, p)
    pad(d) = fill!(similar(x, d), 0)
    p = padtuple(p, size(x))
    for i in 1:length(p)
        pad1 = pad(ntuple(j->(j==i ? p[j][1] : size(x,j)), ndims(x)))
        pad2 = pad(ntuple(j->(j==i ? p[j][2] : size(x,j)), ndims(x)))
        x = cat(pad1, x, pad2; dims=i)
    end
    return x
end

# depthwise convolution: 351ms (x=(56,56,128,128), p=((0,1),(0,1)))
function zeropad3(x, p)
    p = padtuple(p, size(x))
    n = ndims(x)
    d = ntuple(i->(i < n-1 ? 1+sum(p[i]) : i == n-1 ? 1 : size(x,n-1)), n)
    w = fill!(similar(x, d), 0)
    c = ntuple(i->1+p[i][2], n-2)
    w[c...,:,:] .= 1
    conv(w, x; padding=sum.(p), groups=size(x,n-1))
end

# cut-and-paste: med=4.2ms mean=3.8ms (x=(56,56,128,128), p=((0,1),(0,1)))
function zeropad1(x::GPUVal, p)
    p = padtuple(p, size(x))
    n = ndims(x)
    d = ntuple(i->(i < n-1 ? size(x,i)+sum(p[i]) : size(x,i)), n)

    global _y            # assumes p will not change, good for layers
    if !@isdefined(_y) || typeof(_y) !== typeof(x) || size(_y) != d
        _y = similar(x,d)
        fill!(_y, 0)
    end

    y = _y # fill!(similar(x, d), 0)
    c = ntuple(i->(i < n-1 ? (p[i][1]+1:p[i][1]+size(x,i)) : Colon()), n)
    y[c...] .= x
    return y
end

# concatenate: med=18.4ms mean=21.1ms (x=(56,56,128,128), p=((0,1),(0,1)))
function zeropad2(x::GPUVal, p)
    pad(d) = fill!(similar(x, d), 0)
    p = padtuple(p, size(x))

    global _pad                # assume pad/type will not change
    if !@isdefined(_pad) || _pad === nothing
        _pad = Any[Any[nothing,nothing],Any[nothing,nothing]]
    end

    for i in 1:length(p)
        # pad1 = pad(ntuple(j->(j==i ? p[j][1] : size(x,j)), ndims(x)))
        # pad2 = pad(ntuple(j->(j==i ? p[j][2] : size(x,j)), ndims(x)))
        # x = cat(pad1, x, pad2; dims=i)
        if _pad[i][1] === nothing
            _pad[i][1] = pad(ntuple(j->(j==i ? p[j][1] : size(x,j)), ndims(x)))
            _pad[i][2] = pad(ntuple(j->(j==i ? p[j][2] : size(x,j)), ndims(x)))
        end
        x = cat(_pad[i][1], x, _pad[i][2]; dims=i)
    end
    return x
end

# depthwise convolution: med=6.0ms mean=5.4ms (x=(56,56,128,128), p=((0,1),(0,1)))
function zeropad3(x::GPUVal, p)
    p = padtuple(p, size(x))
    n = ndims(x)
    yd = ntuple(i->(i < n-1 ? size(x,i)+sum(p[i]) : size(x,i)), n)

    global _y            # assumes p will not change, good for layers
    if !@isdefined(_y) || typeof(_y) !== typeof(x) || size(_y) != yd
        _y = similar(x,yd)
        fill!(_y, 0)
    end

    global _wpad                # assume padding does not change for testing
    if !@isdefined(_wpad) || _wpad === nothing
        wd = ntuple(i->(i < n-1 ? 1+sum(p[i]) : i == n-1 ? 1 : size(x,n-1)), n)
        _wpad = fill!(similar(x, wd), 0)
        c = ntuple(i->1+p[i][2], n-2)
        _wpad[c...,:,:] .= 1
    end

    conv(_wpad, x; padding=sum.(p), groups=size(x,n-1))
end


text:
@registerN Meixner AbstractCDOP2 γ μ
export Meixner

"""
    Meixner{γ,μ}

References: [Koekoek and Swarttouw §1.9](https://arxiv.org/pdf/math/9602214.pdf)
"""
Meixner

basis_symbol(::Type{<:Meixner{γ,μ}}) where {γ,μ} = "m⁽ᵞᵝ⁾"
Polynomials.domain(::Type{<:Meixner{γ,μ}}) where {γ,μ} = Polynomials.Interval(0, Inf)
abcde(::Type{<:Meixner{γ,μ}}) where {γ,μ} =
    NamedTuple{(:a, :b, :c, :d, :e)}((0, 1, 0, μ - 1, γ * μ))

function kn(P::Type{<:Meixner{γ,μ}}, n::Int) where {γ,μ}
    (1 - 1 / μ)^n
end

function k1k0(P::Type{<:Meixner{γ,μ}}, n::Int) where {γ,μ}
    (1 - 1 / μ)
end
function k1k_1(P::Type{<:Meixner{γ,μ}}, n::Int) where {γ,μ}
    (1 - 1 / μ)^(n > 0 ? 2 : 1)
end

function classical_hypergeometric(P::Type{<:Meixner{γ,μ}}, n::Int, x) where {γ,μ}
    Pochhammer(γ, n) * pFq((-n, -x), γ, 1 - 1 / μ)
end

# Overrides


text: # import TimeSeries


text:
"""
    DigitSet([itr])

Construct a `DigitSet`, an `AbstractSet{Int}` which can only contains numbers 0:63.
A `DigitSet` is stored as a single integer in memory, and is immutable.

See also: [`StackSet`](@ref)
"""
struct DigitSet <: AbstractStackSet
    x::UInt

    DigitSet(x::UInt) = new(x)
end

DigitSet() = DigitSet(zero(UInt))

function DigitSet(itr)
    d = DigitSet()
    for i in itr
        d = push(d, convert(Int, i))
    end
    d
end

function Base.hash(x::DigitSet, h::UInt)
    base = 0x0cabd57b4f53416a % UInt
    return hash(x.x, h ⊻ base)
end

@noinline function throw_DigitSet_digit_err()
    throw(ArgumentError("DigitSet can only contain 0:$(Sys.WORD_SIZE-1)"))
end

function show(io::IO, s::AbstractStackSet)
    print(io, "$(typeof(s))(")
    Base.show_vector(io, s)
    print(io, ')')
end

Base.empty(x::DigitSet) = DigitSet()
Base.isempty(x::DigitSet) = iszero(x.x)
Base.:(==)(x::T, y::T) where T<:AbstractStackSet = x === y
Base.:⊊(x::AbstractStackSet, y::AbstractStackSet) = issubset(x, y) & (x != y)
Base.allunique(x::AbstractStackSet) = true

function Base.iterate(x::DigitSet, state::UInt=x.x)
	iszero(state) ? nothing : (trailing_zeros(state), state & (state - 1))
end

Base.in(x::Int, s::DigitSet) = isodd(s.x >>> unsigned(x))
Base.length(x::DigitSet) = count_ones(x.x)

Base.maximum(x::DigitSet, ::Unsafe) = Sys.WORD_SIZE - 1 - leading_zeros(x.x)
function Base.maximum(x::DigitSet)
    isempty(x) && throw_empty_err()
    return maximum(x, unsafe)
end

Base.minimum(x::DigitSet, ::Unsafe) = trailing_zeros(x.x)
function Base.minimum(x::DigitSet)
    isempty(x) && throw_empty_err()
    return minimum(x, unsafe)
end

function push(s::DigitSet, v::Int, ::Unsafe)
    DigitSet(s.x | (1 << (unsigned(v) & (Sys.WORD_SIZE - 1))))
end

"""
    push(collection, items...)

Return a new collection containing all elements of `collection` and of `items`.
If `collection` is ordered, add the new element to the end.

# Examples

```jldoctest
julia> s = DigitSet([4, 9]);

julia> push(s, 1, 11)
DigitSet with 4 elements:
  1
  4
  9
  11
```
"""
function push end

function push(s::DigitSet, v::Int)
    unsigned(v) ≥ Sys.WORD_SIZE ? throw_DigitSet_digit_err() : push(s, v, unsafe)
end

function push(s::DigitSet, vs...)
    for v in vs
        s = push(s, convert(Int, v))
    end
    s
end

function Base.filter(pred, x::DigitSet)
    r = DigitSet()
    for i in x
        pred(i) && (r = push(r, i, unsafe))
    end
    r
end

"""
    pop(s::AbstractStackSet, v::Int)

Return a copy of `s` without `v`. If `v` is not in `s`, raise an error.
# Examples

```jldoctest
julia> s = DigitSet([4, 41, 9]);

julia> pop(s, 41)
DigitSet with 2 elements:
  4
  9
```

See also: [`delete`](@ref)
"""
function pop end

pop(s::AbstractStackSet, v::Int) = in(v, s) ? delete(s, v, unsafe) : throw(KeyError(v))

"""
    delete(s::AbstractStackSet, v::Int)

Return a copy of `s` that does not contain without `v`.

# Examples

```jldoctest
julia> s = DigitSet([4, 41, 9]);

julia> delete(s, 1)
DigitSet with 3 elements:
  4
  9
  41
```

See also: [`pop`](@ref)
"""
function delete end

delete(s::DigitSet, v::Int) = ifelse((v < 0) | (v ≥ Sys.WORD_SIZE), s, delete(s, v, unsafe))
function delete(s::DigitSet, v::Int, ::Unsafe)
    mask = ~(UInt(1) << (unsigned(v) & (Sys.WORD_SIZE - 1)))
    DigitSet(s.x & mask)
end

Base.issubset(x::AbstractStackSet, y::AbstractStackSet) = isempty(setdiff(x, y))

# isdisjoint was added to Base in 1.5
@static if VERSION < v"1.5"
    """
        isdisjoint(x, y) -> Bool

    Check if `x` and `y` have no elements in common.

    # Examples

    ```jldoctest
    julia> isdisjoint(DigitSet([1,6,4]), DigitSet([0, 61, 44]))
    true

    julia> isdisjoint(DigitSet([1,6,4]), DigitSet([4, 61, 44]))
    false

    julia> isdisjoint(DigitSet(), DigitSet())
    true
    ```
    """
    function isdisjoint end
    isdisjoint(x::AbstractStackSet, y::AbstractStackSet) = isempty(intersect(x, y))
else
    Base.isdisjoint(x::AbstractStackSet, y::AbstractStackSet) = isempty(intersect(x, y))
end

Base.setdiff(x::DigitSet, y::DigitSet) = DigitSet(x.x & ~y.x)

for (func, op) in ((:union, :|), (:symdiff, :⊻), (:intersect, :&))
    @eval begin
        function Base.$(func)(x::DigitSet, y::DigitSet)
            DigitSet($op(x.x, y.x))
        end
    end
end

# These functions work even if the collections in `itrs` contains elements not
# in 0:63.
for func in (:intersect, :setdiff)
    @eval begin
        function Base.$(func)(x::AbstractStackSet, itrs...)
            y = x
            for i in itrs
                z = typeof(x)()
                for j in i
                    if j in y
                        z = push(z, j)
                    end
                end
                y = $func(y, z)
            end
            return y
        end
    end
end

# This functions must fail if any item in a collection in `itrs` are not in
# the applicable domain, e.g. 0:63 for DigitSet
function Base.union(x::AbstractStackSet, itrs...)
    for i in itrs
        y_ = typeof(x)(i)
        x = union(x, y_)
    end
    return x
end

# This function needs to extend past 0:63, because it needs to keep track of
# which elements have been added and/or canceled by the iterables
Base.symdiff(x::AbstractStackSet, itrs...) = typeof(x)(symdiff(Set(x), itrs...))


text:
include("dual_objective_value.jl")

mutable struct Solution
    feasible::Bool
    objective_value::Union{Nothing, Float64}
    # Map between the variable ref
    # and the primal value
    variable_value::Dict{JuMP.VariableRef, Float64}
    # Map between the parameter
    # and the dual value
    parameter_dual::Dict{ParameterRef, Float64}
end

function optimize(model::ParametrizedModel)
    JuMP.optimize!(model.model)
    feasible = JuMP.primal_status(model.model) == MOI.FEASIBLE_POINT
    if feasible
        @assert JuMP.termination_status(model.model) == MOI.OPTIMAL
        objective_value = JuMP.objective_value(model.model)
    else
        @assert JuMP.termination_status(model.model) in (MOI.INFEASIBLE, MOI.ALMOST_INFEASIBLE)
        if isempty(model.parameter_map)
            objective_value = nothing
        else
            if JuMP.termination_status(model.model) == MOI.ALMOST_INFEASIBLE
                @warn("Using infeasibility ray of result with termination status: ALMOST_INFEASIBLE.")
            end
            @assert JuMP.dual_status(model.model) == MOI.INFEASIBILITY_CERTIFICATE
            objective_value = dual_objective_value(JuMP.backend(model.model),
                                                   Float64)
        end
    end
    variable_value = Dict{JuMP.VariableRef, Float64}()
    if feasible
        for vref in values(model.variable_map)
            variable_value[vref] = JuMP.value(vref)
        end
        for θ in values(model.θ)
            variable_value[θ] = JuMP.value(θ)
        end
    end
    parameter_dual = Dict{ParameterRef, Float64}()
    for parameter in values(model.parameter_map)
        parameter_dual[parameter] = JuMP.dual(parameter)
    end
    Solution(feasible, objective_value, variable_value, parameter_dual)
end

function set_parent_solution!(model::ParametrizedModel, parent::ParametrizedModel, parent_solution::Solution)
    for (index, parameter) in model.parameter_map
        vref = parent.variable_map[index]
        value = parent_solution.variable_value[vref]
        JuMP.fix(parameter, value)
    end
end

function add_cutting_planes(master_model, master_solution, sub_models, sub_solutions, TOL)
    cut_added = false
    # add cutting planes, one per scenario
    for (id, sol) = sub_solutions
        sub_model = sub_models[id]
        aff = JuMP.AffExpr(sol.objective_value)
        for (index, parameter) in sub_model.parameter_map
            dual = sol.parameter_dual[parameter]
            vref = master_model.variable_map[index]
            aff.constant -= dual * master_solution.variable_value[vref]
            JuMP.add_to_expression!(aff, dual, vref)
        end
        if sol.feasible
            JuMP.add_to_expression!(aff, -1.0, master_model.θ[id])
            # Check if the cut is useful
            if JuMP.value(aff, vref -> master_solution.variable_value[vref]) - TOL < 0
                continue
            end
        else
            # It is a ray, so `alpha * aff` is also valid for any alpha >= 0.
            # Hence `aff` might have very large coefficients and alter
            # the numerical accuracy of the master's solver.
            # We scale it to avoid this issue
            scaling = abs(aff.constant)
            if iszero(scaling)
                scaling = maximum(term -> term[1], JuMP.linear_terms(aff))
            end
            scaled_aff = typeof(aff)(sign(aff.constant))
            for (coef, var) in JuMP.linear_terms(aff)
                JuMP.add_to_expression!(scaled_aff, coef / scaling, var)
            end
            aff = scaled_aff
        end
        @constraint(master_model.model, aff <= 0)
        cut_added = true
    end
    return cut_added
end

function Benders_pmap(master_model::ParametrizedModel,
                      sub_models::Dict{Int, ParametrizedModel},
                      TOL=1e-5)::Solution
    master_solution = nothing
    cut_added = true
    while cut_added
        master_solution = optimize(master_model)
        if !master_solution.feasible
            break
        end

        for (id, sub_model) in sub_models
            set_parent_solution!(sub_model, master_model, master_solution)
        end

        sub_solutions = Dict(map(x -> (x[1] => optimize(x[2])), collect(sub_models))...)
        cut_added = add_cutting_planes(master_model, master_solution, sub_models, sub_solutions, TOL)
    end
    return master_solution
end


text:
"""
    load_transcript(path_bed12::String)

Loads transcripts information from a BED12-format file

# Arguments
- `path_bed12::String`: the path of the file
"""
function load_transcript(path_bed12::String)
	transcripts = Array{BED.Record,1}()
	open(BED.Reader, path_bed12) do bed12_reader
		for record::BED.Record in bed12_reader
			push!(transcripts, record)
		end
	end
	return transcripts
end


text:
# using Parameters

include("activations.jl")
include("cc.jl")
include("sec.jl")
include("pec.jl")

# torque generator
mutable struct TorqueGenerator
    τ::Float64  # torque
    θ::Float64  # angle
    ω::Float64  # angular velocity
    cc::CC      # contractile component
    sec::SEC    # series elastic component
    # pec::PEC    # parallel elastic component

    function TorqueGenerator(θ₀, ω₀, cc_p, sec_p, α_p)
        # initialise CC with joint angle and angular velocity
        cc = CC(θ₀, ω₀, cc_p, α_p)

        # SEC therefore has zero angle and angular velocity
        sec = SEC(0., 0., sec_p)

        # create new torque generator
        tqgen = new(0., θ₀, ω₀, cc, sec)

        # initialise CC, SEC angle and torque
        initccang!(tqgen)                   # CC angle
        tqgen.sec.θ = tqgen.θ - tqgen.cc.θ  # SEC angle
        torque!(tqgen.cc, 0.0)    # CC torque
        torque!(tqgen.sec) # SEC torque
        tqgen.τ = tqgen.sec.τ               # torque generator torque

        return tqgen
    end
end

# find initial CC angle
function initccang!(tqgen::TorqueGenerator)
    @unpack θ, cc, sec = tqgen

    # initial activation
    α = activation(0.0, cc.α_p)

    # catch if activation zero
    α == 0 && (cc.θ = 0; return θ)

    # closure returns difference in torque between CC and SEC for given CC angle
    f(θcc) = torque(θ - θcc, sec.sec_p) - torque(α, θcc, -cc.ω, cc.cc_p)

    # ensure CC angle +ve - not sure why needed?
    _θ = abs(θ)

    # set up min and max bounds of CC angle
    max_stretch = cc.cc_p.Tmax / sec.sec_p.k
    θmin = _θ - max_stretch
    θmax = _θ

    # solve for CC ang
    tol = 1e-6
    maxiters = 50
    θcc = bisection(f, θmin, θmax, tol, maxiters)

    # assign to CC
    θ < 0.0 && (θcc *= -1)    # resets sign
    cc.θ = θcc

    return cc.θ
end

# bisection method for root finding, stops after maxiters iterations or interval smaller than tol
function bisection(f, a, b, tol, maxiters)
    n = 1
    midpoint = 0.
    while n ≤ maxiters
        midpoint = a + (b - a) / 2
        # abs(b - a) < tol && return midpoint    # return midpoint if interval smaller than tol
        abs(f(midpoint)) < tol && return midpoint # return midpoint if it is smaller than tol
        f(a) * f(midpoint) < 0 ? b = midpoint : a = midpoint  # update interval
        n += 1
    end

    return midpoint
end

# numerically integrate the CC angle independently
function update_torque_generator_ind!(tqgen::TorqueGenerator, dt)
    @unpack θ, ω, cc, sec = tqgen

    # integrate CC angle to next time step
    # solve it here easy but inaccurate, otherwise add CC angles to integrator
    cc.θ += cc.ω * dt
    cc.θ ≥ θ && (cc.θ = θ)

    ### checks on cc angle? ###

    # calculate SEC angle
    sec.θ = θ - cc.θ

    # determine SEC torque and set to torque generator torque
    tqgen.τ = torque!(sec)

    # set SEC and CC torque equal
    cc.τ = sec.τ

    # determine CC angular velocity from torque
    τ_θ = torque_angle(cc.θ)
    τ_θ < 0.01 && (τ_θ = 0.01)    # stop τ_v blowing up if too small
    α < 0.01 && no_activation!(tqgen) # conditions for when activation low
    τ_ω = τ / (cc + p.T0 * α * τ_θ)   # calculate normalised torque from τ-ω relationship
    Tmax_norm = cc_p.Tmax / cc_p.T0
    τ_ω > Tmax_norm && (τ_ω = Tmax_norm)  # constrain normalised torque to realistic values
    τ_ω < 0.0 && (τ_ω = 0.0)
    cc.τ = τ_ω  # set CC torque
    cc.ω = -angvel(cc, cc_p)   # calculate CC angular velocity (opposite sign to function)
    sec.ω = ω - cc.ω    # update SEC angvel

    return tqgen.τ
end

# accepts CC angle as argument, so included in model integration
function update_torque_generator!(tqgen::TorqueGenerator, t, θ, ω, θcc)
    @unpack cc, sec = tqgen
    cc_p, sec_p = cc.cc_p, sec.sec_p

    # update torque generator angle and angular velocity
    tqgen.θ = θ
    tqgen.ω = ω

    # update CC angle
    θcc > tqgen.θ ? cc.θ = θcc : cc.θ = tqgen.θ
    cc.θ = θcc

    # calculate activation level
    α = activation(t, cc.α_p)
    cc.α = α
    α < 0.01 && (α = 0.01) #no_activation!(tqgen)

    ### when rising above α = 0.01 get sudden change in cc angle and angular velocity ###

    ### what if constrained activations to be min of 0.01 ###

    # calculate SEC angle
    θsec_max = cc_p.Tmax / sec_p.k
    sec.θ = θ - cc.θ
    sec.θ > θsec_max && (sec.θ = θsec_max)  # constrain SEC angle to max stretch

    # determine SEC torque and set to torque generator torque
    tqgen.τ = torque!(sec)

    # set SEC and CC torque equal
    cc.τ = sec.τ

    # calculate τ_θ component from CC angle
    torque_angle!(cc)
    cc.τ_θ < 0.01 && (cc.τ_θ = 0.01)    # stop τ_ω blowing up if too small

    # calculate τ-ω component from CC angular velocity
    cc.τ_ω = cc.τ / (cc_p.T0 * α * cc.τ_θ)
    (cc.τ_ω * cc_p.T0) > cc_p.Tmax && (cc.τ_ω = cc_p.Tmax / cc_p.T0)  # constrain τ_ω to realistic values
    cc.τ_ω < 0.0 && (cc.τ_ω = 0.0)

    # calculate CC angular velocity (opposite sign to function)
    angvel_torque_norm!(cc)
    sec.ω = ω - cc.ω    # update SEC angvel for completeness but not needed

    return tqgen.τ
end

# returns dictionary of CC and SEC components for given t, θ, ω and θcc
function get_torque_generator(tqgen, t, θ, ω, θcc)
    _tqgen = deepcopy(tqgen)
    update_torque_generator!(_tqgen, t, θ, ω, θcc)

    return Dict(
        :cc => Dict(
            field => getfield(_tqgen.cc, field) for field in fieldnames(CC)[1:end - 2]
        ),

        :sec => Dict(
            field => getfield(_tqgen.sec, field) for field in fieldnames(SEC)[1:end - 1]
        )
    )
end

function torque(tqgen::TorqueGenerator, θ, θcc)
    θsec = θ - θcc
    # constrain maximum SEC stretch to Tmax/k
    θsec_max = tqgen.cc.cc_p.Tmax / tqgen.sec.sec_p.k
    θsec > θsec_max ? τ = θsec_max : τ = torque(θsec, tqgen.sec.sec_p)

    return τ
end

# state of torque generator when no activation
function no_activation!(tqgen::TorqueGenerator)
    @unpack cc, sec = tqgen
    cc.α = 0.0
    cc.τ = 0.0
    sec.θ = 0.0     # no torque ∴ no SEC stretch
    torque!(sec)
    cc.θ = tqgen.θ  # CC angle must make up torque generator angle
    sec.ω = 0.0     # SEC angle not changing ∴ ang vel = 0
    cc.ω = tqgen.ω  # CC ang vel must make up torque generator ang vel

    # update CC torque components
    torque_angle!(cc)
    torque_angvel!(cc)
end


text:
using Base.Iterators: product

"""
    RectilinearGrid{Tl, Tr} <: AbstractVector{Tuple{Tl, Tr}}

A `RectilinearGrid` is parametrised by `AbstractVector`s of points `xl` and `xr`, whose
element types are `Tl` and `Tr` respectively, comprising `length(xl) * length(xr)`
elements. Linear indexing is the same as `product(eachindex(xl), eachindex(xr))` - `xl`
iterates more quickly than `xr`.
"""
struct RectilinearGrid{
    Tl, Tr, Txl<:AbstractVector{Tl}, Txr<:AbstractVector{Tr},
} <: AbstractVector{Tuple{Tl, Tr}}
    xl::Txl
    xr::Txr
end

Base.size(X::RectilinearGrid) = (length(X.xl) * length(X.xr),)

function Base.collect(X::RectilinearGrid{Tl, Tr}) where{Tl, Tr}
    return vec(
        map(
            ((p, q),) -> (X.xl[p], X.xr[q]),
            product(eachindex(X.xl), eachindex(X.xr)),
        )
    )
end

Base.show(io::IO, x::RectilinearGrid) = Base.show(io::IO, collect(x))

"""
    SpaceTimeGrid{Tr, Tt<:Real}

A `SpaceTimeGrid` is a `RectilinearGrid` in which the left vector corresponds to space, and
the right `time`. The left eltype is arbitrary, but the right must be `Real`.
"""
const SpaceTimeGrid{Tr, Tt<:Real} = RectilinearGrid{
    Tr, Tt, <:AbstractVector{Tr}, <:AbstractVector{Tt},
}

get_space(x::RectilinearGrid) = x.xl

get_time(x::RectilinearGrid) = x.xr


text:
using Random
using Distributions
using StaticArrays
using Dolo

@testset "Tests on outer product functions" begin

    @testset "testing the function outer of ergodic.jl" begin

        λ1, λ2 = rand(Uniform(0, 1), 2)

        @test Dolo.outer(SVector(1 - λ1, λ1), SVector(1 - λ2, λ2)) == @SMatrix [(1-λ1)*(1-λ2) (1-λ1)*λ2; λ1*(1-λ2) λ1*λ2]
        @test Dolo.outer(SVector(1 - λ1, λ1)) == @SVector(1 - λ1, λ1)

    end

    # @testset "testing the function outer2 of ergodic.jl" begin

    # end
end


text:
@testset "Generator" begin

    #
    m = 10
    density = 0.2
    net = WeightedNetwork.rand_network(m, density)
    adj = WeightedNetwork.AdjacencyList(net.weights)
    mat = WeightedNetwork.adjacency_matrix(adj)
    @test mat == net.weights

end

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
# Constant is nice but what if we have to use a variable
# because it needs to by mutated?  Let's solve it by using
# a Ref (scalar) value.

# Initialize constant Ref object with value of 10
const semi_constant = Ref(10)

#=
julia> const semi_constant = Ref(0)
Base.RefValue{Int64}(0)
=#

function add_using_global_semi_constant(x)
    return x + semi_constant[]
end

@btime add_using_global_semi_constant(10);
#= performance
julia> @btime add_using_global_semi_constant(10);
  2.097 ns (0 allocations: 0 bytes)
=#

# how to set Ref object value
semi_constant[] = 20
semi_constant
#=
julia> semi_constant[] = 20
20

julia> semi_constant
Base.RefValue{Int64}(20)
=#

# Option #2.
# How about mutable struct?
mutable struct Slot
    value::Int
end

const constant_slot = Slot(10)

function add_using_global_constant_slot(x)
    return constant_slot.value + x
end

add_using_global_constant_slot(10)

@btime add_using_global_constant_slot(10)
#=
julia> @btime add_using_global_constant_slot(10)
  2.097 ns (0 allocations: 0 bytes)
=#


text:
using JSON

mutable struct Config
    fn::String
    data::Dict{Any, Any}
    mtime::Number
    buffertime::Number
    lastCheck::Number
    checkDelay::Number

    Config(fn::String, data::Dict{K, V}, checkDelay::Number=2.5) where {K, V} = new(fn, data, 0, 0, 0, checkDelay)
    Config(fn::String, buffertime::Number, data::Dict{K, V}, checkDelay::Number=2.5) where {K, V} = new(fn, data, 0, buffertime, 0, checkDelay)
end

function saveConfig(c::Config, force::Bool=true)
    if force || c.buffertime <= 0 || c.mtime + c.buffertime < time()
        path = dirname(c.fn)

        if !isdir(path)
            mkdir(path)
        end

        # Save to $fn.saving first, then delete fn and move the temp file
        # Prevents corruption of data if program is terminated while writing
        open(c.fn * ".saving", "w") do fh
            c.mtime = time()
            JSON.print(fh, c.data, 4)
        end

        mv(c.fn * ".saving", c.fn, force=true)
    end
end

function loadConfig(fn::String, buffertime::Number=0, default::Dict{K, V}=Dict{Any, Any}()) where {K, V}
    tempFn = fn * ".saving"

    # Program terminated before the config was moved
    # Temp file has correct data
    if !isfile(fn) && isfile(tempFn)
        mv(tempFn, fn)
    end

    # Delete any left over temp files
    # Program was terminated while config was saved
    if isfile(tempFn)
        rm(tempFn)
    end

    if isfile(fn)
        return Config(fn, open(JSON.parse, fn))

    else
        config = Config(fn, default)
        saveConfig(config)

        return config
    end
end

setDefaults!(c::Config, d::Dict{K, V}) where {K, V} = c.data = merge(d, c.data)

Base.haskey(c::Config, key::String) = haskey(c.data, key)

function Base.setindex!(c::Config, value::V, key::K) where {K, V}
    prev = haskey(c, key) ? c.data[key] : nothing
    c.data[key] = value

    if value != prev || prev === nothing
        saveConfig(c, false)
    end
end

function Base.getindex(c::Config, key::Any)
    now = time()

    if now > c.lastCheck + c.checkDelay
        c.lastCheck = now

        if mtime(c.fn) > c.mtime
            c.mtime = mtime(c.fn)
            c.data = open(JSON.parse, c.fn)
        end
    end

    return c.data[key]
end

function Base.get(c::Config, key::K, value::V) where {K, V}
    now = time()

    if now > c.lastCheck + c.checkDelay
        c.lastCheck = now

        if mtime(c.fn) > c.mtime
            c.mtime = mtime(c.fn)
            c.data = open(JSON.parse, c.fn)
        end
    end

    # Force a rewrite of the config
    if !haskey(c.data, key)
        c[key] = value
    end

    return c.data[key]
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
# Optimization and solution query

"""
Call Knitro to solve the problem.  The return value indicates
the solution status.

"""
function KN_solve(m::Model)
    ret = @kn_ccall(solve, Cint, (Ptr{Cvoid},), m.env)
    # For KN_solve, we do not return an error if ret is different of 0
    return ret
end

##################################################
# Solution status and info
##################################################
"Return the solution status, objective, primal and dual variables."
function KN_get_solution(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    nx = KN_get_number_vars(m)
    nc = KN_get_number_cons(m)

    x = zeros(Cdouble, nx)
    lambda = zeros(Cdouble, nx + nc)
    status = Cint[0]
    obj = Cdouble[0.]

    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, x, lambda)

    return status[1], obj[1], x, lambda
end

# some wrapper functions for MOI
function get_status(m::Model)
    @assert m.env != C_NULL
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, C_NULL, C_NULL)
    return status[1]
end

function get_objective(m::Model)
    @assert m.env != C_NULL
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, C_NULL, C_NULL)
    return obj[1]
end

function get_solution(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    nx = KN_get_number_vars(m)
    x = zeros(Cdouble, nx)
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, x, C_NULL)
    return x
end

function get_dual(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    nx = KN_get_number_vars(m)
    nc = KN_get_number_cons(m)
    lambda = zeros(Cdouble, nx + nc)
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, C_NULL, lambda)
    return lambda
end


text:
import Base: rand,
             ∈,
             isempty

export EmptySet, ∅,
       an_element,
       linear_map

"""
    EmptySet{N} <: LazySet{N}

Type that represents the empty set, i.e., the set with no elements.
"""
struct EmptySet{N} <: LazySet{N}
    dim::Int
end

isoperationtype(::Type{<:EmptySet}) = false
isconvextype(::Type{<:EmptySet}) = true

# default constructor of type Float64
EmptySet(dim::Int) = EmptySet{Float64}(dim)

"""
    ∅

Alias for `EmptySet{Float64}`.
"""
const ∅ = EmptySet{Float64}

# --- LazySet interface functions ---

"""
    dim(∅::EmptySet)

Return the dimension of an empty set.

### Input

- `∅` -- an empty set

### Output

The dimension of the empty set.
"""
function dim(∅::EmptySet)
    return ∅.dim
end

"""
    σ(d::AbstractVector, ∅::EmptySet)

Return the support vector of an empty set.

### Input

- `d` -- direction
- `∅` -- empty set

### Output

An error.
"""
function σ(d::AbstractVector, ∅::EmptySet)
    error("the support vector of an empty set does not exist")
end

"""
    ρ(d::AbstractVector, ∅::EmptySet)

Evaluate the support function of an empty set in a given direction.

### Input

- `d` -- direction
- `∅` -- empty set

### Output

An error.
"""
function ρ(d::AbstractVector, ∅::EmptySet)
    error("the support function of an empty set does not exist")
end

function isboundedtype(::Type{<:EmptySet})
    return true
end

"""
    isbounded(∅::EmptySet)

Determine whether an empty set is bounded.

### Input

- `∅` -- empty set

### Output

`true`.
"""
function isbounded(::EmptySet)
    return true
end

"""
    isuniversal(∅::EmptySet{N}, [witness]::Bool=false) where {N}

Check whether an empty is universal.

### Input

- `∅`       -- empty set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `false`
* If `witness` option is activated: `(false, v)` where ``v ∉ S``, although
  we currently throw an error
"""
function isuniversal(∅::EmptySet{N}, witness::Bool=false) where {N}
    if witness
        return (false, zeros(N, dim(∅)))
    else
        return false
    end
end

"""
    ∈(x::AbstractVector, ∅::EmptySet)

Check whether a given point is contained in an empty set.

### Input

- `x` -- point/vector
- `∅` -- empty set

### Output

The output is always `false`.

### Examples

```jldoctest
julia> [1.0, 0.0] ∈ ∅(2)
false
```
"""
function ∈(x::AbstractVector, ∅::EmptySet)
    return false
end

"""
    an_element(∅::EmptySet)

Return some element of an empty set.

### Input

- `∅` -- empty set

### Output

An error.
"""
function an_element(∅::EmptySet)
    error("an empty set does not have any element")
end

"""
    rand(::Type{EmptySet}; [N]::Type{<:Real}=Float64, [dim]::Int=0,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create an empty set (note that there is nothing to randomize).

### Input

- `EmptySet` -- type for dispatch
- `N`        -- (optional, default: `Float64`) numeric type
- `dim`      -- (optional, default: 2) dimension
- `rng`      -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`     -- (optional, default: `nothing`) seed for reseeding

### Output

The (only) empty set of the given numeric type and dimension.
"""
function rand(::Type{EmptySet};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    return EmptySet{N}(dim)
end

"""
    isempty(∅::EmptySet)

Return if the empty set is empty or not.

### Input

- `∅` -- empty set

### Output

`true`.
"""
function isempty(∅::EmptySet)
    return true
end

"""
    norm(S::EmptySet, [p]::Real=Inf)

Return the norm of an empty set.
It is the norm of the enclosing ball (of the given ``p``-norm) of minimal volume
that is centered in the origin.

### Input

- `S` -- empty set
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function norm(S::EmptySet, p::Real=Inf)
    error("an empty set does not have a norm")
end

"""
    radius(S::EmptySet, [p]::Real=Inf)

Return the radius of an empty set.
It is the radius of the enclosing ball (of the given ``p``-norm) of minimal
volume with the same center.

### Input

- `S` -- empty set
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function radius(S::EmptySet, p::Real=Inf)
    error("an empty set does not have a radius")
end

"""
    diameter(S::EmptySet, [p]::Real=Inf)

Return the diameter of an empty set.
It is the maximum distance between any two elements of the set, or,
equivalently, the diameter of the enclosing ball (of the given ``p``-norm) of
minimal volume with the same center.

### Input

- `S` -- empty set
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function diameter(S::EmptySet, p::Real=Inf)
    error("an empty set does not have a diameter")
end

"""
    vertices(∅::EmptySet{N}) where {N}

Construct an iterator over the vertices of an empty set.

### Input

- `∅` -- empty set

### Output

The empty iterator, as the empty set does not contain any vertices.
"""
function vertices(∅::EmptySet{N}) where {N}
    return EmptyIterator{Vector{N}}()
end

"""
    vertices_list(∅::EmptySet{N}) where {N}

Return the list of vertices of an empty set.

### Input

- `∅` -- empty set

### Output

The empty list of vertices, as the empty set does not contain any vertices.
"""
function vertices_list(∅::EmptySet{N}) where {N}
    return Vector{Vector{N}}()
end

"""
    linear_map(M::AbstractMatrix{N}, ∅::EmptySet{N}) where {N}

Return the linear map of an empty set.

### Input

- `M` -- matrix
- `∅` -- empty set

### Output

The empty set.
"""
linear_map(M::AbstractMatrix{N}, ∅::EmptySet{N}) where {N} = ∅

"""
    translate(∅::EmptySet, v::AbstractVector)

Translate (i.e., shift) an empty set by a given vector.

### Input

- `∅` -- empty set
- `v` -- translation vector

### Output

The empty set.
"""
function translate(∅::EmptySet, v::AbstractVector)
    return translate!(∅, v)
end

function translate!(∅::EmptySet, v::AbstractVector)
    @assert length(v) == dim(∅) "cannot translate a $(dim(∅))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return ∅
end

"""
    plot_recipe(∅::EmptySet{N}, [ε]=zero(N)) where {N}

Convert an empty set to a sequence of points for plotting.
In the special case of an empty set, we define the sequence as `nothing`.

### Input

- `∅` -- empty set
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

`nothing`.
"""
function plot_recipe(∅::EmptySet{N}, ε=zero(N)) where {N}
    return []
end

"""
    area(∅::EmptySet{N}) where {N}

Return the area of an empty set.

### Input

- `∅` -- empty set

### Output

The zero element of type `N`.
"""
function area(∅::EmptySet{N}) where {N}
    return zero(N)
end

function project(∅::EmptySet{N}, block::AbstractVector{Int}; kwargs...) where {N}
    return EmptySet{N}(length(block))
end

function volume(∅::EmptySet{N}) where {N}
    return zero(N)
end


text:
#==
Logistic
    Regression
==#
@doc """
      Majority class baseline is used to find the most often interpreted
      classification in an array.\n
      --------------------\n
      ==PARAMETERS==\n
     [X] <- Fill with your trainX values. Should be an array of shape (0,1) or (1,0)\n
     [y] <- Fill with your trainy values. Should be an array of shape (0,1) or (1,0)\n
     λ = .0001 <- Lambda Value\n
     fit_intercept = true <- Boolean determines whether to fit an intercept.\n
     max_iter = 1000 <- Determines the maximum number of iterations for the model to perform.\n
     --------------------\n
     ==Functions==\n
     predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)\n
       """
function LogisticRegression(X, y, λ=0.0001, fit_intercept=true, η=0.01, max_iter=1000)
    θ, 𝐉 = logistic_regression_sgd(X, y, 0.0001, true, 0.3, 3000);
    predict(xt) = yhat = predict_class(predict_proba(xt,0))
    cost = 𝐉
    (var) -> (predict;cost)
end
function sigmoid(z)
    return 1 ./ (1 .+ exp.(.-z))
end
function logistic_regression_sgd(X, y, λ, fit_intercept=true, η=0.01, max_iter=1000)

    # Initialize some useful values
    m = length(y); # number of training examples

    if fit_intercept
        # Add a constant of 1s if fit_intercept is specified
        constant = ones(m, 1)
        X = hcat(constant, X)
    else
        X # Assume user added constants
    end

    # Use the number of features to initialise the theta θ vector
    n = size(X)[2]
    θ = zeros(n)

    # Initialise the cost vector based on the number of iterations
    𝐉 = zeros(max_iter)

    for iter in range(1, stop=max_iter)

        # Calcaluate the cost and gradient (∇𝐉) for each iter
        𝐉[iter], ∇𝐉 = regularised_cost(X, y, θ, λ)

        # Update θ using gradients (∇𝐉) for direction and (η) for the magnitude of steps in that direction
        θ = θ - (η * ∇𝐉)
    end

    return (θ, 𝐉)
end
function regularised_cost(X, y, θ, λ)
    m = length(y)
    h = sigmoid(X * θ)
    positive_class_cost = ((-y)' * log.(h))
    negative_class_cost = ((1 .- y)' * log.(1 .- h))
    lambda_regularization = (λ/(2*m) * sum(θ[2 : end] .^ 2))
    𝐉 = (1/m) * (positive_class_cost - negative_class_cost) + lambda_regularization
    ∇𝐉 = (1/m) * (X') * (h-y) + ((1/m) * (λ * θ))
    ∇𝐉[1] = (1/m) * (X[:, 1])' * (h-y)
           return (𝐉, ∇𝐉)
end
function predict_proba(X, θ, fit_intercept=true)
    m = size(X)[1]

    if fit_intercept
        constant = ones(m, 1)
        X = hcat(constant, X)
    else
        X
    end

    h = sigmoid(X * θ)
    return h
end

function predict_class(proba, threshold=0.5)
    return proba .>= threshold
end


text:
@testset "7.3 Inverse hyperbolic tangent" begin
    include("7.3.2 (d x)^m (a+b arctanh(c x^n))^p.jl")
    include("7.3.3 (d+e x)^m (a+b arctanh(c x^n))^p.jl")
    include("7.3.4 u (a+b arctanh(c x))^p.jl")
    include("7.3.5 u (a+b arctanh(c+d x))^p.jl")
    include("7.3.6 Exponentials of inverse hyperbolic tangent functions.jl")
    include("7.3.7 Inverse hyperbolic tangent functions.jl")
end

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/.jl