text:
module HugepageMmap

using Mmap

export hugepage_mmap,
       PageSize4K,
       PageSize2M,
       PageSize1G

const MAP_HUGETLB    = Cint(0x40000)
const MAP_HUGE_SHIFT = Cint(26)
const MAP_HUGE_2MB   = Cint(21 << MAP_HUGE_SHIFT)
const MAP_HUGE_1GB   = Cint(30 << MAP_HUGE_SHIFT)

abstract type AbstractPageSize end
struct PageSize4K <: AbstractPageSize end
struct PageSize2M <: AbstractPageSize end
struct PageSize1G <: AbstractPageSize end

extraflags(::PageSize4K) = Cint(0)
extraflags(::PageSize2M) = MAP_HUGETLB | MAP_HUGE_2MB
extraflags(::PageSize1G) = MAP_HUGETLB | MAP_HUGE_1GB

# Align length for `munmap` to a multiple of page size.
pagesize(::PageSize4K) = 4096
pagesize(::PageSize2M) = 2097152
pagesize(::PageSize1G) = 1073741824

align(x, p::AbstractPageSize) = ceil(Int, x / pagesize(p)) * pagesize(p)
align(x, ::PageSize4K) = x

# This is heavily based on the Mmap stdlib
function hugepage_mmap(::Type{T}, pagesize::AbstractPageSize, dims::Integer...) where {T}
    return hugepage_mmap(T, pagesize, convert.(Int, dims))
end

"""
    hugepage_mmap(::Type{T}, pagesize::AbstractPageSize, dims...)

Allocate a `Array{T}` with dimensions `dims` backed by `pagesize`.
Choices for `pagesize` are `PageSize4k()`, `PageSize2m()`, or `PageSize1G()`.
Make sure the corresponding hugepages are already allocated on your system.
"""
function hugepage_mmap(::Type{T}, pagesize::AbstractPageSize, dims::NTuple{N,Int}) where {N,T}
    mmaplen = sizeof(T) * prod(dims)

    # Build the PROT flags - we want to be able to read and write.
    prot = Mmap.PROT_READ | Mmap.PROT_WRITE
    flags = Mmap.MAP_PRIVATE | Mmap.MAP_ANONYMOUS
    flags |= extraflags(pagesize)

    fd = Base.INVALID_OS_HANDLE
    offset = Cint(0)

    # Fordward this call into the Julia C library.
    ptr = ccall(
        :jl_mmap,
        Ptr{Cvoid},
        (Ptr{Cvoid}, Csize_t, Cint, Cint, RawFD, Int64),
        C_NULL,     # No address we really want.
        mmaplen,
        prot,
        flags,
        fd,
        offset,
    )

    # Wrap this into an Array and attach a finalizer that will unmap the underlying pointer
    # when the Array if GC'd
    A = Base.unsafe_wrap(Array, Ptr{T}(ptr), dims)
    finalizer(A) do x
        systemerror(
            "munmap",
            ccall(:munmap, Cint, (Ptr{Cvoid}, Csize_t), ptr, align(mmaplen, pagesize)) != 0
        )
    end
    return A
end

end # module


text:
using Test, GenericLinearAlgebra, LinearAlgebra
using GenericLinearAlgebra.LAPACK2

@testset "LAPACK wrappers" begin

    n = 100

    T = SymTridiagonal(randn(n), randn(n - 1))
    vals, vecs = eigen(T)
    @testset "steqr" begin
        _vals, _vecs = LAPACK2.steqr!('N', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals

        _vals, _vecs = LAPACK2.steqr!('I', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)

        _vals, _vecs = LAPACK2.steqr!('V', copy(T.dv), copy(T.ev), Matrix{Float64}(I, n, n))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)
    end

    @testset "sterf" begin
        _vals = LAPACK2.sterf!(copy(T.dv), copy(T.ev))
        @test _vals ≈ vals
    end

    @testset "stedc" begin
        _vals, _vecs = LAPACK2.stedc!('N', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals

        _vals, _vecs = LAPACK2.stedc!('I', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)

        _vals, _vecs = LAPACK2.stedc!('V', copy(T.dv), copy(T.ev), Matrix{Float64}(I, n,n))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)
    end

    @testset "stemr" begin
        _vals, _vecs, __ = LAPACK2.stemr!('N', 'A', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals

        _vals, _vecs, __ = LAPACK2.stemr!('N', 'V', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals

        _vals, _vecs, __ = LAPACK2.stemr!('N', 'I', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals

        _vals, _vecs, __ = LAPACK2.stemr!('V', 'A', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)

        _vals, _vecs, __ = LAPACK2.stemr!('V', 'V', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)

        _vals, _vecs, __ = LAPACK2.stemr!('V', 'I', copy(T.dv), copy(T.ev))
        @test vals ≈ _vals
        @test abs.(_vecs'vecs) ≈ Matrix(I, n, n)
    end

    @testset "lahqr" begin
        T = Array(Tridiagonal(fill(0.1, 99), fill(0.0, 100), fill(40.0, 99)))
        _vals = sort(LAPACK2.lahqr!(copy(T))[1])
        @test _vals ≈ sort(real.(GenericLinearAlgebra._eigvals!(copy(T))))
        # LAPACK's multishift algorithm (the default) seems to be broken
        @test !(_vals ≈ sort(eigvals(T)))
    end
end

ext: jl


text:
using Graphs
using DataStructures
#to use connected components, must implement: vertex_list,vertex_map,adjacency_list

type road #start coordinates are always closer to (1,1) than end coordinates
  start_x::Int64
  start_y::Int64
  end_x::Int64
  end_y::Int64
  net_id::Int64
end

type block #upper left corner is always closer to (1,1) than any other corner
  x::Int64
  y::Int64
  dx::Int8
  dy::Int8
end

function gen_test()
  hold = Deque{road}()
  push!(hold,road(5,3,5,3,rand(Int64)))
  push!(hold,road(11,1,11,6,rand(Int64)))
  push!(hold,road(10,3,13,3,rand(Int64)))
  push!(hold,road(5,3,10,3,rand(Int64)))
  push!(hold,road(11,6,11,20,rand(Int64)))
  push!(hold,road(1,3,5,3,rand(Int64)))
  push!(hold,road(20,20,25,20,rand(Int64)))
  push!(hold,road(25,20,25,30,rand(Int64)))
  return hold
end

function intersects(a::road,b::road) #requires that 2 roads overlap
  if (a.start_x <= b.end_x && b.start_x <= a.start_x)||(b.start_x <= a.end_x && a.start_x <= b.start_x)
    if  (a.start_y <= b.end_y && b.start_y <= a.start_y)||(b.start_y <= a.end_y && a.start_y <= b.start_y)
      return true
    else
      return false
    end
  end
  return false
end

function islonely(a::road)
  if (a.start_y == a.end_y) && (a.start_x == a.end_x)
    return true
  end
  return false
end

function ishorizontal(a::road) #returns true if either a is horizontal or a is one tile
  if a.start_y == a.end_y || (a.start_y ==a.end_y && a.start_x == a.end_x)
    return true
  end
  return false
end

function length(road)
  return sqrt((road.start_x-road.end_x)^2+(road.start_y-road.end_y)^2)
end

function merge(a::road,b::road)
  return road(minimum([a.start_x,b.start_x]),minimum([a.start_y,b.start_y]),
              maximum([a.end_x,b.end_x]),maximum([a.end_y,b.end_y]),rand(Int64))
end

function merge_roads!(roads::Deque{road}) #empties roads of all elements
  fixed = Deque{road}()
  hold = Deque{road}()
  while !isempty(hold) || !isempty(roads)
    while !isempty(hold)
      push!(roads,shift!(hold))
    end
    x = shift!(roads)
    merged = false
    while  !isempty(roads)
      y = shift!(roads)
      if intersects(x,y) && (islonely(x)||islonely(y)||(ishorizontal(x) == ishorizontal(y)))
        x = merge(x,y)
        merged = true
      else
        push!(hold,y)
      end
    end
    if !merged
      push!(fixed,x)
    else
      push!(roads,x)
    end
  end
  return fixed
end

function make_graph!(roads::Deque{road}) #assume roads has already been 'fixed'
  hold = Deque{KeyVertex{road}}()
  g = adjlist(KeyVertex{road},is_directed=false)
  while !isempty(roads)
    x = shift!(roads)
    push!(hold,add_vertex!(g,x))
  end
  done = Deque{KeyVertex{road}}()
  osci = Deque{KeyVertex{road}}()
  from_hold = true
  while !isempty(hold) || !isempty(osci)
    if from_hold
      x = shift!(hold)
      while !isempty(hold)
        y = shift!(hold)
        if intersects(x.key,y.key)
          add_edge!(g,x,y)
        end
        push!(osci,y)
      end
      from_hold = false
    else
      x = shift!(osci)
      while !isempty(osci)
        y = shift!(osci)
        if intersects(x.key,y.key)
          add_edge!(g,x,y)
        end
        push!(hold,y)
      end
      from_hold = true
    end
  end
  return g
end

function set_net_id!(connected_components)
  for i in connected_components
    x = rand(Int64)
    for j in i
      j.key.net_id = x
    end
  end
end

max_stars_repo_name: JKSano/netgameproj

max_issues_repo_name: JKSano/netgameproj

max_forks_repo_name: JKSano/netgameproj


text:
module GitHubActionsUtils

import GitHub

function get_env(x)
    get(ENV, x, nothing)
end

is_github_actions() = get_env("CI") == "true" && get_env("GITHUB_ACTIONS") == "true"
event_name() = something(get_env("GITHUB_EVENT_NAME"))
is_push() = event_name() == "push"
is_pull_request() = event_name() == "pull_request"
head_ref() = something(get_env("GITHUB_HEAD_REF"))
github_ref() = something(get_env("GITHUB_REF"))
is_branch() = match(r"^refs/heads/(.*)", github_ref()) !== nothing
branch_name() = is_branch() ? match(r"^refs/heads/(.*)", github_ref())[1] : nothing
is_tag() = match(r"^refs/tags/(.*)", github_ref()) !== nothing
tag_name() = is_tag() ? match(r"^refs/tags/(.*)", github_ref())[1] : nothing
pull_request_number() = is_pull_request() ? parse(Int, match(r"^refs/pull/(\d+)/.*", github_ref())[1]) : nothing
pull_request_source() = is_pull_request() ? head_ref() : nothing
repository() = something(get_env("GITHUB_REPOSITORY"))

const _auth = Ref{Any}()

function auth()
    if !isassigned(_auth)
        _auth[] = GitHub.authenticate(ENV["GITHUB_TOKEN"])
    end
    return _auth[]
end

function comment_on_pr(pr_id, comment)
    GitHub.create_comment(repository(), pr_id; auth = auth(), params = :body => comment)
end

is_local_git_branch(branch) = success(`git show-ref refs/heads/$branch`)

function switch_to_or_create_branch(branch; orphan = false)
    if is_local_git_branch(branch)
        run(`git switch $branch`)
    else
        if orphan
            run(`git switch --orphan $branch`)
        else
            run(`git switch -c $branch`)
        end
    end
end

push_git_branch(branch; remote = "origin") = run(`git push $remote $branch`)

function set_github_actions_bot_as_git_user()
    bot_address = "41898282+github-actions[bot]@users.noreply.github.com"
    run(`git config --local user.email $bot_address`)
    run(`git config --local user.name "github-actions"`)
end

end


text:
#' ---
#' title: Sudoku
#' ---

#' **Originally Contributed by**: Iain Dunning

#'
#'
A partially solved Sudoku puzzle

#'
#' [Sudoku](http://en.wikipedia.org/wiki/Sudoku) is a popular number puzzle. The goal is to place the digits 1,...,9 on a nine-by-nine grid, with some of the digits already filled in. Your solution must satisfy the following rules:
#'
#' * The numbers 1 to 9 must appear in each 3x3 square
#' * The numbers 1 to 9 must appear in each row
#' * The numbers 1 to 9 must appear in each column
#'
#' This isn't an optimization problem, its actually a *feasibility* problem: we wish to find a feasible solution that satsifies these rules. You can think of it as an optimization problem with an objective of 0.
#'
#' We can model this problem using 0-1 integer programming: a problem where all the decision variables are binary. We'll use JuMP to create the model, and then we can solve it with any integer programming solver.

using JuMP
using GLPK

#' We will define a binary variable (a variable that is either 0 or 1) for each possible number in each possible cell. The meaning of each variable is as follows:
#'
#'     x[i,j,k] = 1  if and only if cell (i,j) has number k
#'
#' where `i` is the row and `j` is the column.

# Create a model
sudoku = Model(with_optimizer(GLPK.Optimizer))

# Create our variables
@variable(sudoku, x[i=1:9, j=1:9, k=1:9], Bin);

#' Now we can begin to add our constraints. We'll actually start with something obvious to us as humans, but what we need to enforce: that there can be only one number per cell.

for i = 1:9, j = 1:9  # Each row and each column
    # Sum across all the possible digits
    # One and only one of the digits can be in this cell,
    # so the sum must be equal to one
    @constraint(sudoku, sum(x[i,j,k] for k in 1:9) == 1)
end

#' Next we'll add the constraints for the rows and the columns. These constraints are all very similar, so much so that we can actually add them at the same time.

for ind = 1:9  # Each row, OR each column
    for k = 1:9  # Each digit
        # Sum across columns (j) - row constraint
        @constraint(sudoku, sum(x[ind,j,k] for j in 1:9) == 1)
        # Sum across rows (i) - column constraint
            @constraint(sudoku, sum(x[i,ind,k] for i in 1:9) == 1)
    end
end

#' Finally, we have the to enforce the constraint that each digit appears once in each of the nine 3x3 sub-grids. Our strategy will be to index over the top-left corners of each 3x3 square with `for` loops, then sum over the squares.

for i = 1:3:7, j = 1:3:7, k = 1:9
    # i is the top left row, j is the top left column
    # We'll sum from i to i+2, e.g. i=4, r=4, 5, 6
    @constraint(sudoku, sum(x[r,c,k] for  r in i:i+2, c in j:j+2) == 1)
end

#' The final step is to add the initial solution as a set of constraints. We'll solve the problem that is in the picture at the start of the notebook. We'll put a `0` if there is no digit in that location.

# The given digits
init_sol = [ 5 3 0 0 7 0 0 0 0;
             6 0 0 1 9 5 0 0 0;
             0 9 8 0 0 0 0 6 0;
             8 0 0 0 6 0 0 0 3;
             4 0 0 8 0 3 0 0 1;
             7 0 0 0 2 0 0 0 6;
             0 6 0 0 0 0 2 8 0;
             0 0 0 4 1 9 0 0 5;
             0 0 0 0 8 0 0 7 9]
for i = 1:9, j = 1:9
    # If the space isn't empty
    if init_sol[i,j] != 0
        # Then the corresponding variable for that digit
        # and location must be 1
        @constraint(sudoku, x[i,j,init_sol[i,j]] == 1)
    end
end

# solve problem
optimize!(sudoku)

# Extract the values of x
x_val = value.(x)
# Create a matrix to store the solution
sol = zeros(Int,9,9)  # 9x9 matrix of integers
for i in 1:9, j in 1:9, k in 1:9
    # Integer programs are solved as a series of linear programs
    # so the values might not be precisely 0 and 1. We can just
    # round them to the nearest integer to make it easier
    if round(Int,x_val[i,j,k]) == 1
        sol[i,j] = k
    end
end
# Display the solution
sol

#' Which is the correct solution:
#'


text:
@doc raw"""
    simple_operation(num1::Float64, num2::Float64)

Function for adding two numbers
"""
function simple_operation(num1::Float64, num2::Float64)
    num3 = num1 + num2
    return num3
end


text:
mutable struct SingleLayer <: Architecture
  wts::AbstractVecOrMat
end

mutable struct MultiLayer <: Architecture
  wts::AbstractMatrix
end

# mutable struct ANN <: Architecture
#
# end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1ba17e2e-a9d2-11eb-3e83-bd0873b44b78
md"# Data Structures and Algorithms Final Exam Review

- CIS 3223 Section 002 (Spring 2021)

- Instructor: Yu Wang

- Temple University"

# ╔═╡ ed49dc2f-7845-456d-9c41-481996339d98
md"
### Format

- Close-book and close-notes
    - Allow 2 sheets (4 sides) of notes
    - Submit your notes to get credits
- In class, via Zoom + Canvas
    - Join **regular zoom meeting** for announcement & questions
    - Sign out after submission
- The final will be delivered as an assignment on canvas
    - Upload your solutions via the assignment in Canvas
- Please work on your own, do NOT discuss anything with your classmates
    - No messaging or group chatting allowed
"

# ╔═╡ fe320c90-5d5d-410f-8715-c1ef4c8ab9fb
md"
### Time & Submission

- Monday (May 3) 3:30-5:50pm
    - 120 minutes for exam and the remaining for downloading & uploading
- Make sure you have stable internet during the exam period
"

# ╔═╡ 92ce6a5e-0ad3-4905-9d40-f8048acc10af
md"
### Number of Problems and Grade

- 9 Problems total.
- Full mark is 120.
- Points will be deducted for difficult to read answers.
- 30% of your final score comes from your final exam.
"

# ╔═╡ ca1cc38c-da0e-4d64-8fe2-c195f46f88cc
md"
### Style: Quiz, Midterm, homework

- Calculation or run algorithm on examples
    - e.g. running DFS, Dijkstra, MST, Max-flow, calculating RSA or gcd(x,y)
- Design or analyze algorithms
    - e.g. design/analyze divide-and-conquer, greedy, DP algorithms
- Prove certain claims/lemmas or analyze certain functions
    - e.g. proving Big-O relation, proving using *induction* or *contradiction*, solving recursive function using Master theorem.
"

# ╔═╡ d8a8ee9c-b782-4413-b665-b344332216f2
md"
### Review material (Ch0-Ch7):

- Homework #1-#8 + Quiz #1-#3 + Midterm
- Lecture Notes
- Textbook
"

# ╔═╡ a6bbe6a3-145b-498f-ac53-bb9bdf2cd72b
md"## Important Topics

### 0 Prologue

- big-O notation
- enter Fibonacci

### 1 Algorithms with numbers

- modular arithmetic
- cryptography

### 2 Divide and conquer

- mergesort
- master theorem

### 3 Decompositions of graphs

- Depth-first search in undirected graphs
- Depth-first search in directed graphs
- Strongly connected components

### 5 Paths in graphs

- breadth-first search
- Dijkstsra's algorithm

### 5 Greedy algorithms

- minimum spanning trees
- Kruskal, Prim

### 6 Dynamic programming

- weighted interval scheduling
- knapsack
- independent sets in trees

### 7.2 Flows in networks

- max-flow and min-cut problems
- Ford-Fulkerson algorithm"

# ╔═╡ 4b240b0f-aeb2-43d4-8957-acc77bc1b044
md"# Algorithms"

# ╔═╡ 22507365-4ef7-4400-a9ba-7e93dc4a5203
# Exponential fibonacci algorithm
function fib1(n)
	if n == 0 return 0 end
	if n == 1 return 1 end
	return fib1(n - 1) + fib1(n - 2)
end;

# ╔═╡ c51f5256-a278-4ce3-b0cf-7274d174f454
# Polynomial fibonacci algorithm
function fib2(n)
	if n == 0 return 0 end
	f = zeros(Int, n+1)
	f[2] = 1
	for i = 3:n
		f[i] = f[i - 1] + f[i - 2]
	end
	return f[n]
end;

# ╔═╡ 2174fcff-76a9-47ab-99fd-662679baf33f
# Decimal to binary
function dec2bin(N)
	s = ""
	while N > 0
		s = "$(N % 2)$s"
		N >>= 1
	end
	return s
end;

# ╔═╡ 7b96f661-a77b-4a40-8269-d14ec887a46a
dec2bin(245)

# ╔═╡ 9db70690-4f7b-49a1-b259-a490c0f6787a
# Decimal to hexadecimal
function dec2hex(N)
	s = ""
	while N > 0
		d = let
			r = N % 16
			if r > 9
				'a' + r - 10
			else
				r
			end
		end
		s = "$d$s"
		N >>= 4
	end
	return s
end;

# ╔═╡ dc49581c-2e2d-48c2-85f4-3aeeddbae8f0
dec2hex(249231)

# ╔═╡ 103dcea9-522b-45e0-bbc4-851450a755aa
# Convert decimal to any base b < 10
function dec2base(N, b)
	s = ""
	while N > 0
		s = "$(N % b)$s"
		N = N ÷ b
	end
	return s
end;

# ╔═╡ dfdf2ab2-797d-4744-a2ee-379320e120a5
dec2base(245, 3)

# ╔═╡ 0c60b870-e122-4592-b42f-d5b64941d908
function dec2base_steps(N, b)
	s = [b N -1]
	while N > 0
		(q, r) = divrem(N, b)
		s = [s; [b q r]]
		N = N ÷ b
	end
	return s, join(reverse(s[2:end,3]))
end;

# ╔═╡ 37c60ab8-f2ef-41c3-8fdd-ecbf6189d625
dec2base_steps(245, 3)

# ╔═╡ 1754fc92-1036-4ce6-b209-f4f5bc75efba
# Euclid's gcd algorithm
function Euclid(a, b)
	if b == 0 return a end
	return Euclid(b, a % b)
end;

# ╔═╡ 4bc44e0c-b2e1-4866-86df-da59331a620e
# Extended Euclidean algorithm
function extEuclid(a, b)
	if b == 0 return (1, 0, a) end
	(x′, y′, d) = extEuclid(b, a % b)
	return (y′, x′ - a ÷ b * y′, d)
end

# ╔═╡ c311334c-eee4-4171-93e3-56aa6e822c9f
extEuclid(210, 588)

# ╔═╡ 5481ab40-7e2c-4f35-94a1-90df39e2558f
let
	# Prove that T(n) = n log₂n, assuming n is a power of 2
	function T(n)
		if n == 1
			return 0
		else
			return 2T(n ÷ 2) + n
		end
	end

	T(2^11), (n -> n * log2(n))(2^11)
end

# ╔═╡ da5b0b05-dc8c-4bca-8b8c-c50b89f0018c
md"# Homework Problems"

# ╔═╡ d097fbac-2ab5-4e44-b9cc-98a694c2a7fa
md"## 0 Prologue"

# ╔═╡ 3918986d-8fcf-41ff-91bc-52e5eb6346c0
md"### Problem 0.1

In each of the following situations, indicate whether $f = O(g)$, or $f = \Omega(g)$, or both (in which case $f = \Theta(g)$).

$
\begin{array}{cll}
   & f(n) & g(n) \\
	\mathbf{(a)} & n - 100 & n - 200 \\
	\mathbf{(b)} & n^{1/2} & n^{2/3} \\
	\mathbf{(c)} & 100n + \log{n} & n + (\log{n})^2 \\
	\mathbf{(d)} & n \log{n} & 10n \log{10n} \\
	\mathbf{(e)} & \log{2n} & \log{3n} \\
	\mathbf{(f)} & 10 \log{n} & \log{(n^2)} \\
	\mathbf{(g)} & n^{1.01} & n \log^2{n} \\
	\mathbf{(h)} & n^2 / \log{n} & n (\log{n})^2 \\
	\mathbf{(i)} & n^{0.1} & (\log{n})^{10} \\
	\mathbf{(j)} & (\log{n})^{\log{n}} & n / \log{n} \\
	\mathbf{(k)} & \sqrt{n} & (\log{n})^3 \\
	\mathbf{(l)} & n^{1/2} & 5^{\log_2{n}} \\
	\mathbf{(m)} & n2^n & 3^n \\
	\mathbf{(n)} & 2^n & 2^{n+1} \\
	\mathbf{(o)} & n! & 2^n \\
	\mathbf{(p)} & (\log{n})^{\log{n}} & 2^{(\log_2{n})^2} \\
	\mathbf{(q)} & \sum_{i=1}^n i^k & n^{k+1}
\end{array}$
"

# ╔═╡ f471306e-733c-421f-a243-1a0079989dfb
md"#### Part (a)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n - 100}{n - 200} = \lim_{n \to \infty} \frac{1 - \frac{100}{n}}{1 - \frac{200}{n}} = 1 \implies f = \Theta(g)$"

# ╔═╡ 33a9b77b-2c02-40de-8072-1164b6614e25
md"#### Part (b)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n^{1/2}}{n^{2/3}} = \frac{n^{1/2 - 2/3}}{1} = \lim_{n \to \infty} \frac{1}{n^{1/3}} = 0 \implies f = O(g)$"

# ╔═╡ be37a3a8-0a69-4b9c-93eb-d4025b2fccea
md"#### Part (c)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{100n + \log{n}}{n + (\log{n})^2} = \lim_{n \to \infty} \frac{100 + \frac{\log{n}}{n}}{1 + \frac{(\log{n})^2}{n}} = 100 \implies f = \Theta(g)$"

# ╔═╡ 4ee7fe91-2397-45af-865a-c252920bcf59
md"#### Part (d)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n \log{n}}{10n \log{10n}} = \lim_{n \to \infty} \frac{1}{\frac{10n\log{10}}{n \log{n}} + \frac{10n\log{n}}{n \log{n}}} = \frac{1}{10} \implies f = \Theta(g)$
"

# ╔═╡ 02641c72-1bd8-4855-b169-1e67b7f840c8
md"#### Part (e)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{\log{2n}}{\log{3n}} = \lim_{n \to \infty} \frac{\log{2} + \log{n}}{\log{3} + \log{n}} = \lim_{n \to \infty} \frac{\frac{\log{2}}{\log{n}} + 1}{\frac{\log{3}}{\log{n}} + 1} = 1 \implies f = \Theta(g)$
"

# ╔═╡ e09e3174-4b62-4d48-8640-202d1da7fb3f
md"#### Part (f)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{10 \log{n}}{\log{(n^2)}} = \lim_{n \to \infty} \frac{10 \log{n}}{2 \log{n}} = 5 \implies f = \Theta(g)$"

# ╔═╡ 510a5fb5-e96b-44f7-9257-40fbfd4ea191
md"#### Part (g)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n^{1.01}}{n \log^2{n}} = \lim_{n \to \infty} \frac{n^{0.01}}{\log^2{n}} = \infty \implies f = \Omega(g)$"

# ╔═╡ 5e4d7f29-974b-4d25-aa7b-8bf827ab1589
md"#### Part (h)

$\lim_{n \to \infty} \frac{\frac{n^2}{\log{n}}}{n (\log{n})^2} = \lim_{n \to \infty} \frac{n}{(\log{n})^3} = \infty \implies f = \Omega(g)$"

# ╔═╡ 5272ccd3-959d-4d29-8626-3d83cdb830a7
md"#### Part (i)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n^{0.1}}{(\log{n})^{10}} = \infty \implies f = \Omega(g)$"

# ╔═╡ 93046ae6-acf3-46fe-b90f-b3731115030d
md"#### Part (j)

$\begin{align*}
\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{(\log{n})^{\log{n}}}{\frac{n}{\log{n}}} &= \lim_{n \to \infty} 10^{\log{\left(\frac{(\log{n})^{\log{n}}}{\frac{n}{\log{n}}}\right)}} \\
&= \lim_{n \to \infty} 10^{\left(\left(\log n\right)\log\left(\log\left(n\right)\right)-\log\left(n\right)+\log\left(\log n\right)\right)} \\
&= \lim_{n \to \infty} \frac{\left(10^{\log{n}}\right)^{\log{(\log{n})}} 10^{\log{(\log{n})}}}{10^{\log{n}}} \\
&= \lim_{n \to \infty} \frac{n^{\log\left(\log n\right)}\log n}{n} \\
&= \lim_{n \to \infty} n^{\log{(\log{n})} - 1} \log{n} \\
&= \infty \implies f = \Omega(g)
\end{align*}$"

# ╔═╡ f01338f8-23f0-4c81-9dbc-c1afc3c006f9
md"#### Part (k)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{\sqrt{n}}{(\log{n})^3} = \infty \implies f = \Omega(g)$"

# ╔═╡ 98bcd301-6c2c-4d0b-8ec3-ca891de91a05
md"#### Part (l)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n^{1/2}}{5^{\log_2{n}}} = 0 \implies f = O(g)$"

# ╔═╡ 66de4f55-61ef-4ef1-bba0-c8733ad6f008
md"#### Part (m)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n2^n}{3^n} = \lim_{n \to \infty} n \left(\frac{2}{3}\right)^n = 0 \implies f = O(g)$"

# ╔═╡ dd202daf-e918-4306-a38e-f263ea10a225
md"#### Part (n)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{2^n}{2^{n+1}} = \lim_{n \to \infty} \frac{2^n}{2^n \cdot 2} = \frac{1}{2} \implies f = \Theta(g)$"

# ╔═╡ 29a65c5a-479b-4e5b-8471-6d852b11d43b
md"#### Part (o)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{n!}{2^n} = \infty \implies f = \Omega(g)$"

# ╔═╡ 96745057-1360-4bc4-a4af-7fa6734aac26
md"#### Part (p)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{(\log{n})^{\log{n}}}{2^{(\log_2{n})^2}} = 0 \implies f = O(g)$"

# ╔═╡ 17d6abd2-3f11-4534-b031-8b934c59c583
md"#### Part (q)

$\lim_{n \to \infty} \frac{f(n)}{g(n)} = \lim_{n \to \infty} \frac{\sum_{i = 1}^n i^k}{n^{k + 1}} = \frac{n^k}{n^k \cdot n} = 0 \implies f = O(g)$"

# ╔═╡ 5db0b92d-e1dc-48ad-b434-14f359b58807
md"### Problem 0.2

Show that, if $c$ is a positive real number, then $g(n) = 1 + c + c^2 + \cdots + c^n$ is:

(a) $\Theta(1)$ if $c < 1$.

(b) $\Theta(n)$ if $c = 1$.

(c) $\Theta(c^n)$ if $c > 1$.

The moral: in big-Θ terms, the sum of a geometric series is simply the first term if the series is strictly decreasing, the last term if the series is strictly increasing, or the number of terms if the series is unchanging."

# ╔═╡ e434304a-88f6-4ad8-b41a-d44fe38094d3
md"#### Part (a)

$g(n) = 1 + c + c^2 + \cdots + c^n = \frac{1}{1 - c}$

$\lim_{n \to \infty} \frac{\frac{1}{1 - c}}{1} = \frac{1}{1 - c} \implies g = \Theta(1)$"

# ╔═╡ b78ef83e-f960-463f-8981-ef57a10906a0
md"#### Part (b)

$g(n) = 1 + c + c^2 + \cdots + c^n = 1 + 1 + 1 + \cdots + 1 = 1 + n$

$\lim_{n \to \infty} \frac{1+n}{n} = 1 \implies g = \Theta(n)$
"

# ╔═╡ 07038e25-aa70-4426-95e8-150868d10c30
md"#### Part (c)

$g(n) = 1 + c + c^2 + \cdots + c^n = \frac{1}{c - 1}$

Each term $c^i$ is dominated by $c^n$ for $i < n$, so $g(n) = \Theta(c^n)$.

**Source:** [https://web.cs.dal.ca/~nauzerk/csci3110/Assignments/a1\_solution\_f12.pdf](https://web.cs.dal.ca/~nauzerk/csci3110/Assignments/a1_solution_f12.pdf)
"

# ╔═╡ 5e48f296-0fe9-40b0-a238-420907425765
md"## 1 Algorithms with numbers"

# ╔═╡ 6be04797-3732-4034-9a21-bc3cb79ebcde
md"### Problem 1.2

Show that any binary integer is at most four times as long as the corresponding decimal integer.
For very large numbers, what is the ratio of these two lengths, approximately?"

# ╔═╡ 3bce9942-38c0-431c-be82-93ddead7b0c7
md"**Answer:**
An the size of an integer $N$ in base $a$ is the same as its size in base $b$, times a constant factor $\log_a{b}$.

Then the size of an integer $N$ represented in binary form is at most $\log_{2}{10} \approx 3.32$ times as long as its decimal representation.
For very large numbers, the ratio of these two lengths is approximately $\log_{2}{10} \approx 3.32$"

# ╔═╡ ca3fd00f-e17a-47d7-920b-df6fde825c90
log(2, 10)

# ╔═╡ 69227b69-01d5-47ba-84cb-50871aa7b056
md"### Problem 1.7

How long does the recursive multiplication algorithm (page 25) take to multiply an $n$-bit number by an $m$-bit number? Justify your answer."

# ╔═╡ da1776df-0827-483f-b5cf-dc77b279ea0e
function multiply(x, y)
	if y == 0 return 0 end
	z = multiply(x, y >> 1)
	if iseven(y)
		return z << 1
	else
		return x + (z << 1)
	end
end;

# ╔═╡ b79687ef-ad98-4e1f-880a-86b72349f423
md"**Answer:** The recursive multiplication algorithm takes $O(mn)$. It takes $m$ recursive calls since for each call $m$ is shifted to the right. Then it takes $O(n)$ bit operations for each call."

# ╔═╡ de5c451d-5af3-4b83-bb1d-7cc1c3b98048
md"### Problem 1.8

Justify the correctness of the recursive division algorithm given in page 26, and show that it takes time $O(n^2)$ on $n$-bit inputs."

# ╔═╡ f34b3872-9825-4bad-bdca-7da97906a212
function divide(x, y)
	if x == 0 return (q, r) = (0, 0) end
	(q, r) = divide(x >> 1, y)
	q = q << 1
	r = r << 1
	if isodd(x)
		r += 1
	end
	if r >= y
		r -= y
		q += 1
	end
	return (q, r)
end;

# ╔═╡ ea58b610-804b-4c4b-af04-c5411b52c350
md"**Answer:** The recursive division algorithm takes $n$ recursive calls since x is shifted to the right for each call and each call takes $O(n)$ bit operations, thus the algorithm takes $O(n^2)$ time.

We can justify its correctness using induction.

- Base case: $n = 2$, $x = 0$

- Inductive step: $n > 2$, $x > 1$, $x = yq + r$

- Case 1: $x$ is even, $q = 2q$, $r = 2r$

- Case 2: $x$ is odd, $q = 2q$, $r = 2r + 1$"

# ╔═╡ 9dbd6d07-870c-423d-90e0-63214ed1beb9
md"### Problem 1.18

Compute $\gcd(210, 588)$ two different ways: by finding the factorization of each number, and by using Euclid's algorithm.
"

# ╔═╡ d686c4e1-e9ba-4a34-8eca-03fd531a2cc0
md"#### Factorization

$\begin{align*}
210 &= 2 \cdot 105 = 2 \cdot 5 \cdot 21 = 2 \cdot 5 \cdot 3 \cdot 7 \\
588 &= 2 \cdot 294 = 2 \cdot 2 \cdot 147 = 2^2 \cdot 7 \cdot 21 = 2^2 \cdot 3 \cdot 7^2
\end{align*}$

$\gcd(210, 588) = 2 \cdot 3 \cdot 7 = 42$
"

# ╔═╡ f64a5f44-0db6-45f2-91ea-f95035bc0149
md"#### Euclid's algorithm
$\begin{align*}
\gcd(210, 588) &= \text{Euclid}(588, 210) \\
&= \text{Euclid}(210, 168) \\
&= \text{Euclid}(168, 42) \\
&= \text{Euclid}(42, 0) \\
&= 42
\end{align*}$"

# ╔═╡ 5318005b-e739-4617-abbf-88e4b69a4428
let
	function Euclid(a, b)
		if b == 0 return a end
		return Euclid(b, a % b)
	end

	Euclid(210, 588)
end

# ╔═╡ bf530c92-b3ea-4c70-9f08-2c8bc9aca846
md"### Problem 1.27

Consider an RSA key set with $p = 17$, $q = 23$, $N = 391$, and $e = 3$ (as in Figure 1.9).
What value of $d$ should be used for the secret key?
What is the encryption of the message $M = 41$?"

# ╔═╡ 462047c0-8891-438b-beda-900f7f25dbc5
md"**Answer:** First, verify that $e$ and $(p - 1)(q - 1)$ are relatively prime, i.e., $\gcd(e, (p - 1)(q - 1)) = 1$."

# ╔═╡ 442dab3f-9e22-49c1-a5d3-1c2f7fd592fd
let
	(p, q, e) = 17, 23, 3
	gcd(e, (p - 1) * (q - 1))
end

# ╔═╡ ddd61960-90a8-4d90-b151-eba80ce79569
md"The value of $d$ is the inverse of $e$ modulo $(p - 1)(q - 1)$"

# ╔═╡ cb0740b1-d802-47b5-92b0-c12a79c53026
let
	(p, q, e) = 17, 23, 3
	invmod(e, (p - 1) * (q - 1))
end

# ╔═╡ 653860b7-79aa-4802-994c-3eb0cdc3b672
md"The encryption of the message $M$ is

$M^e \bmod N$"

# ╔═╡ 5434c444-a214-410a-8b01-e31c1a5765a6
let
	(e, N, M) = 3, 391, 41
	powermod(M, e, N)
end

# ╔═╡ 0c775540-e00a-460a-9d7d-04f2b783981a
md"So $d = 235$ and the encryption of the message $M$ is 105."

# ╔═╡ 3cf4dc67-0d27-4198-97bd-0d50d73b7b19
md"## 2 Divide and conquer"

# ╔═╡ af2a2fa0-5ecd-4fab-8651-9ea8c6fa5fb4
md"### Problem 2.3

Section 2.2 describes a method for solving recurrence relationship which is based on analyzing the recursion tree and deriving a formula for the work done at each level.
Another (closely related) method is to expand out the recurrence a few times, until a pattern emerges.
For instance, let's start with the familiar $T(n) = 2T(n/2) + O(n)$.
Think of $O(n)$ as being $\leq cn$ for some constant $c$, so $T(n) \leq 2T(n/2) + cn$.
By repeatedly applying this rule, we can bound $T(n)$ in terms of $T(n/2)$, then $T(n/4)$, then $T(n/8)$, and so on, at each step getting closer to the value of $T(\cdot)$ we do know, namely $T(1) = O(1)$.

$\begin{align*}
T(n) &\leq 2T(n/2) + cn \\
&\leq 2[2T(n/4) + cn/2] + cn = 4T(n/4) + 2cn \\
&\leq 4[2T(n/8) + cn/4] + 2cn = 8T(n/8) + 3cn \\
&\leq 8[2T(n/16) + cn/8] + 3cn = 16T(n/16) + 4cn \\
&\vdots
\end{align*}$

A pattern is emerging... the general term is

$T(n) \leq 2^k T(n/2^k) + kcn.$

Plugging in $k = \log_2{n}$, we get $T(n) \leq nT(1) + cn \log_2{n} = O(n \log{n})$.

(a) Do the same thing for the recurrence $T(n) = 3T(n/2) + O(n)$.
    What is the general $k$th term in this case?
    And what value of $k$ should be plugged in to get the answer?

(b) Now try the recurrence $T(n) = T(n - 1) + O(1)$, a case which is not covered by the master theorem.
    Can you solve this too?
"

# ╔═╡ b750fc45-9547-467f-b37c-fcae7a78aacf
md"#### Part (a)

$\begin{align*}
T(n) &\leq 3T(n/2) + cn \\
&\leq 3[3T(n/4) + cn/2] + cn = 9T(n/4) + 3cn/2 \\
&\leq 9[3T(n/8) + cn/4] + 3cn/2 = 27T(n/8) + 7cn/4 \\
&\leq 27[3T(n/16) + cn/8] + 7cn/4 = 81T(n/16) + 15cn/8
\end{align*}$"

# ╔═╡ b279c66e-92d6-4453-82f1-cfb9acd306c1
md"### Problem 2.5

Solve the following recurrence relations and give a $\Theta$ bound for each of them.

$\begin{array}{ll}
\mathbf{(a)} & T(n) = 2T(n/3) + 1 \\
\mathbf{(b)} & T(n) = 5T(n/4) + n \\
\mathbf{(c)} & T(n) = 7T(n/7) + n \\
\mathbf{(d)} & T(n) = 9T(n/3) + n^2 \\
\mathbf{(e)} & T(n) = 8T(n/2) + n^3 \\
\mathbf{(f)} & T(n) = 49T(n/25) + n^{3/2} \log{n} \\
\mathbf{(g)} & T(n) = T(n-1) + 2 \\
\mathbf{(h)} & T(n) = T(n-1) + n^c, \text{ where } c \geq 1 \text{ is a constant } \\
\mathbf{(i)} & T(n) = T(n-1) + c^n, \text{ where } c > 1 \text{ is some constant } \\
\mathbf{(j)} & T(n) = 2T(n-1) + 1 \\
\mathbf{(k)} & T(n) = T(\sqrt{n}) + 1
\end{array}$"

# ╔═╡ 4db102db-2b38-478a-89b1-80afce55c78b
md"Master theorem. Suppose that $T(n)$ is a function on the nonnegative integers that satisfies the recurrence

$T(n) = aT\left(\frac{n}{b}\right) + f(n)$ where $n/b$ is integer division. Let $k = \log_b{a}$. Then,

- Case 1: If $f(n) = O(n^{k - \epsilon})$ for some constant $\epsilon > 0$, then $T(n) = \Theta(n^k)$.

- Case 2: If $f(n) = \Theta(n^k \log^p{n})$ for some constant $p \geq 0$, then $T(n) = \Theta(n^k \log^{p + 1} n$).

- Case 3: If $f(n) = \Omega(n^{k + \epsilon})$ for some constant $\epsilon > 0$, then $T(n) = \Theta(f(n))$."

# ╔═╡ f51e4978-db41-43c1-abc5-635f91b8d618
md"#### Part (a)

$\begin{gather*}
T(n) = 2T(n/3) + 1 \\
k = \log_b{a} = \log_3{2} \\
\text{Case 1} \quad f(n) = O(n^{\log_3{2} - \epsilon}) \\
T(n) = \Theta(n^{\log_3{2}})
\end{gather*}$
"

# ╔═╡ 5147b4b4-6445-4623-a049-7dc69e5c5779
md"#### Part (b)

$\begin{gather*}
T(n) = 5T(n/4) + n \\
k = \log_b{a} = \log_4{5} \\
\text{Case 1} \quad f(n) = O(n^{\log_4{5} - \epsilon}) \\
T(n) = \Theta(n^{\log_4{5}})
\end{gather*}$"

# ╔═╡ 79a6bcd6-fc17-49c0-b8bb-d4c262d2006a
md"#### Part (c)

$\begin{gather*}
T(n) = 7T(n/7) + n \\
k = \log_b{a} = \log_7{7} = 1 \\
\text{Case 1} \quad f(n) \neq O(n^{1 - \epsilon}) \\
\text{Case 2} \quad f(n) = \Theta(n \log^p{n}) \\
T(n) = \Theta(n \log{n})
\end{gather*}$"

# ╔═╡ 51e514c2-c633-4c11-95a4-f918fcdb0993
md"#### Part (d)

$\begin{gather*}
T(n) = 9T(n/3) + n^2 \\
k = \log_b{a} = \log_3{9} = 2 \\
\text{Case 1} \quad f(n) \neq O(n^{2 - \epsilon}) \\
\text{Case 2} \quad f(n) = \Theta(n^{2} \log^p{n}) \\
T(n) = \Theta(n^2 \log{n})
\end{gather*}$"

# ╔═╡ d55157cb-8476-452c-9766-2e50f1694b7b
md"#### Part (e)

$\begin{gather*}
T(n) = 8T(n/2) + n^3 \\
k = \log_b{a} = \log_2{8} = 3 \\
\text{Case 1} \quad f(n) \neq O(n^{3 - \epsilon}) \\
\text{Case 2} \quad f(n) = \Theta(n^{3 \log^p{n}}) \\
T(n) = \Theta(n^3 \log{n})
\end{gather*}$
"

# ╔═╡ 5e80011f-c44e-4a62-b905-3cc81b1fc1c7
md"#### Part (f)

$\begin{gather*}
T(n) = 49T(n/25) + n^{3/2} \log{n} \\
k = \log_{25}{49} \\
\text{Case 1} \quad f(n) \neq O(n^{\log_{25}{49} - \epsilon}) \\
\text{Case 2} \quad f(n) \neq \Theta(n^{\log_{25}{49}} \log^p{n}) \\
\text{Case 3} \quad f(n) = \Omega(n^{\log_{25}{49} + \epsilon}) \\
T(n) = \Theta(n^{3/2} \log{n})
\end{gather*}$"

# ╔═╡ 1d7d37a2-62f5-4433-a8cc-63e3a392bf99
md"#### Part (g)

$\begin{gather*}
T(n) = T(n - 1) + 2 \\
T(n - 1) = \Theta(n) \text{ dominates } 2 = \Theta(1) \\
T(n) = \Theta(n)
\end{gather*}$"

# ╔═╡ 47f5d93b-afb5-484f-91ff-736072cf8166
md"#### Part (h)

$\begin{gather*}
T(n) = T(n - 1) + n^c, \text{ where } c \geq 1 \text{ is a constant} \\
n^c = \Theta(n^c) \text{ dominates } T(n - 1) = \Theta(n) \\
T(n) = \Theta(n^{c+1})
\end{gather*}$"

# ╔═╡ ac43746a-c516-48b5-915c-05df5443c6ac
md"#### Part (i)

$\begin{gather*}
T(n) = T(n - 1) + c^n, \text{ where } c > 1 \text{ is some constant} \\
c^n = \Theta(c^n) \text{ dominates } T(n - 1) = \Theta(n) \\
T(n) = \Theta(c^n)
\end{gather*}$"

# ╔═╡ 88e903d1-e6e2-4527-b2b1-596c52de9428
md"#### Part (j)

$\begin{gather*}
T(n) = 2T(n - 1) + 1 \\
2T(n - 1) = \Theta(2^n) \text{ dominates } 1 = \Theta(1) \\
T(n) = \Theta(2^n)
\end{gather*}$"

# ╔═╡ e4c7a6f5-5fc9-4647-a56f-6aee6bd82783
md"#### Part (k)

$\begin{gather*}
T(n) = T(\sqrt{n}) + 1 \\
T(n) = \Theta(\log{(\log{n})})
\end{gather*}$"

# ╔═╡ dc597861-e2fc-4744-8c00-624cdad2cc2e
md"
$n^{(\log{n})^{10}}= 10^{\log{\left(n^{(\log{n})^{10}}\right)}} = 10^{(\log{n})^{11}} = n^{11} \implies T(n) = \Theta(n^{(\log{n})^{10}})$
"

# ╔═╡ 1b632b7e-97ea-4f81-9c37-b2915c026295
md"### Problem 2.12

How many lines, as a function of $n$ (in $\Theta(\cdot)$ form), does the following program print?
Write a recurrence and solve it.
You may assume $n$ is a power of 2."

# ╔═╡ c9faa850-9381-4f28-aff5-87e47ce03f16
function f(n)
	if n > 1
		println("still going")
		f(n/2)
		f(n/2)
	end
end;

# ╔═╡ 9da50846-fb1a-4b0a-b215-f0eba1274249
md"**Answer:**

$\begin{gather*}
T(n) = 2T(n/2) + 1\\
k = \log_2{2} = 1\\
T(n) = \Theta(n)
\end{gather*}$
"

# ╔═╡ Cell order:
# ╟─1ba17e2e-a9d2-11eb-3e83-bd0873b44b78
# ╟─ed49dc2f-7845-456d-9c41-481996339d98
# ╟─fe320c90-5d5d-410f-8715-c1ef4c8ab9fb
# ╟─92ce6a5e-0ad3-4905-9d40-f8048acc10af
# ╟─ca1cc38c-da0e-4d64-8fe2-c195f46f88cc
# ╟─d8a8ee9c-b782-4413-b665-b344332216f2
# ╟─a6bbe6a3-145b-498f-ac53-bb9bdf2cd72b
# ╟─4b240b0f-aeb2-43d4-8957-acc77bc1b044
# ╠═22507365-4ef7-4400-a9ba-7e93dc4a5203
# ╠═c51f5256-a278-4ce3-b0cf-7274d174f454
# ╠═2174fcff-76a9-47ab-99fd-662679baf33f
# ╠═7b96f661-a77b-4a40-8269-d14ec887a46a
# ╠═9db70690-4f7b-49a1-b259-a490c0f6787a
# ╠═dc49581c-2e2d-48c2-85f4-3aeeddbae8f0
# ╠═103dcea9-522b-45e0-bbc4-851450a755aa
# ╠═dfdf2ab2-797d-4744-a2ee-379320e120a5
# ╠═0c60b870-e122-4592-b42f-d5b64941d908
# ╠═37c60ab8-f2ef-41c3-8fdd-ecbf6189d625
# ╠═1754fc92-1036-4ce6-b209-f4f5bc75efba
# ╠═4bc44e0c-b2e1-4866-86df-da59331a620e
# ╠═c311334c-eee4-4171-93e3-56aa6e822c9f
# ╠═5481ab40-7e2c-4f35-94a1-90df39e2558f
# ╟─da5b0b05-dc8c-4bca-8b8c-c50b89f0018c
# ╟─d097fbac-2ab5-4e44-b9cc-98a694c2a7fa
# ╟─3918986d-8fcf-41ff-91bc-52e5eb6346c0
# ╟─f471306e-733c-421f-a243-1a0079989dfb
# ╟─33a9b77b-2c02-40de-8072-1164b6614e25
# ╟─be37a3a8-0a69-4b9c-93eb-d4025b2fccea
# ╟─4ee7fe91-2397-45af-865a-c252920bcf59
# ╟─02641c72-1bd8-4855-b169-1e67b7f840c8
# ╟─e09e3174-4b62-4d48-8640-202d1da7fb3f
# ╟─510a5fb5-e96b-44f7-9257-40fbfd4ea191
# ╟─5e4d7f29-974b-4d25-aa7b-8bf827ab1589
# ╟─5272ccd3-959d-4d29-8626-3d83cdb830a7
# ╟─93046ae6-acf3-46fe-b90f-b3731115030d
# ╟─f01338f8-23f0-4c81-9dbc-c1afc3c006f9
# ╟─98bcd301-6c2c-4d0b-8ec3-ca891de91a05
# ╟─66de4f55-61ef-4ef1-bba0-c8733ad6f008
# ╟─dd202daf-e918-4306-a38e-f263ea10a225
# ╟─29a65c5a-479b-4e5b-8471-6d852b11d43b
# ╟─96745057-1360-4bc4-a4af-7fa6734aac26
# ╟─17d6abd2-3f11-4534-b031-8b934c59c583
# ╟─5db0b92d-e1dc-48ad-b434-14f359b58807
# ╟─e434304a-88f6-4ad8-b41a-d44fe38094d3
# ╟─b78ef83e-f960-463f-8981-ef57a10906a0
# ╟─07038e25-aa70-4426-95e8-150868d10c30
# ╟─5e48f296-0fe9-40b0-a238-420907425765
# ╟─6be04797-3732-4034-9a21-bc3cb79ebcde
# ╟─3bce9942-38c0-431c-be82-93ddead7b0c7
# ╠═ca3fd00f-e17a-47d7-920b-df6fde825c90
# ╟─69227b69-01d5-47ba-84cb-50871aa7b056
# ╠═da1776df-0827-483f-b5cf-dc77b279ea0e
# ╟─b79687ef-ad98-4e1f-880a-86b72349f423
# ╟─de5c451d-5af3-4b83-bb1d-7cc1c3b98048
# ╠═f34b3872-9825-4bad-bdca-7da97906a212
# ╟─ea58b610-804b-4c4b-af04-c5411b52c350
# ╟─9dbd6d07-870c-423d-90e0-63214ed1beb9
# ╟─d686c4e1-e9ba-4a34-8eca-03fd531a2cc0
# ╟─f64a5f44-0db6-45f2-91ea-f95035bc0149
# ╠═5318005b-e739-4617-abbf-88e4b69a4428
# ╟─bf530c92-b3ea-4c70-9f08-2c8bc9aca846
# ╟─462047c0-8891-438b-beda-900f7f25dbc5
# ╠═442dab3f-9e22-49c1-a5d3-1c2f7fd592fd
# ╟─ddd61960-90a8-4d90-b151-eba80ce79569
# ╠═cb0740b1-d802-47b5-92b0-c12a79c53026
# ╟─653860b7-79aa-4802-994c-3eb0cdc3b672
# ╠═5434c444-a214-410a-8b01-e31c1a5765a6
# ╟─0c775540-e00a-460a-9d7d-04f2b783981a
# ╟─3cf4dc67-0d27-4198-97bd-0d50d73b7b19
# ╟─af2a2fa0-5ecd-4fab-8651-9ea8c6fa5fb4
# ╟─b750fc45-9547-467f-b37c-fcae7a78aacf
# ╟─b279c66e-92d6-4453-82f1-cfb9acd306c1
# ╟─4db102db-2b38-478a-89b1-80afce55c78b
# ╟─f51e4978-db41-43c1-abc5-635f91b8d618
# ╟─5147b4b4-6445-4623-a049-7dc69e5c5779
# ╟─79a6bcd6-fc17-49c0-b8bb-d4c262d2006a
# ╟─51e514c2-c633-4c11-95a4-f918fcdb0993
# ╟─d55157cb-8476-452c-9766-2e50f1694b7b
# ╟─5e80011f-c44e-4a62-b905-3cc81b1fc1c7
# ╟─1d7d37a2-62f5-4433-a8cc-63e3a392bf99
# ╟─47f5d93b-afb5-484f-91ff-736072cf8166
# ╟─ac43746a-c516-48b5-915c-05df5443c6ac
# ╟─88e903d1-e6e2-4527-b2b1-596c52de9428
# ╟─e4c7a6f5-5fc9-4647-a56f-6aee6bd82783
# ╟─dc597861-e2fc-4744-8c00-624cdad2cc2e
# ╟─1b632b7e-97ea-4f81-9c37-b2915c026295
# ╠═c9faa850-9381-4f28-aff5-87e47ce03f16
# ╟─9da50846-fb1a-4b0a-b215-f0eba1274249

max_stars_repo_name: airicbear/finals-review

max_issues_repo_name: airicbear/finals-review

max_forks_repo_name: airicbear/finals-review


text:
"""
    readlabels(io::IO, index::Integer)

Jumps to the position of `io` where the byte for the `index`'th
label is located and returns the byte at that position as `UInt8`
"""
function readlabels(io::IO, index::Integer)
    seek(io, LABELOFFSET + (index - 1))
    read(io, UInt8)::UInt8
end

"""
    readlabels(io::IO, indices::AbstractVector)

Reads the byte for each label-index in `indices` and stores them
in a `Vector{UInt8}` of length `length(indices)` in the same order
as denoted by `indices`.
"""
function readlabels(io::IO, indices::AbstractVector)
    labels = Array{UInt8}(undef, length(indices))
    dst_index = 1
    for src_index in indices
        labels[dst_index] = readlabels(io, src_index)
        dst_index += 1
    end
    labels::Vector{UInt8}
end

"""
    readlabels(file::AbstractString, [indices])

Reads the label denoted by `indices` from `file`. The given `file`
is assumed to be in the MNIST label-file format, as it is described
on the official homepage at http://yann.lecun.com/exdb/mnist/

- if `indices` is an `Integer`, the single label is returned as `UInt8`.

- if `indices` is a `AbstractVector`, the labels are returned as
  a `Vector{UInt8}`, length `length(indices)` in the same order as
  denoted by `indices`.

- if `indices` is ommited all all are returned
  (as `Vector{UInt8}` as described above)
"""
function readlabels(file::AbstractString, index::Integer)
    gzopen(file, "r") do io
        _, nlabels = readlabelheader(io)
        @assert minimum(index) >= 1 && maximum(index) <= nlabels
        readlabels(io, index)
    end::UInt8
end

function readlabels(file::AbstractString, indices::AbstractVector)
    gzopen(file, "r") do io
        _, nlabels = readlabelheader(io)
        @assert minimum(indices) >= 1 && maximum(indices) <= nlabels
        readlabels(io, indices)
    end::Vector{UInt8}
end

function readlabels(file::AbstractString)
    gzopen(file, "r") do io
        _, nlabels = readlabelheader(io)
        readlabels(io, 1:nlabels)
    end::Vector{UInt8}
end


text:
module SplittingOrderConditions

import Base: start, next, done

export Lyndon, MultiFor, multinomial_coeff, generate_equations, call_form
export generate_equations_ABC, get_ABC_coeffs
export compute_residua, compute_residua_ABC, get_lem, get_p_lem

using Combinatorics

immutable Lyndon
    s::Int
    n::Int
end

Base.start(::Lyndon) = Int[-1]

function Base.next(L::Lyndon, w::Vector{Int})
    if w == [-1]
        w[end] += 1
        return (copy(w), w)
    end
    m = length(w)
    while length(w) < L.n
        push!(w, w[end-m+1])
    end
    while length(w) > 0 && w[end] == L.s - 1
        pop!(w)
    end
    w[end] += 1
    (copy(w), w)
end

Base.done(L::Lyndon, w::Vector{Int}) = w == [L.s-1]

immutable MultiFor
    k::Array{Int,1}
end

Base.start(MF::MultiFor) = Int[]

Base.done(MF::MultiFor, k::Array{Int,1}) = MF.k==k

function Base.next(MF::MultiFor, k::Array{Int,1})
    if k==Int[]
        k = zeros(Int, length(MF.k))
        return(copy(k), k)
    end
    for i=1:length(k)
        if k[i]<MF.k[i]
            k[i] += 1
            for j = 1:i-1
                k[j] = 0
            end
            return (copy(k), k)
        end
    end
end

multinomial_coeff(q::Int, k::Array{Int,1}) = div(factorial(q), prod([factorial(i) for i in k]))

function generate_equations(q::Int, s::Int; symmetric::Bool=false, palindromic::Bool=false,
                            last_b_is_zero::Bool=false, with_brackets::Bool=true)
    eqs = Dict{Array{Int,1}, AbstractString}()
    if symmetric && iseven(q)
        return eqs
    end
    if symmetric
        palindromic = false
        last_b_is_zero = true
        s2 = ceil(Int,s/2)
        s3 = isodd(s)?s2-1:s2
        if with_brackets
            a = vcat([string("a[", j, "]") for j=1:s2], [string("a[", j, "]") for j=s3:-1:1])
        else
            a = vcat([string("a", j) for j=1:s2], [string("a", j) for j=s3:-1:1])
        end
        s2 = ceil(Int,(s-1)/2)
        s3 = iseven(s)?s2-1:s2
        if with_brackets
            b = vcat([string("b[", j, "]") for j=1:s2], [string("b[", j, "]") for j=s3:-1:1],["0"])
        else
            b = vcat([string("b", j) for j=1:s2], [string("b", j) for j=s3:-1:1],["0"])
        end
    elseif palindromic
        last_b_is_zero = false
        s2 = ceil(Int,s/2)
        s3 = isodd(s)?s2-1:s2
        if with_brackets
            a = vcat([string("a[", j, "]") for j=1:s2], [string("b[", j, "]") for j=s3:-1:1])
        else
            a = vcat([string("a", j) for j=1:s2], [string("b", j) for j=s3:-1:1])
        end
        b = reverse(a)
    else
        if with_brackets
            a = [string("a[", j, "]") for j=1:s]
            b = [string("b[", j, "]") for j=1:s]
        else
            a = [string("a", j) for j=1:s]
            b = [string("b", j) for j=1:s]
        end
    end
    if palindromic
        for j in Lyndon(2, q)
            if length(j)==q && !haskey(eqs, [x==0?1:0 for x in reverse(j)])
                eqs[j] = ""
            end
        end
    else
        for j in Lyndon(2, q)
            if length(j)==q
                eqs[j] = ""
            end
        end
    end

    for k0 in Combinatorics.WithReplacementCombinations(1:s,q)
        k = zeros(Int,s)
        for j in k0
            k[j]+=1
        end
        for l = MultiFor(k)
            W=Int[]
            for j=s:-1:1
                append!(W, ones(Int, l[j]))
                append!(W, zeros(Int, k[j]-l[j]))
            end
            if W in keys(eqs) && !(last_b_is_zero && l[s]>0)
                coeff =  multinomial_coeff(q, k)
                for j=1:s
                    coeff = coeff * binomial(k[j], l[j])
                end
                h = string(coeff)
                for j=1:s
                    if k[j]-l[j]>0
                        h = string(h,"*",a[j])
                    end
                    if k[j]-l[j]>1
                        h = string(h,"^",k[j]-l[j])
                    end
                    if l[j]>0
                        h = string(h,"*",b[j])
                    end
                    if l[j]>1
                        h = string(h,"^",l[j])
                    end
                end
                eqs[W] = string(eqs[W],"+",h)
            end
        end
    end
    for (W, h) in eqs
        eqs[W] = string(h,"-1")
    end
    eqs
end

function generate_equations_ABC(q::Int, s::Int; symmetric::Bool=false, palindromic::Bool=false,
                            last_c_is_zero::Bool=false, last_b_is_zero::Bool=false, with_brackets::Bool=true)
    eqs = Dict{Array{Int,1}, AbstractString}()
    if symmetric && iseven(q)
        return eqs
    end
    if last_b_is_zero
       last_c_is_zero = true
    end
    if symmetric
        error("symmertric not implemented.")
        #palindromic = false
        #last_c_is_zero = true
        #last_b_is_zero = true
        #s2 = ceil(Int,s/2)
        #s3 = isodd(s)?s2-1:s2
        #if with_brackets
        #    a = vcat([string("a[", j, "]") for j=1:s2], [string("a[", j, "]") for j=s3:-1:1])
        #else
        #    a = vcat([string("a", j) for j=1:s2], [string("a", j) for j=s3:-1:1])
        #end
        #s2 = ceil(Int,(s-1)/2)
        #s3 = iseven(s)?s2-1:s2
        #if with_brackets
        #    b = vcat([string("b[", j, "]") for j=1:s2], [string("b[", j, "]") for j=s3:-1:1],["0"])
        #else
        #    b = vcat([string("b", j) for j=1:s2], [string("b", j) for j=s3:-1:1],["0"])
        #end
    elseif palindromic
        last_c_is_zero = false
        last_b_is_zero = false
        s2 = ceil(Int,s/2)
        s3 = isodd(s)?s2-1:s2
        if with_brackets
            a = vcat([string("a[", j, "]") for j=1:s2], [string("c[", j, "]") for j=s3:-1:1])
            b = vcat([string("b[", j, "]") for j=1:s2], [string("b[", j, "]") for j=s3:-1:1])
        else
            a = vcat([string("a", j) for j=1:s2], [string("c", j) for j=s3:-1:1])
            b = vcat([string("b", j) for j=1:s2], [string("b", j) for j=s3:-1:1])
        end
        c = reverse(a)
    else
        if with_brackets
            a = [string("a[", j, "]") for j=1:s]
            b = [string("b[", j, "]") for j=1:s]
            c = [string("c[", j, "]") for j=1:s]
        else
            a = [string("a", j) for j=1:s]
            b = [string("b", j) for j=1:s]
            c = [string("c", j) for j=1:s]
        end
    end
    if palindromic
        rr = (2, 1, 0)
        for j in Lyndon(3, q)
            if length(j)==q && !haskey(eqs, [rr[x+1] for x in reverse(j)])
                eqs[j] = ""
            end
        end
    else
        for j in Lyndon(3, q)
            if length(j)==q
                eqs[j] = ""
            end
        end
    end

    for k0 in Combinatorics.WithReplacementCombinations(1:s,q)
        k = zeros(Int,s)
        for j in k0
            k[j]+=1
        end
        for lA = MultiFor(k)
        for lB = MultiFor(k-lA)
        for lC = MultiFor(k-lA-lB)
            W=Int[]
            for j=s:-1:1
                append!(W, 2*ones(Int, lC[j]))
                append!(W, ones(Int, lB[j]))
                append!(W, zeros(Int, lA[j]))
            end
            if W in keys(eqs) && !(last_b_is_zero && lB[s]>0)  && !(last_c_is_zero && lC[s]>0)
                coeff =  multinomial_coeff(q, k)
                for j=1:s
                    coeff = coeff * multinomial_coeff(k[j], [lA[j], lB[j], lC[j]])
                end
                h = string(coeff)
                for j=1:s
                    if lA[j]>0
                        h = string(h,"*",a[j])
                    end
                    if lA[j]>1
                        h = string(h,"^",lA[j])
                    end
                    if lB[j]>0
                        h = string(h,"*",b[j])
                    end
                    if lB[j]>1
                        h = string(h,"^",lB[j])
                    end
                    if lC[j]>0
                        h = string(h,"*",c[j])
                    end
                    if lC[j]>1
                        h = string(h,"^",lC[j])
                    end
                end
                eqs[W] = string(eqs[W],"+",h)
            end
        end
        end
        end
    end
    for (W, h) in eqs
        eqs[W] = string(h,"-1")
    end
    eqs
end

function compute_residua(q::Int, a::Vector, b::Vector)
    s = length(a)
    T = typeof(a[1])
    eqs = Dict{Array{Int,1}, T}()

    aa = zeros(T, s, q+1)
    bb = zeros(T, s, q+1)
    aa[:,1] = ones(T, s)
    bb[:,1] = ones(T, s)
    for k=1:q
        aa[:,k+1] = aa[:,k].*a
        bb[:,k+1] = bb[:,k].*b
    end

    for j in Lyndon(2, q)
        if length(j)==q
            eqs[j] = -one(T)
        end
    end

    for k0 in Combinatorics.WithReplacementCombinations(1:s,q)
        k = zeros(Int,s)
        for j in k0
            k[j]+=1
        end
        for l = MultiFor(k)
            W=Int[]
            for j=s:-1:1
                append!(W, ones(Int, l[j]))
                append!(W, zeros(Int, k[j]-l[j]))
            end
            if W in keys(eqs)
                h =  multinomial_coeff(q, k)
                for j=1:s
                    h *= binomial(k[j], l[j])
                end
                for j=1:s
                    h *= aa[j, k[j]-l[j]+1] * bb[j, l[j]+1]
                end
                eqs[W] += h
            end
        end
    end
    eqs
end

get_lem(q, a, b) = norm(collect(values(compute_residua(q, a, b))))

function get_p_lem(a,b)
    q = 1
    while true
        lem = get_lem(q,a,b)
        if lem>1e-11
            return q-1, lem
        end
        q += 1
    end
end

function get_ABC_coeffs(a,b)
    s=length(a)
    T = typeof(a[1])
    S = 2*s-1
    A = zeros(T, S)
    B = zeros(T, S)
    C = zeros(T, S)
    h = a[1]
    for k=1:s
        K = 2*k-1
        A[K] = a[k]
        B[K] = h
        C[K] = b[k]
        h = b[k]-h
        if k<s
            A[K+1] = 0.0
            B[K+1] = h
            C[K+1] = 0.0
            h = a[k+1]-h
        end
    end
    A, B, C
end

function compute_residua_ABC(q::Int, a::Vector, b::Vector, c::Vector)
    s = length(a)
    T = typeof(a[1])
    eqs = Dict{Array{Int,1}, T}()

    aa = zeros(T, s, q+1)
    bb = zeros(T, s, q+1)
    cc = zeros(T, s, q+1)
    aa[:,1] = ones(T, s)
    bb[:,1] = ones(T, s)
    cc[:,1] = ones(T, s)
    for k=1:q
        aa[:,k+1] = aa[:,k].*a
        bb[:,k+1] = bb[:,k].*b
        cc[:,k+1] = cc[:,k].*c
    end

    for j in Lyndon(3, q)
        if length(j)==q
            eqs[j] = -one(T)
        end
    end

    for k0 in Combinatorics.WithReplacementCombinations(1:s,q)
        k = zeros(Int,s)
        for j in k0
            k[j]+=1
        end
        for lA = MultiFor(k)
        for lB = MultiFor(k-lA)
        for lC = MultiFor(k-lA-lB)
            W=Int[]
            for j=s:-1:1
                append!(W, 2*ones(Int, lC[j]))
                append!(W, ones(Int, lB[j]))
                append!(W, zeros(Int, lA[j]))
            end
            if W in keys(eqs)
                h =  multinomial_coeff(q, k)
                for j=1:s
                    h *= multinomial_coeff(k[j], [lA[j], lB[j], lC[j]])
                end
                for j=1:s
                    h *= aa[j, lA[j]+1] * bb[j, lB[j]+1] * cc[j, lC[j]+1]
                end
                eqs[W] += h
            end
        end
        end
        end
    end
    eqs
end

get_lem(q, a, b, c) = norm(collect(values(compute_residua_ABC(q, a, b, c))))

function get_p_lem(a,b,c)
    q = 1
    while true
        lem = get_lem(q,a,b,c)
        if lem>1e-11
            return q-1, lem
        end
        q += 1
    end
end

function call_form(input::ASCIIString; threads=1)
    attempts = 0
    FORM_PID=0
    path = joinpath(dirname(@__FILE__), "../deps/bin")
    out = "Error initializing preset external channels\n"
    while out == "Error initializing preset external channels\n"
        (so,si,pr) = readandwrite(`$path/tform -w$(threads) -t /tmp -M  -q -pipe 0,1 $path/pipe.frm`)
        FORM_PID = readuntil(so,'\n')[1:end-1];
        print(si, FORM_PID,',',getpid(),"\n\n", input,"\n\n")
        close(si)
        out = readall(so)
        close(si)
        close(so)
        close(pr)
        attempts +=1
    end
    println(STDERR, "# of attemps to call form = ",attempts)
    try
        rm("/tmp/xform$FORM_PID.str") # delete generated temporary file
    end
    out
end

function compile_fg(fun::AbstractString, n::Integer; threads=1)
    fun = replace(replace(fun,"[", "("), "]", ")")
    input = string("""
Off Statistics;
V x;
S u;
""","Local f0 = ",
    fun,
    ";\n",
"""
#write <> "begin"
#write <> "  if length(g)==0"
Print;
Format O3;
Format maple;
.sort;
#write <> "    return f0"
#write <> "  else"
I i;
S m, xx, u;
Hide f0;
#do i=1,$(n)
    Local f'i' = f0;
    id x('i') = xx;
    id xx^m? = m*xx^(m-1);
    id xx = x('i');
    .sort
    Hide f'i';
#enddo
.sort;
""","Local H = ",
    join(["u^$(j+1)*f$(j)" for j=0:n],:+),
    ";\n",
"""
B u;
.sort
#optimize H
B u;
.sort
""",
    join(["Local F$(j) = H[u^$(j+1)];\n" for j=0:n], ""),
"""
.sort
#write <> "%4O",
""", "#write <> \"\n   f=%e   ",
join(["g($(j))=%e" for j=1:n],"    "),
"\", ",
join(["F$(j)" for j=0:n],","),"\n",
"""
#write <> "    return f"
#write <> "  end"
#write <> "end"
.end
""")
    out = call_form(input, threads=threads)
    out = replace(replace(out,"(", "["), ")", "]")
    out = replace(out,"length[g]", "length(g)")
    eval(parse(string("(x,g)->",out)))
end

end


text:
using Combinatorics

# Part 2
checksum = 0;
for line in eachline(STDIN)
  input = parse.(Int, split(line, ' '))

  for (x, y) in combinations(input, 2)
    if x % y == 0
      checksum += Int(x / y)
      break
    end

    if y % x == 0
      checksum += Int(y / x)
      break
    end
  end
end

println("Part 2: $checksum")


text:
module SamplesTest

using PopGenCore
using PopGenSims
using Test

cats = @nancycats ;
d = Dict("1" => 5, "4" => 7, "11" => 3)

@testset "Simulate Samples" begin
    @testset "sample_locus" begin
        d = Dict(133 => 0.125, 135 => 0.5625, 143 => 0.25, 137 => 0.0625)
        @test PopGenSims.sample_locus(d, 3, 2) |> length == 3
        @test PopGenSims.sample_locus(d, 3, 2) |> first |> length == 2
        @test PopGenSims.sample_locus(d, 3, 3) |> length == 3
        @test PopGenSims.sample_locus(d, 3, 3) |> first |> length == 3
    end

    @testset "simulate errors" begin
        @test_throws ArgumentError simulate(cats)
        @test_throws ArgumentError simulate(cats, n = 1, scale = 1)
        @test_throws ArgumentError simulate(cats, n = d, scale = 2)
    end

    @testset "simulate flat" begin
        sims = simulate(cats , n = 100)
        @test sims isa PopData
        @test sims.metadata.samples == 1700
        @test length(sims.genodata.name) == 15300
    end

    @testset "simulate proportional" begin
        sims = simulate(cats , scale = 1)
        @test sims isa PopData
        @test sims.metadata.populations == cats.metadata.populations
        @test sims.metadata.samples == cats.metadata.samples
        @test length(sims.genodata.name) == length(cats.genodata.name)
        sims = simulate(cats , scale = 4)
        @test sims isa PopData
        @test sims.metadata.populations == cats.metadata.populations
        @test sims.metadata.samples == cats.metadata.samples * 4
        @test length(sims.genodata.name) == length(cats.genodata.name) * 4
    end

    @testset "simulate arbitrary" begin
        sims = simulate(cats , n = d)
        @test sims isa PopData
        @test sims.metadata.samples == sum(values(d))
        @test length(sims.genodata.name) == sum(values(d)) * cats.metadata.loci
        @test length(d) == sims.metadata.populations
    end
end

end # module


text:
"""
    FeatureTransforms.TestUtils

Provides fake [`Transform`](@ref)s and utilities for testing purposes only.

Each fake [`Transform`](@ref) has different a different `cardinality`: `OneToOne`, OneToMany`,
`ManyToOne`, or `ManyToMany`. So when users extend FeatureTransforms.jl for new data types
they only need to test against these 4 fakes to guarantee their type can support any
[`Transform`](@ref) in the package.

Similarly, [`is_transformable`](@ref) is used to check that the output of a `transform`
pipeline is a transformable type.
"""

module TestUtils

using ..FeatureTransforms
using ..FeatureTransforms: OneToOne, OneToMany, ManyToOne, ManyToMany
using InteractiveUtils: methodswith
using Tables

export FakeOneToOneTransform, FakeOneToManyTransform
export FakeManyToOneTransform, FakeManyToManyTransform
export is_transformable

for C in (:OneToOne, :OneToMany, :ManyToOne, :ManyToMany)
    FT = Symbol(:Fake, C, :Transform)
    @eval begin
        """
            $($FT) <: Transform

        A fake `$($C)` transform for test purposes. Calling `apply` will return an
        array of ones with a size and dimension matching the `cardinality` of the transform.
        """
        struct $FT <: Transform end
        FeatureTransforms.cardinality(::$FT) = $C()
    end
end

function FeatureTransforms._apply(A, ::FakeOneToOneTransform; kwargs...)
    return replace(one, A)
end

function FeatureTransforms._apply(A, ::FakeOneToManyTransform; kwargs...)
    return hcat(replace(one, A), replace(one, A))
end

function FeatureTransforms._apply(A, ::FakeManyToOneTransform; kwargs...)
    return replace(one, first(A))
end

function FeatureTransforms._apply(A, ::FakeManyToManyTransform; kwargs...)
    return hcat(replace(one, A), replace(one, A))
end

"""
    is_transformable(x)

Determine if `x` is both a valid input and output of any [`Transform`](@ref), i.e. that it
has an `apply` method defined and therefore follows the [`transform`](@ref) interface.
"""
function is_transformable(T::Type)
    Tables.istable(T) && return true  # cannot directly check against the method with no type
    return !isempty(methodswith(T, FeatureTransforms.apply; supertypes=true))
end
is_transformable(::T) where T = is_transformable(T)
# Need this to get around using methodswith, which would otherwise return true
is_transformable(::Type{<:Transform}) = false

end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function getpotstr(pot)

if pot==0
    sp0 = "(T *u, T *rho, T *mu, T *eta, int *kappa, int nrho, int nmu, int neta, int nkappa, int ng)\n";
    sp1 = "(double *, double *, double *, double *, int*, int, int, int, int, int);\n";
    sp2 = "(float *, float *, float *, float *, int *, int, int, int, int, int);\n";
    sp3 = "(u, rho, mu, eta, kappa, nrho, nmu, neta, nkappa, ng);\n";
elseif pot==1
    sp0 = "(T *u, T *xi, T *qi, int *ti, int *ai, T *mu, T *eta, int *kappa, int dim, int ncq, int nmu, int neta, int nkappa, int ng)\n";
    sp1 = "(double *, double *, double *, int *, int *, double *, double *, int*, int, int, int, int, int, int);\n";
    sp2 = "(float *, float *, float *, int *, int *, float *, float *, int *, int, int, int, int, int, int);\n";
    sp3 = "(u, xi, qi, ti, ai, mu, eta, kappa, dim, ncq, nmu, neta, nkappa, ng);\n";
elseif pot==2
    sp0 = "(T *u, T *xij, T *qi, T *qj, int *ti, int *tj, int *ai, int *aj, T *mu, T *eta, int *kappa, int dim, int ncq, int nmu, int neta, int nkappa, int ng)\n";
    sp1 = "(double *, double *, double *, double *, int *, int *, int *, int *, double *, double *, int*, int, int, int, int, int, int);\n";
    sp2 = "(float *, float *, float *, float *, int *, int *, int *, int *, float *, float *, int *, int, int, int, int, int, int);\n";
    sp3 = "(u, xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, dim, ncq, nmu, neta, nkappa, ng);\n";
elseif pot==3
    sp0 = "(T *u, T *xij, T *xik, T *qi, T *qj, T *qk, int *ti, int *tj, int *tk, int *ai, int *aj, int *ak, T *mu, T *eta, int *kappa, int dim, int ncq, int nmu, int neta, int nkappa, int ng)\n";
    sp1 = "(double *, double *, double *, double *, double *, double *, int *, int *, int *, int *, int *, int *, double *, double *, int*, int, int, int, int, int, int);\n";
    sp2 = "(float *, float *, float *, float *, float *, float *, int *, int *, int *, int *, int *, int *, float *, float *, int *, int, int, int, int, int, int);\n";
    sp3 = "(u, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, mu, eta, kappa, dim, ncq, nmu, neta, nkappa, ng);\n";
elseif pot==4
    sp0 = "(T *u, T *xij, T *xik, T *xil, T *qi, T *qj, T *qk, T *ql, int *ti, int *tj, int *tk, int *tl, int *ai, int *aj, int *ak, int *al, T *mu, T *eta, int *kappa, int dim, int ncq, int nmu, int neta, int nkappa, int ng)\n";
    sp1 = "(double *, double *, double *, double *, double *, double *, double *, double *, int *, int *, int *, int *, int *, int *, int *, int *, double *, double *, int*, int, int, int, int, int, int);\n";
    sp2 = "(float *, float *, float *, float *, float *, float *, float *, float *, int *, int *, int *, int *, int *, int *, int *, int *, float *, float *, int *, int, int, int, int, int, int);\n";
    sp3 = "(u, xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, dim, ncq, nmu, neta, nkappa, ng);\n";
end

return sp0, sp1, sp2, sp3

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.05"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
data1 = CSV.read(datadir("purchaseData.csv"))

# ╔═╡ 3d63b2d0-243e-11eb-1312-ab095ba05734
data2 = CSV.read(datadir("purchaseData.csv"), copycols=true)

# ╔═╡ 3d63db0c-243e-11eb-2aaf-6b546a467ef9
try data1[1, :Name] = "YARDEN" catch; @warn "Cannot: data1 is immutable" end

# ╔═╡ 3d6454ec-243e-11eb-27ef-0f63bc8e3e5d
begin
	data2[1, :Name] = "YARDEN"
	first(data2, 3)
end

# ╔═╡ 3d6aa5cc-243e-11eb-1a5b-0b128c7a8002
begin
	data1[!, :Price] ./= 1000
	rename!(data1, :Price=>Symbol("Price(000's)"))
	first(data1, 3)
end

# ╔═╡ 3d7655d4-243e-11eb-1e29-ebfc5d9a583a
begin
	replace!(data2[!,:Grade], ["E"=>"F", "D"=>"E"]...)
	first(data2, 3)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.05"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═3d63b2d0-243e-11eb-1312-ab095ba05734
# ╠═3d63db0c-243e-11eb-2aaf-6b546a467ef9
# ╠═3d6454ec-243e-11eb-27ef-0f63bc8e3e5d
# ╠═3d6aa5cc-243e-11eb-1a5b-0b128c7a8002
# ╠═3d7655d4-243e-11eb-1e29-ebfc5d9a583a
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
function solve_corner!(cost::CostParameters, mp::RLSModel, ic::Int64, τ::Array{Stage,1}, ESS)
    c = cost.C[ic]
    P1 = 0.
    vN1 = mp.β*((1-cost.π[ic])*Φ(τ[ic+1].EQs[ic, ic, ESS[ic,ic,ic+1]].vN1, τ[ic+1].EQs[ic, ic, ESS[ic,ic,ic+1]].vI1)+Φ(0, -K(c, cost)))/(1 - cost.π[ic]*mp.β)
    vI1 = vN1 - K(c, cost)
    vN2 = mp.β*( (1-cost.π[ic])*Φ(τ[ic+1].EQs[ic, ic, ESS[ic,ic,ic+1]].vN2, τ[ic+1].EQs[ic, ic, ESS[ic,ic,ic+1]].vI2)+Φ(0, -K(c, cost)))/(1 - cost.π[ic]*mp.β)
    vI2 = vN2 - K(c, cost)
    P2 = 0.
    @inbounds τ[ic].EQs[ic,ic,1] = EQ(P1, vN1, vI1, P2, vN2, vI2)
end

function solve_edges!(cost::CostParameters, mp::RLSModel, ic::Int64, τ::Array{Stage,1}, ess::ESS)
    c = cost.C[ic]
    pi = cost.π[ic]

    # if at (c1,c,c) edge
    @inbounds for ic1 = 1:ic-1
        c1 = cost.C[ic1]
        vI1 = -K(c,cost) + mp.β*H1(ic,ic,ic, pi, τ, ess)

        a = vI1
        A = mp.β*(1-pi)*Φ(τ[ic+1].EQs[ic1,ic,ess[ic1,ic,ic+1]].vN1, τ[ic+1].EQs[ic1, ic, ess[ic1,ic,ic+1]].vI1)
        B = mp.β*(1-pi)
        vN1 = 0
        if a == 0
            if  A >= 0
                vN1 = -A/(B-1)
            elseif A < 0
                vN1 = A
            end
        elseif a > 0
            if A >= a
                vN1 = -A/(B-1)
            elseif 0 < A < a
                if B <= (a-A)/a
                    vN1 = a*B+A
                elseif (a-A)/a < B
                    vN1 = -A/(B-1)
                end
            elseif A <= 0
                x = a*B+A
            end
        elseif a < 0
            if A >= 0
                vN1 = -A/(B-1)
            elseif a < A < 0
                if B <= (a-A)/a
                    vN1 = -A/(B-1)
                elseif (a-A)/a < B
                    vN1 = a*B+A
                end
            elseif A <= a
                vN1 = a*B+A
            end
        end

        P1 = convert(Float64,vI1>vN1)

        vN2 = (r2(c1,c) + mp.β*(P1*((1-pi)*Φ(τ[ic+1].EQs[ic, ic, ess[ic,ic,ic+1]].vN2,τ[ic+1].EQs[ic, ic, ess[ic,ic,ic+1]].vI2) + pi*Φ(τ[ic].EQs[ic, ic, ess[ic,ic,ic]].vN2,τ[ic].EQs[ic, ic, ess[ic,ic,ic]].vI2))+(1-P1)*((1-pi)*Φ(τ[ic+1].EQs[ic1, ic, ess[ic1,ic,ic+1]].vN2,τ[ic+1].EQs[ic1, ic, ess[ic1,ic,ic+1]].vI2)+pi*Φ(0,-K(c, cost)))))/(1-mp.β*(1-P1)*pi)
        vI2 = vN2 - K(c, cost)
        P2 = 0.

        @inbounds τ[ic].EQs[ic1,ic,1] = EQ(P1,
                                           vN1,
                                           vI1,
                                           P2,
                                           vN2,
                                           vI2)
    end

    # if at (c,c2,c) edge
    @inbounds for ic2 = 1:ic-1
        c2 = cost.C[ic2]
        vI2 = -K(c,cost) + mp.β*((1-pi)*Φ(τ[ic+1].EQs[ic, ic, ess[ic,ic,ic+1]].vN2,τ[ic+1].EQs[ic, ic, ess[ic,ic,ic+1]].vI2) + pi*Φ(τ[ic].EQs[ic, ic, ess[ic,ic,ic]].vN2,τ[ic].EQs[ic, ic, ess[ic,ic,ic]].vI2))
        a = vI2
        A = mp.β*(1-pi)*Φ(τ[ic+1].EQs[ic,ic2,ess[ic,ic2,ic+1]].vN2, τ[ic+1].EQs[ic, ic2, ess[ic,ic2,ic+1]].vI2)
        B = mp.β*(1-pi)
        vN2 = 0
        if a == 0
            if  A >= 0
                vN2 = -A/(B-1)
            elseif A < 0
                vN2 = A
            end
        elseif a > 0
            if A >= a
                vN2 = -A/(B-1)
            elseif 0 < A < a
                if B <= (a-A)/a
                    vN2 = a*B+A
                elseif (a-A)/a < B
                    vN2 = -A/(B-1)
                end
            elseif A <= 0
                x = a*B+A
            end
        elseif a < 0
            if A >= 0
                vN2 = -A/(B-1)
            elseif a < A < 0
                if B <= (a-A)/a
                    vN2 = -A/(B-1)
                elseif (a-A)/a < B
                    vN2 = a*B+A
                end
            elseif A <= a
                vN2 = a*B+A
            end
        end

        P2 = convert(Float64,vI2>vN2)

        vN1 = (r1(c,c2) + mp.β*(P2*H1(ic,ic,ic, pi, τ, ess)+(1-P2)*((1-pi)*Φ(τ[ic+1].EQs[ic, ic2, ess[ic,ic2,ic+1]].vN1,τ[ic+1].EQs[ic, ic2, ess[ic,ic2,ic+1]].vI1)+pi*Φ(0,-K(c, cost)))))/(1-mp.β*(1-P2)*pi)
        vI1 = vN1-K(c, cost)
        P1 = 0.

        @inbounds τ[ic].EQs[ic,ic2,1] = EQ(P1,
                                           vN1,
                                           vI1,
                                           P2,
                                           vN2,
                                           vI2)
    end
end

solve_interior!(cost::CostParameters, mp::RLSModel, ic::Int64, ic1::Int64, ic2::Int64, τ::Array{Stage,1}, ess::ESS) = find_interior!(mp, cost, ic1, ic2, ic, cost.C[ic], τ, ess)

function solve_interior!(cost::CostParameters, mp::RLSModel, ic::Int64, τ::Array{Stage,1}, ess::ESS)
    for ic2 in 1:ic-1
        for ic1 in 1:ic-1
            find_interior!(mp, cost, ic1, ic2, ic, cost.C[ic], τ, ess)
        end
    end
end

function find_interior!(mp, cost, ic1, ic2, ic, c, τ::Array{Stage,1}, ess::ESS)
    pi = cost.π[ic]

    ## First player 1 -> leads to P2 candidates
    a = r1(cost.C[ic1], cost.C[ic2]) - K(c, cost) +mp.β*H1(ic, ic2, ic, pi, τ, ess)
    b = mp.β*(H1(ic, ic, ic, pi, τ, ess)-H1(ic, ic2, ic, pi, τ, ess))
    d = r1(cost.C[ic1],cost.C[ic2]) + mp.β*(1-pi)*Φ(τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vN1,τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vI1)
    e = mp.β*H1(ic1, ic, ic, pi, τ, ess) - mp.β*(1-pi)*Φ(τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vN1,τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vI1)

    pa = - mp.β * pi * b
    pb = e + (mp.β*pi-1) * b - mp.β*pi*a
    pc = d + (mp.β*pi-1) * a

    pstar2 = quad(pa, pb, pc)

    ## Then player 2 -> leads to P1 candidates
    A = r2(cost.C[ic1], cost.C[ic2]) - K(c, cost) + mp.β*H2(ic1, ic, ic, pi, τ, ess)
    B = mp.β*(H2(ic, ic, ic, pi, τ, ess)-H2(ic1, ic, ic, pi, τ, ess))
    D = r2(cost.C[ic1],cost.C[ic2]) + mp.β*(1-pi)*Φ(τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vN2,τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vI2)
    E = mp.β*H2(ic, ic2, ic, pi, τ, ess)- mp.β*(1-pi)*Φ(τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vN2,τ[ic+1].EQs[ic1, ic2, ess[ic1,ic2,ic+1]].vI2)

    qa = - mp.β * pi * B
    qb = E + (mp.β*pi-1) * B - mp.β*pi*A
    qc = D + (mp.β*pi-1) * A

    pstar1 = quad(qa, qb, qc)

    count = 0
    @inbounds for i = 1:length(pstar1)
        for j = 1:length(pstar2)
            if i in 1:2 && j in 1:2
                exP1 = pc + pb * pstar2[j] + pa * pstar2[j]^2 < 0. ? 1. : 0.
                exP2 = qc + qb * pstar1[i] + qa * pstar1[i]^2 < 0. ? 1. : 0.

                if abs(exP1 - pstar1[i]) < 1e-7 && abs(exP2-pstar2[j]) < 1e-7
                    count += 1
                    @inbounds τ[ic].EQs[ic1, ic2, count] = EQ(pstar1[i],
                                                              (d + e*pstar2[j] + mp.β*pi*(1-pstar2[j])*(a+b*pstar2[j]))*pstar1[i]+(1-pstar1[i])*(d+e*pstar2[j])/(1-mp.β*pi*(1-pstar2[j])),
                                                              a + b*pstar2[j],
                                                              pstar2[j],
                                                              (D + E*pstar1[i] + mp.β*pi*(1-pstar1[i])*(A+B*pstar1[i]))*pstar2[j]+(1-pstar2[j])*(D+E*pstar1[i])/(1-mp.β*pi*(1-pstar1[i])),
                                                              A + B*pstar1[i])
                end
            elseif i > 2 && j > 2 && pstar1[i] >= 0 && pstar2[j] >= 0 && pstar1[i] <= 1 && pstar2[j] <= 1
                count += 1
                v1 = a + b * pstar2[j]
                v2 = A + B * pstar1[i]
                @inbounds τ[ic].EQs[ic1, ic2, count] = EQ(pstar1[i],
                                                          v1,
                                                          v1,
                                                          pstar2[j],
                                                          v2,
                                                          v2)
            end
        end
    end
#    #τ[ic].ne[ic1, ic2] = count
    @inbounds ess.bases[1 + div(ess.nC*(ess.nC+1)*(2*ess.nC+1),6) - div(ic*(ic+1)*(2*ic+1),6) + 2*(ic - 1) + sub2ind((ic-1,ic-1),ic1,ic2)] = count
end


text:
# check the Harmonic numbers
#    this is mostly just a check that code is correct
#    most errors are arising because of errors in the functions being used
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot

H = Symbol("real(H_n)")

##############################################
# integer arguments
data1 = CSV.read("../data/harmonic_test_data_1.csv", DataFrame)
m = size(data1,1)
H1 = zeros(Float64, m)
error1 = zeros(Float64, m)
rel_error1 = zeros(Float64, m)
for i=1:m
    H1[i] = harmonic( data1[i,:n] )
    error1[i] = H1[i] -  data1[i,H]
    rel_error1[i] = error1[i]/data1[i,H]
end
figure(1)
clf()
semilogy( data1[!,:n], abs.(rel_error1) )
xlabel("n")
ylabel("relative absolute error")
ylim( 1.0e-17, 1.0e-15)
title("Harmonic numbers")
savefig("plots/harmonic_bench_1.pdf")

##############################################
# real arguments
data2 = CSV.read("../data/harmonic_test_data_2.csv", DataFrame)
m = size(data2,1)
H2 = zeros(Float64, m)
error2 = zeros(Float64, m)
rel_error2 = zeros(Float64, m)
for i=1:m
    H2[i] = harmonic( data2[i,:x] )
    error2[i] = H2[i] -  data2[i,H]
    rel_error2[i] = error2[i]/data2[i,H]
end
figure(2)
clf()
# plot( data2[!,:x], rel_error2 )
semilogy( data2[!,:x], abs.(rel_error2) )
k2 = findall( abs.(rel_error2) .> 1.0e-13 )
plot( data2[k2,:x], abs.(rel_error2[k2]), "r" )
xlabel("x")
ylabel("relative absolute error")
title("Harmonic numbers (real arguments)")
savefig("plots/harmonic_bench_2.pdf")

##############################################
# complex arguments
zr = Symbol("Re(z)")
zi = Symbol("Im(z)")
Hr = Symbol("Re(H(z))")
Hi = Symbol("Im(H(z))")
data3 = CSV.read("../data/harmonic_test_data_3.csv", DataFrame)
m = size(data3,1)
H3 = zeros(Complex{Float64}, m)
error3 = zeros(Complex{Float64}, m)
rel_error3 = zeros(Float64, m)
for i=1:m
    z = data3[i,zr] + im*data3[i,zi]
    Hz = data3[i,Hr] + im*data3[i,Hi]
    # println(" i = $i, z = $z, typeof(z)=$(typeof(z)), H = $(Hz)" )
    H3[i] = harmonic( z )
    error3[i] = H3[i] - Hz
    rel_error3[i] = abs(error3[i])/ abs( Hz )
end
figure(3)
clf()
# plot( data3[!,:x], rel_error3 )
semilogy( atan.(data3[!,zi],data3[!,zr])/pi, abs.(rel_error3) )
k3 = findall( abs.(rel_error3) .> 1.0e-13 )
plot( data3[k3,:x], abs.(rel_error3[k3]), "r" )
xlabel("θ/π")
ylabel("relative absolute error")
title("Harmonic numbers (complex arguments)")
savefig("plots/harmonic_bench_3.pdf")

##############################################
# integer arguments for generalize H
H = Symbol("real(H_{n,r})")
data4 = CSV.read("../data/harmonic_test_data_4.csv", DataFrame)
m = size(data4,1)
H4 = zeros(Float64, m)
error4 = zeros(Float64, m)
rel_error4 = zeros(Float64, m)
for i=1:m
    H4[i] = harmonic( data4[i,:n], data4[i,:r] )
    error4[i] = H4[i] -  data4[i,H]
    rel_error4[i] = error4[i]/data4[i,H]
end
figure(4)
clf()
semilogy( data4[!,:n] + 0.1*data4[!,:r], abs.(rel_error4) )
xlabel("n + 0.1*r")
ylabel("relative absolute error")
ylim( 1.0e-17, 1.0e-15)
title("Generalised harmonic numbers")
savefig("plots/harmonic_bench_4.pdf")

##############################################
# integer n, real r for generalized H
H = Symbol("real(H_{n,r})")
data5 = CSV.read("../data/harmonic_test_data_5.csv", DataFrame)
m = size(data5,1)
H5 = zeros(Float64, m)
error5 = zeros(Float64, m)
rel_error5 = zeros(Float64, m)
for i=1:m
    H5[i] = harmonic( data5[i,:n], data5[i,:r] )
    error5[i] = H5[i] -  data5[i,H]
    rel_error5[i] = error5[i]/data5[i,H]
end
figure(5)
clf()
semilogy( data5[!,:n] + 0.1*data5[!,:r], abs.(rel_error5) )
xlabel("n + 0.1*r")
ylabel("relative absolute error")
ylim( 1.0e-17, 1.0e-15)
title("Generalised harmonic numbers")
savefig("plots/harmonic_bench_5.pdf")


text:
export generate

"""
    put_online(repo_name)

put a repository online: create a github and travis repository (if they don't
already exist) and connect them with a key.
"""
function put_online(repo_name; github_time = 60, travis_time = 60)
    if !endswith(repo_name, ".jl")
        ArgumentError("repo_name $repo_name must end with .jl")
    end
    github = GitHub(repo_name)
    travis = Travis(repo_name)
    if exists(github)
        @info "github already exists"
    else
        create(github)
        @info "Waiting $github_time seconds for github creation"
        sleep(github_time)
    end

    if exists(travis)
        @info "travis already exists"
    else
        user!(travis)
        sync(travis)
        @info "Waiting $travis_time seconds for travis syncing"
        sleep(travis_time)
        if !exists(travis)
            error("travis doesn't exist, likely due to incomplete syncing")
        end
    end
    repo!(travis)
    activate(travis)

    public_key, private_key = make_keys()

    delete_keys(github, ".documenter")
    add_key(github, ".documenter", public_key)
    delete_keys(travis, "DOCUMENTER_KEY")
    add_key(travis, "DOCUMENTER_KEY", private_key)

    github, travis
end


text:
using Documenter, Chatter

makedocs(
    modules=[Chatter],
    format=:html,
    pages=[
        "Home" => "index.md",
        "ChatterClient" => "pages/client.md",
        "Sources" => "pages/sources.md",
        "Index" => "pages/index.md",
    ],
    repo="https://github.com/christopher-dG/Chatter.jl/blob/{commit}{path}#L{line}",
    sitename="Chatter.jl",
    authors="Chris de Graaf",
    assets=[],
)

deploydocs(
    repo="github.com/christopher-dG/Chatter.jl",
    target="build",
    julia="0.6",
    deps=nothing,
    make=nothing,
)


text:
using ImageFeatures

function random_circle_example()
    dims = rand(100:200,2)
    img = zeros(dims...)
    r = rand(10:20)
    x,y = rand(20:80,2)
    for i in axes(img)[1]
        for j in axes(img)[2]
            if (x-i)^2 + (y-j)^2 <= r^2
                img[i,j] = 1
            end
        end
    end
    img, [x,y], r
end

@testset "Hough_Transform" begin
    @testset "Hough Line Transform" begin

    #For images containing a straight line parallel to axes
        for i in 1:9
            img = zeros(Bool,9,9)
            img[i,:] .= true
            h = hough_transform_standard(img)
            @test length(h) == 1
            line = first(h)
            @test line == (i, pi/2)
        end

    #For images with diagonal line
        img = Matrix(Diagonal([true, true ,true]))
        h = hough_transform_standard(img, angles=range(0,stop=pi,length=100))
        @test length(h) == 1
        @test h[1][1] == 0

    #For a square image
        img = zeros(Bool,10,10)
        for i in 1:10
            img[2,i] = img[i,2] = img[7,i] = img[i,9] = true
        end
        # h = hough_transform_standard(img,1,
        h = hough_transform_standard(img, angles=range(0,stop=π/2,length=100))
        @test length(h) == 4
        r = [h[i][1] for i in CartesianIndices(size(h))]
        @test all(r .== [2,2,7,9])
        theta = [h[i][2] for i in CartesianIndices(size(h))]
        er = sum(map((t1,t2) -> abs(t1-t2), theta, [0, π/2, π/2, 0]))
        @test er <= 0.1
    end

    @testset "Hough Circle Gradient" begin
        for _ in 1:10
            img, c_truth, r_truth = random_circle_example()
            img_edges = canny(img, (Percentile(80), Percentile(20)))
            dx, dy=imgradients(img, KernelFactors.ando5)
            img_phase = phase(dx, dy)
            radii = 10:20
            centers, rs = hough_circle_gradient(img_edges, img_phase, radii)
            @test length(centers) == length(rs) == 1
            c_hough = [Tuple(first(centers))...]
            r_hough = first(rs)
            @test r_hough ≈ r_truth atol=4
            @test c_hough ≈ c_truth atol=4
        end
    end
end


text:
using EntropicTree
using Test

@testset "EntropicTree.jl" begin
    # Write your tests here.
    @test  2+2 == 4
end


text:
######################################################################
# ergodic_manager/r2t.jl
#
# Euclidean (2d), with time-dependence
######################################################################

export ErgodicManagerR2T

# kpixl is a (K+1 x bins) matrix, each entry storing cos(k*pi*x / L)
#  this assumes some discretization
"""
`ErgodicManagerR2T(d::Domain, phi::Matrix{Float64}, K::Int)`
"""
mutable struct ErgodicManagerR2T <: ErgodicManager
	domain::Domain				# spatial domain
	K::Int						# number of Fourier coefficients
	phi::Array{Float64,3}		# spatial distribution
	phik::Array{Float64,3}		# distribution's Fourier coefficients

	# constant regardless of phi (depend on k1,k2)
	Lambda::Array{Float64,3}
	hk::Array{Float64,3}

	# to speed up computation
	kpixl::Matrix{Float64}
	kpiyl::Matrix{Float64}
	kpitl::Matrix{Float64}

	function ErgodicManagerR2T(d::Domain, phi::Array{Float64,3}, K::Int=5)
		em = new()
		em.domain = deepcopy(d)
		em.K = K
		em.hk = zeros(K+1,K+1,K+1)
		em.phi = deepcopy(phi)
		em.phik = zeros(K+1,K+1,K+1)
		em.Lambda = zeros(K+1,K+1,K+1)
		em.kpixl = zeros(K+1, d.cells[1])
		em.kpiyl = zeros(K+1, d.cells[2])
		N = size(phi, 3) - 1
		em.kpitl = zeros(K+1, N+1)

		Lambda!(em)
		kpixl!(em, N)
		hk!(em, N)
		decompose!(em)
		return em
	end

end

# fills each entry Lambda[k1,k2] in the Lambda matrix
# TODO: this is just repeated stuff, I should make this better
function Lambda!(em::ErgodicManagerR2T)
	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K
				den = (1.0 + k1*k1 + k2*k2 + k3*k3) ^ 2.0
				em.Lambda[k1+1, k2+1, k3+1] = 1.0 / den
			end
		end
	end
end

# TODO: check if the row/col ordering of this is julia-efficient
# TODO: also, I could probably speed up by iterating over k first?
function kpixl!(em::ErgodicManagerR2T, N::Int)
	Lx = em.domain.lengths[1]
	xmin = x_min(em)
	for xi = 1:x_cells(em)
		# TODO: seems like I just cancel xmin out. Ignore it maybe?
		x = xmin + (xi-0.5)*x_size(em)
		for k = 0:em.K
			em.kpixl[k+1,xi] = cos(k*pi*(x-xmin) / Lx)
		end
	end

	Ly = em.domain.lengths[2]
	ymin = y_min(em)
	for yi = 1:y_cells(em)
		y = ymin + (yi-0.5)*y_size(em)
		for k = 0:em.K
			em.kpiyl[k+1,yi] = cos(k*pi*(y-ymin) / Ly)
		end
	end

	# cos(k*pi*t/T) = cos(k*pi*(n*dt) / (N*dt)) = cos(k*pi*n / N)
	for n = 0:N
		for k = 0:em.K
			em.kpitl[k+1,n+1] = cos(k*pi*n / N)
		end
	end
end

# generates the hk coefficients for the ergodic manager
# these coefficients only need to be computed once
function hk!(em::ErgodicManagerR2T, N::Int)
	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K
				em.hk[k1+1,k2+1,k3+1] = hk_ij(em, k1, k2, k3, N)
			end
		end
	end
end

# computes the coefficients for a specific value of k1 and k2
# called by hk!
function hk_ij(em::ErgodicManagerR2T, k1::Int, k2::Int, k3::Int, N::Int)
	val = 0.0
	for xi = 1:x_cells(em)
		cx = em.kpixl[k1+1,xi]
		cx2 = cx * cx
		for yi = 1:y_cells(em)
			cy = em.kpiyl[k2+1,yi]
			cy2 = cy * cy
			for n = 0:N
				ct = em.kpitl[k3+1, n+1]
				ct2 = ct * ct
				val += cx2 * cy2 * ct2 * em.domain.cell_size
			end
		end
	end

	return sqrt(val / (N+1))
end

######################################################################
# Computing Fourier coefficients
######################################################################
# update the Fourier coefficients based on some distribution
# Here, I assume it is discrete, but I should change this...
# TODO: maybe I should do some bounds checking?
function decompose!(em::ErgodicManagerR2T, d::Array{Float64,3})
	N = size(d, 3) - 1
	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K
				em.phik[k1+1,k2+1,k3+1] = phi_ij(em, k1, k2, k3, d,N)/(N+1)
			end
		end
	end
	em.phi = d
end

# iterate over the state space
function phi_ij(em::ErgodicManagerR2T, k1::Int, k2::Int, k3::Int, d::Array{Float64,3}, N::Int)
	val = 0.0
	for xi = 1:x_cells(em)
		cx = em.kpixl[k1+1,xi]
		for yi = 1:y_cells(em)
			cy = em.kpiyl[k2+1,yi]
			for n = 0:N
				ct = em.kpitl[k3+1,n+1]
				val += d[xi,yi,n+1] * cx * cy * ct * em.domain.cell_size
			end
		end
	end
	return val / em.hk[k1+1,k2+1,k3+1]
end

function decompose(em::ErgodicManagerR2T, traj::VVF)
	N = length(traj)-1

	# Array to hold trajectory's Fourier coefficients
	ck = zeros(em.K+1, em.K+1, em.K+1)

	# length of each dimension
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]

	# minimum values in each dimension
	xmin = x_min(em)
	ymin = y_min(em)

	for k1 = 0:em.K
		kpiL1 = k1 * pi / Lx
		for k2 = 0:em.K
			kpiL2 = k2 * pi / Ly
			for k3 = 0:em.K
				kpiL3 = k3 * pi / N
				hk = em.hk[k1+1, k2+1, k3+1]
				fk_sum = 0.0
				# now loop over time
				for n = 0:N
					xn = traj[n + 1]
					c1 = cos(kpiL1 * (xn[1] - xmin))
					c2 = cos(kpiL2 * (xn[2] - ymin))
					c3 = cos(kpiL3 * n)
					fk_sum += c1*c2*c3
				end
				ck[k1+1, k2+1, k3+1] = fk_sum / (hk * (N+1))
			end
		end
	end
	return ck
end


text:
#!/bin/bash
#=
exec julia --color=yes --startup-file=no -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
=#

# @show ARGS  # put any Julia code here

##################################################
# 2021_08_17, 10_16, 10_29-30
# Compile spacer information for each genome, including level of self-targeting

# on RONIN, pwd: crispr-genome-analysis
# in Dropbox, pwd: genome-analysis

# usage: julia spacerinfo.jl
##################################################

##################################################
### Read spacer, dr, and Vink data ###
##################################################
using DelimitedFiles

spacerdata = readdlm("../CRISPRCasdb/2021/spacerdata_4.csv",',')
drdata = readdlm("../CRISPRCasdb/2021/drdata_4.csv",',')
Vinkdata = readdlm("Vinkdata.csv",',')[2:end,1:13]

#=
Parse Vinkdata
=#
seqIDs = String[]
spinds_seqs = Array{Int}[]

for i in 1:size(Vinkdata,1)
    IDs = strip.(split(Vinkdata[i,4],'|'))
    for id in IDs
        if in(id,seqIDs) # accession no. already present
            ind = findfirst(x->x==id,seqIDs)
            push!(spinds_seqs[ind],i)
        else # new accession no.
            push!(seqIDs,id)
            push!(spinds_seqs,[i])
        end
    end
end

##################################################
### Functions ###
##################################################
#=
Find the BLAST output row containing the closest match outside a CRISPR array to save
=#
function find_closestmatch(blastout, inds_s, arrayinfo, seqinfo)
    ind_tosave = 0
    len_max = 0
    for ind in inds_s
        m = blastout[ind,:]

        # 1) get no. of matches in alignment
        len_aln = m[4] - m[5]

        # 2) check if it's longer than the currently saved longest alignment
        if len_aln > len_max

            # 3) check if it's outside a CRISPR array
            if inarray(string(m[2]),sort([m[9],m[10]]),arrayinfo,seqinfo) == false
                len_max = len_aln
                ind_tosave = ind
            end
        end
    end

    return len_max, ind_tosave
end

#=
Returns true if alignment overlaps with *any* CRISPR array, false otherwise
Input:
    seq : accession no. of matching sequence
    pos : [start,end] of alignment on matching sequence (start < end)
    arrayinfo : information about CRISPR arrays in genome (rows of crisprdata_4.csv)
    seqinfo : information about sequences in genome
=#
function inarray(seq::String, pos::Array{Int}, arrayinfo, seqinfo)
    seqid = seqinfo[findfirst(x->x==seq,seqinfo[:,1]),4] # for identifying accession no. with sequence in CRISPRCasdb
    for i in 1:size(arrayinfo,1)
        if seqid == arrayinfo[i,2] && # array is on matching sequence
            length(intersect(pos[1]:pos[2],arrayinfo[i,3]:arrayinfo[i,3]+arrayinfo[i,4]-1)) > 0 # alignment overlaps with array
            ### println("spacer match to CRISPR array")
            return true
        end
    end
    return false
end

#=
Compute spacer category (0, 1, 2, or 3)
=#
function get_spcat(l_s::Int,len_max::Int)::Int
    if len_max == 0 || l_s - len_max > 10 # no self-target outside of CRISPR arrays (or with >10 mismatches)
        return 0
    elseif l_s - len_max > 3 # self-target with 4-10 mismatches
        return 1
    elseif l_s - len_max > 0 # self-target with 1-3 mismatches
        return 2
    else # exact self-target
        return 3
    end
end

##################################################
### Run through genomes and read array, sequence, and blast output info ###
##################################################

# # read genome
# g = "GCA_000006175.2"
# g = ARGS[1]

# read genomes
genomes = readdlm("dl.txt")

for g in genomes
    arrayinfo = readdlm("crisprs/$(g)_CRISPR.csv",',') # information about CRISPR arrays in genome (rows of crisprdata_4.csv)
    seqinfo = readdlm("seqinfo/$(g)_seqinfo.tsv",'\t') # information about sequences in genome
    # blastout = readdlm("test-STS/test.csv",',')
    blastout = readdlm("blast/$(g)_blast.csv",',')

    #=
    Parse spacer info (BLAST output col. 1); sid = [genbank, seqid, locus_start, spacer_start, spacer_length]
    =#
    sids_blastout = [split(m,'|') for m in blastout[:,1]]

    ##################################################
    ### Run through CRISPR arrays and compile spacer info ###
    ##################################################
    for i in 1:size(arrayinfo,1)
        # array information
        seqid_a = strip(arrayinfo[i,2],'"') # seqid
        start_a = arrayinfo[i,3] # array start
        drcons_a = strip(arrayinfo[i,6],'"') # dr consensus
        accnum_a = seqinfo[findfirst(x->x==seqid_a,seqinfo[:,4]),1] # sequence accession no.
        Vind_a = findfirst(x->x==accnum_a, seqIDs) # index in Vink seqIDs

        # all BLAST output rows corresponding to array
        inds_ainblastout = findall(x->x[2]==seqid_a && parse(Int,x[3])==start_a, sids_blastout)

        # all spacer IDs in array
        sids_a = unique(sids_blastout[inds_ainblastout])

        # find corresponding rows in spacerdata & drdata
        inds_s0 = findall(i->spacerdata[i,1]==g && spacerdata[i,2]==seqid_a && spacerdata[i,3]==start_a, 1:size(spacerdata,1))
        inds_d0 = findall(i->drdata[i,1]==g && drdata[i,2]==seqid_a && drdata[i,3]==start_a, 1:size(drdata,1))

        # run through each spacer in array
        for j in 1:length(sids_a)
            sid = sids_a[j]

            # spacer information
            start_s = parse(Int,sid[4])
            l_s = parse(Int,sid[5])
            seq_s = spacerdata[inds_s0[j],5]
            dr1_s = drdata[inds_d0[j],5]
            dr2_s = drdata[inds_d0[j+1],5]

            # obtain spacer information from Vinkdata
            # Vind_s = spinds_seqs[Vind_a][findfirst(x->x==seq_s, Vinkdata[spinds_seqs[Vind_a],2])]
            # subtype_s = Vinkdata[Vind_s,5]
            # orient_s = Vinkdata[Vind_s,12]
            # PAM_s = Vinkdata[Vind_s,13]

            if Vind_a != nothing && (Vind = findfirst(x->x==seq_s, Vinkdata[spinds_seqs[Vind_a],2])) != nothing
                Vind_s = spinds_seqs[Vind_a][Vind]
                subtype_s = Vinkdata[Vind_s,5]
                orient_s = Vinkdata[Vind_s,12]
                PAM_s = Vinkdata[Vind_s,13]
            else
                orient_s = ""
                subtype_s = ""
                PAM_s = ""
            end

            # all BLAST output rows corresponding to spacer
            inds_s = findall(x->x[2]==sid[2] && x[3]==sid[3] && x[4]==sid[4], sids_blastout)

            len_max, ind_tosave = find_closestmatch(blastout, inds_s, arrayinfo, seqinfo)

            c = get_spcat(l_s,len_max)

            ##################################################
            ### Output files ###
            ##################################################
            # append to spacer info file (1 line per spacer)
            open("spacerinfo/$(g)_spacerinfo.csv", "a") do io
                write(io, "$g,$i,$accnum_a,$start_a,$drcons_a,$subtype_s,$orient_s,$PAM_s,$j,$start_s,$l_s,$seq_s,$dr1_s,$dr2_s,$len_max,$c\n")
            end

            # append to new BLAST output file (1 line per spacer)
            open("blast1/$(g)_blast1.csv", "a") do io
                if c > 0
                    arr = string.(blastout[ind_tosave,:])
                    str = string([arr[i]*"," for i in 1:11]...,arr[end])
                    write(io, "$str\n")
                else
                    write(io, "$(blastout[inds_s[1],1]),,,,,,,,,,,\n")
                end
            end
        end
    end
end


text:
@testset "Simplicial Complex" begin
    cplx = SimplicialComplex(Simplex(1,2,3), Simplex(2,4), Simplex(4))

    @test celltype(cplx) == Simplex{Int}
    @test length(cells(cplx)) == 3
    cs = cells(cplx,2)
    @test cs !== nothing
    @test length(cs) == 1
    cs = cells(cplx,3)
    @test cs === nothing
    @test size(cplx) == (4,4,1)
    @test size(cplx,0) == 4
    @test size(cplx,1) == 4
    @test size(cplx,2) == 1
    @test size(cplx,3) == 0
    @test cplx[Simplex(4), 0] == 4
    @test cplx[Simplex(10), 0] > size(cplx,0)
    @test cplx[Simplex(10), 0] == 5

    @test length(simplices(cplx)) == 9
    scitr = simplices(cplx, 1)
    @test length(scitr) == 4
    for s in simplices(cplx, 2)
        @test s == Simplex(1,2,3)
    end

    cplx = SimplicialComplex(Simplex(1,2,3),
                       Simplex(2,4),
                       Simplex(3,4),
                       Simplex(5,4),
                       Simplex(6))

    qch = Chain(2, [1], [1])
    resch = Chain(1, [1, -1, 1], [3, 2, 1])
    @testset "Complex boundary chain" for (a,b) in zip(boundary(cplx, qch), resch)
        @test a == b
    end

    qch = Chain(1, [1], [1])
    resch = Chain(2, [-1, -1], [1, 2])
    @testset "Complex coboundary chain" for (a,b) in zip(coboundary(cplx, qch), resch)
        @test a == b
    end

    @test convert(Matrix, boundary_matrix(cplx, 2)) == [1 -1 1 0 0 0]'

    cplx = SimplicialComplex(Int)

    splx = push!(cplx, Simplex(1))
    @test splx[1][:index] == 1
    @test dim(splx[1]) == 0
    @test sum(size(cplx)) == 1
    @test cplx.cells[0][1][:values] == [1]

    splx = push!(cplx, Simplex(1,2))
    @test splx[1][:index] == 1
    @test dim(splx[1]) == 1
    @test sum(size(cplx)) == 2
    @test cplx.cells[1][1][:values] == [1,2]

    cplx = SimplicialComplex(Int)
    splxs = push!(cplx, Simplex(1,2), recursive=true)
    @test splxs[1][:index] == 1
    @test splxs[1][:values] == [1,2]
    @test splxs[2][:index] == 1
    @test splxs[2][:values] == [1]
    @test splxs[3][:index] == 2
    @test splxs[3][:values] == [2]
    @test sum(size(cplx)) == 3
    @test size(cplx, 0) == 2
    @test size(cplx, 1) == 1

    splxs = push!(cplx, Simplex(1,3), recursive=true)
    @test splxs[1][:index] == 2
    @test splxs[1][:values] == [1,3]
    @test splxs[2][:index] == 3
    @test splxs[2][:values] == [3]
    @test sum(size(cplx)) == 5
    @test size(cplx, 0) == 3
    @test size(cplx, 1) == 2

    io = IOBuffer()
    write(io, cplx)
    seekstart(io)
    cplx2 = read(io, SimplicialComplex{Int})
    @test size(cplx2, 0) == 3
    @test size(cplx2, 1) == 2

    ch = boundary(cplx2, 1, 0, Int)
    @test dim(ch) == -1
end


text:
"""
    NeXLMatrixCorrection.quantify(
      ffr::FitResult;
      strip::AbstractVector{Element} = [],
      mc::Type{<:MatrixCorrection} = XPP,
      fl::Type{<:FluorescenceCorrection} = ReedFluorescence,
      kro::KRatioOptimizer = SimpleKRatioOptimizer(1.5),
    )::IterationResult

Facilitates converting `FilterFitResult` or `BasicFitResult` objects into estimates of composition by extracting
k-ratios from measured spectra and applying matrix correction algorithms.
"""
function NeXLMatrixCorrection.quantify(
    ffr::FitResult;
    strip::AbstractVector{Element} = Element[],
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
    kro::KRatioOptimizer = SimpleKRatioOptimizer(1.5),
    unmeasured::UnmeasuredElementRule = NullUnmeasuredRule(),
    coating::Union{Nothing, Pair{CharXRay, <:Material}}=nothing
)::IterationResult
    iter = Iteration(mc, fc, cc, unmeasured = unmeasured)
    krs = optimizeks(kro, filter(kr -> !(element(kr) in strip), kratios(ffr)))
    return quantify(iter, ffr.label, krs, coating=coating)
end

"""
    NeXLMatrixCorrection.quantify(
        spec::Union{Spectrum,AbstractVector{<:Spectrum}},
        ffp::FilterFitPacket;
        kwargs...
    )::IterationResult

Failitates quantifying spectra.  First filter fits and then matrix corrects.
"""
function NeXLMatrixCorrection.quantify(
    spec::Spectrum,
    ffp::FilterFitPacket;
    kwargs...,
)::IterationResult
    return quantify(fit_spectrum(spec, ffp); kwargs...)
end
function NeXLMatrixCorrection.quantify(
    specs::AbstractVector{<:Spectrum},
    ffp::FilterFitPacket;
    kwargs...,
)::Vector{IterationResult}
    return map(spec->quantify(fit_spectrum(spec, ffp); kwargs...), specs)
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum