text:
# Test 2D modeling
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, JUDI.TimeModeling, SeisIO

## Set up model structure
n = (120,100)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v0 = ones(Float32,n) .+ 0.4f0
v[:,Int(round(end/2)):end] .= 4.f0

rho = ones(Float32,n)
rho[:,Int(round(end/2)):end] .= 1.2f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 2	# number of sources
model = Model(n,d,o,m; rho=rho)	# to include density call Model(n,d,o,m,rho)
model0 = Model(n,d,o,m0; rho=rho)

## Set up receiver geometry
nxrec = 60
xrec = range(50f0,stop=500f0,length=nxrec)
yrec = 0f0
zrec = range(50f0,stop=50f0,length=nxrec)

# receiver sampling and recording time
timeR = 1000f0	# receiver recording time [ms]
dtR = 2f0	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(100f0,stop=20f0,length=nsrc))
ysrc = convertToCell(range(0f0,stop=0f0,length=nsrc))
zsrc = convertToCell(range(20f0,stop=20f0,length=nsrc))

# source sampling and number of time steps
timeS = 1000.
dtS = 2.	# receiver sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.01
wavelet = ricker_wavelet(timeS,dtS,f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)
info = Info(prod(n),nsrc,ntComp)

######################## WITH DENSITY ############################################

# Keep data in memory
println("Test modeling without density, data in RAM")

opt = Options(limit_m = true,
              buffer_size = 100f0
			  )

# Setup operators
Pr = judiProjection(info,recGeometry)
F = judiModeling(info,model;options=opt)
F0 = judiModeling(info,model0;options=opt)
Ps = judiProjection(info,srcGeometry)
q = judiVector(srcGeometry,wavelet)

# Combined operator Pr*F*adjoint(Ps)
Ffull = judiModeling(info,model,srcGeometry,recGeometry)

# Nonlinear modeling
d1 = Pr*F*adjoint(Ps)*q	# equivalent to d = Ffull*q
qad = Ps*adjoint(F)*adjoint(Pr)*d1

# Linearized modeling
J = judiJacobian(Pr*F0*adjoint(Ps),q)	# equivalent to J = judiJacobian(Ffull,q)
dD = J*dm
rtm = adjoint(J)*dD

# fwi objective function
f,g = fwi_objective(model0,q,d1)
f,g = fwi_objective(model0,subsample(q,1),subsample(d1,1))

# Subsampling
dsub1 = subsample(d1,1)
dsub2 = subsample(d1,[1,2])
Fsub1 = subsample(F,1)
Fsub2 = subsample(F,[1,2])
Jsub1 = subsample(J,1)
Jsub2 = subsample(J,[1,2])
Ffullsub1 = subsample(Ffull,1)
Ffullsub2 = subsample(Ffull,[1,2])
Psub1 = subsample(Pr,1)
Psub2 = subsample(Pr,[1,2])

# vcat, norms, dot
dcat = [d1,d1]
norm(d1)
dot(d1,d1)

#########################################################################
# Save data to segy files
println("Test modeling, save data to disk")

# Options structures
opt = Options(limit_m = true,
              buffer_size = 100f0,
              save_data_to_disk=true,
			  file_path=pwd(),	# path to files
			  file_name="shot_record"	# saves files as file_name_xsrc_ysrc.segy
			  )

opt0 = Options(limit_m = true,
              buffer_size = 100f0,
              save_data_to_disk=true,
			  file_path=pwd(),	# path to files
			  file_name="linearized_shot_record"	# saves files as file_name_xsrc_ysrc.segy
			  )

# Setup operators
F = judiModeling(info,model;options=opt)
F0 = judiModeling(info,model0;options=opt0)

# Nonlinear modeling
dsave = Pr*F*adjoint(Ps)*q	# equivalent to d = Ffull*q
qad = Ps*adjoint(F)*adjoint(Pr)*d1

# Linearized modeling
J = judiJacobian(Pr*F0*adjoint(Ps),q)	# equivalent to J = judiJacobian(Ffull,q)
dDsave = J*dm
rtm = adjoint(J)*dD

# Subsampling
dsub1 = subsample(d1,1)
dsub2 = subsample(d1,[1,2])
Fsub1 = subsample(F,1)
Fsub2 = subsample(F,[1,2])
Jsub1 = subsample(J,1)
Jsub2 = subsample(J,[1,2])
Ffullsub1 = subsample(Ffull,1)
Ffullsub2 = subsample(Ffull,[1,2])
Psub1 = subsample(Pr,1)
Psub2 = subsample(Pr,[1,2])

# fwi objective function
f,g = fwi_objective(model0,q,dsave)
f,g = fwi_objective(model0,subsample(q,1),subsample(dsave,1))

# vcat, norms, dot
dcat = [d1,d1]
norm(d1)
dot(d1,d1)

rm("shot_record_100.0_0.0.segy")
rm("shot_record_20.0_0.0.segy")
rm("linearized_shot_record_100.0_0.0.segy")
rm("linearized_shot_record_20.0_0.0.segy")


text:
# Utility -------

function github_obj_from_type(data::Dict)
    t = get(data, "type", nothing)

    if t == "User"
        return User(data)
    elseif t == "Organization"
        return Organization(data)
    end
end

ext: jl


text:
"""
    readMapFromFITS{T <: Number}(f::FITSIO.FITSFILE, column, t::Type{T})
    readMapFromFITS{T <: Number}(fileName::String, column, t::Type{T})

Read a Healpix map from the specified (1-base indexed) column in a
FITS file. The values will be read as numbers of type T. If the code
fails, FITSIO will raise an exception. (Refer to the FITSIO library
for more information.)
"""

function readMapFromFITS(f::FITSIO.FITSFile, column, t::Type{T}) where
                        {T <: Number}
    value, comment = FITSIO.fits_read_keyword(f, "NSIDE")
    nside = parse(Int, value)

    value, comment = FITSIO.fits_read_keyword(f, "ORDERING")
    ringOrdering = uppercase(strip(value[2:end-1])) == "RING"

    repeat = (FITSIO.fits_get_coltype(f, column))[2]
    nrows = FITSIO.fits_get_num_rows(f)

    if repeat * nrows != nside2npix(nside)
        error("Wrong number of pixels in column $column of FITS file (NSIDE=$nside)")
    end

    if ringOrdering
        result = Map{T, Ring}(Array{T}(undef, nside2npix(nside)))
    else
        result = Map{T, Nest}(Array{T}(undef, nside2npix(nside)))
    end
    FITSIO.fits_read_col(f, column, 1, 1, result.pixels)

    result
end

function readMapFromFITS(fileName::AbstractString,
                         column,
                         t::Type{T}) where {T <: Number}
    f = FITSIO.fits_open_table(fileName)
    result = readMapFromFITS(f, column, t)
    FITSIO.fits_close_file(f)

    result
end

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

"""
    savePixelsToFITS(map::Map{T}, f::FITSIO.FITSFile, column) where {T <: Number}

Save the pixels of `map` into the column with index/name `column` in the FITS
file, which must have been already opened.
"""
function savePixelsToFITS(map::AbstractMap{T, O},
                          f::FITSIO.FITSFile,
                          column) where {T <: Number, O}

    FITSIO.fits_update_key(f, "PIXTYPE", "HEALPIX",
                           "HEALPIX pixelisation")
    FITSIO.fits_update_key(f, "NSIDE", map.resolution.nside,
                           "Value of NSIDE")
    FITSIO.fits_update_key(f, "FIRSTPIX", 1,
                           "First pixel (1 based)")
    FITSIO.fits_update_key(f, "LASTPIX", map.resolution.numOfPixels,
                           "Last pixel (1 based)")
    FITSIO.fits_update_key(f, "INDXSCHM", "IMPLICIT",
                           "Indexing: IMPLICIT or EXPLICIT")
    FITSIO.fits_write_col(f, column, 1, 1, map.pixels)

end

"""
    saveToFITS{T <: Number, O <: Order}(map::Map{T, O},
                                        f::FITSIO.FITSFile,
                                        column)
    saveToFITS{T <: Number, O <: Order}(map::Map{T, O},
                                        fileName::String,
                                        typechar="D",
                                        unit="",
                                        extname="MAP")

Save a Healpix map in the specified (1-based index) column in a FITS
file. If the code fails, FITSIO will raise an exception. (Refer to the
FITSIO library for more information.)
"""
function saveToFITS(map::OrderedMap{T, Ring},
                    f::FITSIO.FITSFile,
                    column) where {T <: Number}

    FITSIO.fits_update_key(f, "ORDERING", "RING")
    savePixelsToFITS(map, f, column)

end

function saveToFITS(map::OrderedMap{T, Nest},
                    f::FITSIO.FITSFile,
                    column) where {T <: Number}

    FITSIO.fits_update_key(f, "ORDERING", "NEST")
    savePixelsToFITS(map, f, column)

end

"""
    saveToFITS(map::Map{T, O}, filename::AbstractString, typechar="D", unit="", extname="MAP") where {T <: Number, O <: Order}

Save a map into a FITS file. The name of the file is specified in
`filename`; if it begins with `!`, existing files will be overwritten
without warning. The parameter `typechar` specifies the data type to
be used in the FITS file: the default (`D`) will save 64-bit
floating-point values. See the CFITSIO documentation for other
values. The keyword `unit` specifies the measure unit used for the
pixels in the map. The keyword `extname` specifies the name of the HDU
where the map pixels will be written.

"""
function saveToFITS(map::OrderedMap{T, O},
                    fileName::AbstractString;
                    typechar="D",
                    unit="",
                    extname="MAP") where {T <: Number, O <: Order}

    f = FITSIO.fits_create_file(fileName)
    try
        FITSIO.fits_create_binary_tbl(f, 0, [("PIXELS", "1$typechar", unit)], extname)
        saveToFITS(map, f, 1)
    finally
        FITSIO.fits_close_file(f)
    end

end

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


text:
function gradient_penalty(D, x; target::Float32=1f0)
	B = size(x, ndims(x))
	l, b = Flux.pullback(D, x)
	ones_vec = Zygote.ignore() do
					ones(Float32, size(l)) |> device(x)
				end
	grads = b(ones_vec)
    Flux.mean((sqrt.(sum(reshape(grads[1], :, B).^2, dims = 1)) .- target).^2)
end

function gradient_penalty(D, x, xtilde; target::Float32=1f0)
	B = size(x, ndims(x))
	ϵ = Zygote.ignore(() -> Float32.(rand(Uniform(), 1, B)) |> device(x))
	ϵx = reshape(ϵ, ones(Int,  ndims(x)-1)..., B)
	xhat = ϵx .* xtilde .+ (1f0 .- ϵx) .* x
	gradient_penalty(D, xhat, target=target)
end


text:
using MixedModels, SelectMixedModels, RDatasets

ds = dataset("lme4", "Dyestuff");
fm1 = fit(lmm(Yield ~ 1|Batch, ds))
@time conAIC(fm1)
marAIC(fm1)
marBIC(fm1)

sleep = dataset("lme4", "sleepstudy");
fm2 = fit(lmm(Reaction ~ Days + (Days | Subject), sleep))
@time conAIC(fm2)
marAIC(fm2)
marBIC(fm2)

inst = dataset("lme4", "InstEval");
fm3 = fit(lmm(Y ~ Dept*Service + (1|S) , inst))
@time conAIC(fm3)
marAIC(fm3)
marBIC(fm3)

chem = dataset("mlmRev", "Chem97");
fm4 = fit(lmm(Score ~  (1|Lea), chem))
@time conAIC(fm4)
marAIC(fm4)
marBIC(fm4)

pen = dataset("lme4", "Penicillin");
fm5 = fit(lmm(Diameter ~ (1|Plate) + (1|Sample), pen))
@time conAIC(fm5)
marAIC(fm5)
marBIC(fm5)

psts = dataset("lme4", "Pastes");
fm6 = fit(lmm(Strength ~ (1|Sample) + (1|Batch), psts))
@time conAIC(fm6)
marAIC(fm6)
marBIC(fm6)

fm7 = fit(lmm(Score ~  (1|School) + (1|Lea), chem))
@time conAIC(fm7)
marAIC(fm7)
marBIC(fm7)

fm8 = fit(lmm(Y ~ Dept*Service + (1|S) + (1|D), inst))
@time conAIC(fm8)
marAIC(fm8)
marBIC(fm8)

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
function mktestseis()
  L0 = 30
  L1 = 10
  os = 5
  tt = time()
  t1 = round(Int64, tt/μs)
  t2 = round(Int64, (L0+os)/μs) + t1

  S = SeisData(5)
  S.name = ["Channel 1", "Channel 2", "Channel 3", "Longmire", "September Lobe"]
  S.id = ["XX.TMP01.00.BHZ","XX.TMP01.00.BHN","XX.TMP01.00.BHE","CC.LON..BHZ","UW.SEP..EHZ"]
  S.fs = collect(Main.Base.Iterators.repeated(100.0, S.n))
  S.fs[4] = 20.0
  for i = 1:S.n
    os1 = round(Int64, 1/(S.fs[i]*μs))
    S.x[i] = randn(Int(L0*S.fs[i]))
    S.t[i] = [1 t1+os1; length(S.x[i]) 0]
  end

  T = SeisData(4)
  T.name = ["Channel 6", "Channel 7", "Longmire", "September Lobe"]
  T.id = ["XX.TMP02.00.EHZ","XX.TMP03.00.EHN","CC.LON..BHZ","UW.SEP..EHZ"]
  T.fs = collect(Main.Base.Iterators.repeated(100.0, T.n))
  T.fs[3] = 20.0
  for i = 1:T.n
    T.x[i] = randn(Int(L1*T.fs[i]))
    T.t[i] = [1 t2; length(T.x[i]) 0]
  end
  return (S,T)
end


text:
include("../src/hypergraph-affinity-functions.jl")
using MAT

using Plots
M = matread("Total_Ns.mat")
Ng = M["Ng"]
Nf = M["Nf"]
Nw = M["Nw"]
r = 10

## Run all
for j = 1:3
for k = 2:4
if  j == 1
    N = Nf
elseif j == 2
    N = Nw
else
    N = Ng
end

alpha = 1/2
Hw, Hm = affinities_from_N(N)
alpha = 1/2
Bw = baselines(r,alpha)
Bm = baselines(r,1-alpha)
R1 = Hw./Bw
R2 = Hm./Bm

## Counting balance between men and women

N = Nw
T = 2*sum(N[2,:]) + 3*sum(N[3,:]) + 4*sum(N[4,:])
women = N[2,2] + 2*N[2,3] + N[3,2] + 2*N[3,3] + 3*N[3,4] + N[4,2] + 2*N[4,3] + 3*N[4,4] + 4*N[4]
@show women/T

## Plot
s1 = 375
s2 = 300
gfs = 18
tfs = 18
titlesize = 18
color1 = :green
color2 = :blue
title = ""
leg = false
if j == 3
    xlab = "Affinity type t"
else
    xlab = ""
end

if k == 2
    yl = [.2,1.7]
elseif k ==3
    yl = [.2,2]
elseif k == 4
    yl = [.1,2.5]
end

if k == 2
    ylab = "Affinity / Baseline"
else
    ylab = ""
end
p = Plots.plot([1; k],[1; 1],yscale = :identity,fontfamily = "helvetica",linewidth = 2,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = 1:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs,ylim = yl)
plot!(p,1:k,R1[k,1:k],color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)
plot!(p,1:k,R2[k,1:k],color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)

savefig("Figures/Ratios_$(k)_$(j).pdf")

end

end

## Compute pairs in the clique projection
N = Nw
In1 = 0
In2 = 0
Cross = 0
r = 10
for k = 2:r
    global In1, In2, Cross
    for j = 1:k+1
        i = j-1         # number of women in this type hyperedge
        Hik = N[k,j]    # number of hyperedges of this type
        In1 += binomial(i,2)*Hik    # number of extra edges of this type
        In2 += binomial(k-i,2)*Hik
        Cross += Hik*i*(k-i)
        # println("$k, $i")
    end
end

hw = round(2*In1/(2*In1 + Cross),digits = 3)
hm = round(2*In2/(2*In2 + Cross),digits = 3)
println("Clique expansion homophily scores: women = $hw, men = $hm")
println("$hw & $hm")

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


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

using Markdown
using InteractiveUtils

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

# ╔═╡ ae24c8b2-f60b-11ea-2c7a-03857d1217b2
using Pkg

# ╔═╡ bc14fc1a-f60b-11ea-207a-91b967f28076
begin
	pkg"add Colors ColorSchemes Images ImageMagick PlutoUI Suppressor InteractiveUtils"
	using Colors, ColorSchemes, Images, ImageMagick
	using Suppressor, InteractiveUtils, PlutoUI
end

# ╔═╡ 0feb5674-f5d5-11ea-0714-7379a7d381a3
using LinearAlgebra

# ╔═╡ 20125236-f5d6-11ea-3877-6b332497be62
using SparseArrays

# ╔═╡ 0f6ecba6-f5da-11ea-22c2-2929e562f413
using Statistics

# ╔═╡ b0ba5b8c-f5d1-11ea-1304-3f0e47f935fe
md"# What do data structures, machine learning and the singular-value decomposition have in common?"

# ╔═╡ ffa95430-f5d1-11ea-3cb7-5bb8d8f13701
md"""Blurry Fons magic: Exploiting Structure!

ADD A HINT BOX HERE WITH THE ANSWER "STRUCTURE" """

# ╔═╡ 261c4df2-f5d2-11ea-2c72-7d4b09c46098
md"One-hot vectors:  Numbers and images

0s and 1s and colors"

# ╔═╡ b5177f70-f60b-11ea-14a9-f5a574cc5185
Pkg.activate(mktempdir())

# ╔═╡ ef8f44b2-f5fc-11ea-1e4d-bd873cd39d6c
@bind nn Slider(1:20, show_value=true)

# ╔═╡ fd9211c0-f5fc-11ea-1745-7f2dae88af9e
@bind kk Slider(1:nn, default=nn, show_value=true)

# ╔═╡ 3a1ed5b8-f5fd-11ea-2ecd-b9d08349651f
colors = [colorant"black", colorant"red"]

# ╔═╡ 77ae1146-f5d2-11ea-1226-27d15c90a8df
begin
	v = fill(1, nn)
	v[kk] = 2

	colors[v]'
end

# ╔═╡ 676f6e3a-f5fd-11ea-3552-e7e7c6382276
v

# ╔═╡ 3592648a-f5fd-11ea-2ed0-a1d4a4c14d07

# ╔═╡ 7ffa02f6-f5d2-11ea-1b41-15c21f92e9ee
Int.((1:10) .== 3)

# ╔═╡ 21a15906-f5d3-11ea-3a71-f53eabc31acc
md"Colors and sliders"

# ╔═╡ 3c5f0ba8-f5d3-11ea-22d3-4f7c513144bc
md"Data structure for a one-hot vector"

# ╔═╡ 4d2db24e-f5d4-11ea-10d4-b126f643abee

# ╔═╡ 4624cd26-f5d3-11ea-1cf8-7d6555eae5fa
struct OneHot <: AbstractVector{Int}
	n::Int
	k::Int
end

# ╔═╡ 397ac764-f5fe-11ea-20cc-8d7cab19d410
Base.size(x::OneHot) = (x.n, )

# ╔═╡ 82c7046c-f5d3-11ea-04e2-ef7c0f4db5da
Base.getindex(x::OneHot, i::Int) = Int(x.k == i)

# ╔═╡ 7a98f708-f5d3-11ea-2301-cde71ca469ff
x = OneHot(5, 3)

# ╔═╡ 60214c1a-f5fe-11ea-2d08-c59715bcedd0
x[1:2]

# ╔═╡ cf77d83c-f5fd-11ea-136e-8951de64e05e
show_vector(x::OneHot) = colors[x .+ 1]'

# ╔═╡ 5813e1b2-f5ff-11ea-2849-a1def74fc065
begin
	imshow(M) = get.(Ref(ColorSchemes.rainbow), M ./ maximum(M))
	imshow(x::AbstractVector) = imshow(x')
end

# ╔═╡ e51d4ef6-f5fd-11ea-28f1-a1616e17f715
imshow(rand(5, 5))

# ╔═╡ 982590d4-f5ff-11ea-3802-73292c75ad6c
imshow(x)

# ╔═╡ c3f70e18-f5ff-11ea-35aa-31b22ca506b8
get.([ColorSchemes.rainbow], x)'

# ╔═╡ c200deb0-f5fd-11ea-32e1-d96596b16ebc
x[3]

# ╔═╡ a5989a76-f5d3-11ea-3fe5-5f9959199fe8
typeof(x)

# ╔═╡ b0bd6c6a-f5d3-11ea-00ee-b7c155f23481
x[3:5]

# ╔═╡ bd94f82c-f5d3-11ea-2a1e-77ddaadbfeb9
collect(x)

# ╔═╡ 8d2c6910-f5d4-11ea-1928-1baf09815687
md"""How much "information" (numbers) do you need to represent a one-hot vector? Is it $n$ or is it two?
"""

# ╔═╡ 81c35324-f5d4-11ea-2338-9f982d38732c
md"# Diagonal matrices"

# ╔═╡ 2cfda0dc-f5d5-11ea-16c4-b5a33b90e37f
md"As you might see in high school:"

# ╔═╡ 150432d4-f5d5-11ea-32b2-19a2a91d9637
M = [5 0 0
	 0 6 0
	 0 0 -10]

# ╔═╡ 21328d1c-f5d5-11ea-288e-4171ad35326d
Diagonal(M)

# ╔═╡ d8f36278-f5d5-11ea-3338-8573ce40e65e
md"How much information do you need for a diagonal matrix?"

# ╔═╡ e90c55fc-f5d5-11ea-10f1-470ff772985d
md"""We should always look for *structure* where it exists"""

# ╔═╡ 19775c3c-f5d6-11ea-15c2-89618e654a1e
md"## Sparse matrices"

# ╔═╡ 232d1dcc-f5d6-11ea-2710-658c75b9c7a4
sparse(M)

# ╔═╡ 2b2feb9c-f5d6-11ea-191c-df20360b12a9
M2 = spzeros(10^6, 10^6)

# ╔═╡ 50c74f6c-f5d6-11ea-29f1-5de9997d5d9f
M2[1, 2] = 1

# ╔═╡ 5de72b7c-f5d6-11ea-1b6f-35b830b5fb34
M3 = sparse([1, 2, 10^6], [4, 10, 10^6], [7, 8, 9])

# ╔═╡ 2fd7e52e-f5d7-11ea-3b5a-1f338e2451e0
M4 = [1 0 2 0 10; 0 3 4 0 9; 0 0 0 5 8; 0 0 0 0 7] .* 0

# ╔═╡ aa09c008-f5d8-11ea-1bdc-b51ee6eb2478
sparse(M4)

# ╔═╡ d941cd66-f5d8-11ea-26ff-47ba7779ab20
md"Sparse matrix dictionaries? Before CSC"

# ╔═╡ 62a6ec62-f5d9-11ea-071e-ed33c5dea0cd
md"## Rand: Where is the structure?"

# ╔═╡ 67274c3c-f5d9-11ea-3475-c9d228e3bd5a
# join(rand("ACGT", 100))

vv = rand(1:9, 1000000)

# ╔═╡ b6c7a918-f600-11ea-18ff-6521507358c6
md"Mention lossless compression e.g. run-length encoding"

# ╔═╡ 765c6552-f5d9-11ea-29d3-bfe7b4b04612
md"""Some might guess that there is "no structure"

Take mean and standard deviation -- some would say that's the structure

"""

# ╔═╡ 126fb3ea-f5da-11ea-2f7d-0b3259a296ce
mean(vv), std(vv), 5, sqrt(10 * 2/3)

# ╔═╡ 5f79e8f4-f5da-11ea-2b55-ef344b8a3ba2
var(vv)

# ╔═╡ 9b9e2c2a-f5da-11ea-369b-b513b196515b
md"Statisticians (and professors who've just graded exams) might say that under certain circumstances the mean and the variance give you the necessary structure, and the rest can be thrown away"

# ╔═╡ e68b98ea-f5da-11ea-1a9d-db45e4f80241
m = sum(vv) / length(vv)  # mean

# ╔═╡ f20ccac4-f5da-11ea-0e69-413b5e49f423
σ² = sum( (vv .- m) .^ 2 ) / (length(vv) - 1)

# ╔═╡ 0bc792e8-f5db-11ea-0b7a-1502ddc8008e
var(vv)

# ╔═╡ 12a2e96c-f5db-11ea-1c3e-494ae7446886
σ = sqrt(σ²)

# ╔═╡ 22487ce2-f5db-11ea-32e9-6f70ab2c0353
std(vv)

# ╔═╡ 389ae62e-f5db-11ea-1557-c3adbbee0e5c
md"Sometimes the summary statistics are all you want. (But sometimes not)"

# ╔═╡ 0c2b6408-f5d9-11ea-2b7f-7fece2eecc1f
md"## Multiplication tables"

# ╔═╡ 542a9556-f5db-11ea-0375-99f52416f6e4
md"How do you make a multiplication table?"

# ╔═╡ 165788b2-f601-11ea-3e69-cdbbb6558e54
md"Do you recognise this?"

# ╔═╡ 22941bb8-f601-11ea-1d6e-0d955297bc2e
rand(3) .* rand(4)'

# ╔═╡ 2f75df7e-f601-11ea-2fc2-aff4f335af33
imshow(rand(3) .* rand(4)')

# ╔═╡ 53e6b612-f601-11ea-05a9-5395e69b3c41
svd(rand(3) .* rand(4)')

# ╔═╡ 3e919766-f601-11ea-0485-05f45484bf8d
md"It's not easy to see the structure!"

# ╔═╡ 1052993e-f601-11ea-2c55-0d67e31b670e
M5 =
[1  2  3   4
 2  4  6   8
 3  6  9  12
]

# ╔═╡ 68190822-f5db-11ea-117f-d10a161208c3
md"Comprehension:"

# ╔═╡ 71b44874-f5db-11ea-1f67-47bad9295e03
md"Slider"

# ╔═╡ 6c51eddc-f5db-11ea-1235-332cdbb072fa
[i * j for i in 1:3, j in 1:4]

# ╔═╡ 86fb49ee-f5db-11ea-3bfa-c95c3b8775a3
md"Explain '"

# ╔═╡ 173cfab4-f5d9-11ea-0c7c-bf8b0888f6e7
(1:3) .* (1:4)'

# ╔═╡ 8bab3e36-f5db-11ea-187a-f31fa8cf357d
reshape([1, 2, 3], 3, 1)

# ╔═╡ e64291e8-f5db-11ea-0cab-8567b781408f
[1
 2
 3]

# ╔═╡ 0111a124-f5dc-11ea-0904-fdd88d7acac4
MM = rand(3, 3)

# ╔═╡ 06406158-f5dc-11ea-02b1-2519a0176993
MM[:, 1]

# ╔═╡ 12051f06-f5dc-11ea-2f6c-0fdc50eeff01
vvv = [1, 2, 3]

# ╔═╡ 195854e6-f5dc-11ea-114b-2333f87173f7
a = reshape([1, 2, 3], 3, 1)

# ╔═╡ 768116f6-f5dc-11ea-1cfe-b3016c574725
b = reshape([4, 5, 6], 3, 1)

# ╔═╡ 9050426e-f5dc-11ea-373c-65456732bd34
a * b'

# ╔═╡ 965bd07e-f5dc-11ea-2e85-d34996cf2fae
v4 = [1, 2, 3]

# ╔═╡ a9daf4cc-f5dc-11ea-270b-2566f89f168c
v4 * v4'

# ╔═╡ dc0c8b72-f5dc-11ea-3e6f-0f43cbf58f56
v4 .* v4'

# ╔═╡ 0ed78e76-f5dd-11ea-2ad8-a35a69c0ef9a
v4

# ╔═╡ 1648a0fa-f5dd-11ea-0292-495207e83de9
[v4 v4]

# ╔═╡ 1e9cefea-f5dd-11ea-3e5f-a189fd41c42e
[1*v4 5*v4 6*v4]

# ╔═╡ 24664f20-f5dd-11ea-2a69-cd3e0ebd5c39
v4 * [1, 5, 6]'

# ╔═╡ 43a4920c-f5dd-11ea-1ab1-0b3d673c0f1e
v4 .* [1, 5, 6]'

# ╔═╡ 5d767290-f5dd-11ea-2189-81198fd216ce
outer(v, w) = [x * y for x in v, y in w]  # just a multiplication table

# ╔═╡ 5a493052-f601-11ea-2f5f-f940412905f2
begin
	v6 = rand(3)
	w6 = rand(4)

	U6, Σ6, V6 = svd( outer(v6, w6) )
end

# ╔═╡ 8633afb2-f601-11ea-206b-e9c4b9621c2a
outer(U6[:,1], V6[:, 1]) .* Σ6[1]

# ╔═╡ 9918c4fa-f601-11ea-3bf1-3506dcb437f7
outer(v6, w6)

# ╔═╡ 6aae805e-f5dd-11ea-108c-733daae313dc
outer(v4, [1, 5, 6])

# ╔═╡ 9a023cf8-f5dd-11ea-3016-f95d433e6df0
outer(1:10, 1:10)  # works with things other than Vectors

# ╔═╡ b4c82246-f5dd-11ea-068f-2f63a5a382e2
md"Did you memorize this in third grade?"

# ╔═╡ d1f87b22-f5dd-11ea-3bc3-471d5b3a5202
md"Slider: 12 is the standard in the UK"

# ╔═╡ d790281e-f5dd-11ea-0d1c-f57da5018a6b
md"How much information do I need to store a multiplication table?"

# ╔═╡ d1578d4c-f601-11ea-2983-27dc131d39b8
md"### Scaled multiplication table"

# ╔═╡ d9556f32-f601-11ea-3dd8-1bc876b7b719
10 .* outer(1:10, 1:10)

# ╔═╡ e36e4ec2-f5dd-11ea-34ea-1bcf5fd7c16d
md"In the context of 1:n times 1:n, just one number n is needed.

But given arbitrary vectors v and w,  we need to store the whole of v and w"

# ╔═╡ 52e857ca-f5de-11ea-14bb-bdc0ac24ab90
md"Toeplitz in hw"

# ╔═╡ 98f08990-f5de-11ea-1f56-1f2d73649773
"Add mult tables to make flag"

# ╔═╡ 21bbb60a-f5df-11ea-2c1b-dd716a657df8
cs = distinguishable_colors(100)

# ╔═╡ a5d637ea-f5de-11ea-3b70-877e876bc9c9
flag = outer([1, 1, 1, 2, 2, 2, 1, 1, 1], ones(Int, 9))

# ╔═╡ 2668e100-f5df-11ea-12b0-073a578a5edb
cs[flag]

# ╔═╡ e8d727f2-f5de-11ea-1456-f72602e81e0d
cs[flag + flag']

# ╔═╡ 4c80c786-f5df-11ea-31ec-318439349648
cs[outer(rand(1:10, 3), rand(1:10, 5))]

# ╔═╡ 8d2bae22-f5df-11ea-10d3-859f4c3aa6c7
Gray.((100 .- outer(1:10, 1:10)) ./ 100)

# ╔═╡ e23debc8-f5df-11ea-2c1e-b58a64f9acd3
ColorSchemes.rainbow[floor.(Int, (1:10, 1:10) ./ 10)]

# ╔═╡ 70b8918e-f5e0-11ea-3c86-6fa72df5a28e
ColorSchemes.rainbow.colors[floor.(Int, outer(1:10, 1:10) ./ 10)]

# ╔═╡ a0934122-f5e0-11ea-1f3b-ab0021ac6906
ColorSchemes.rainbow.colors[ceil.(Int, outer(1:10, 1:10) ./ 10)]

# ╔═╡ b9381b8a-f5e0-11ea-1f84-e39325203038

# ╔═╡ 4cf96558-f5e0-11ea-19be-db4c59a41120
ColorSchemes.rainbow

# ╔═╡ 11de523c-f5e0-11ea-2f3d-c981c1b6a1fe
outer(1:10, 1:10) ./ 100

# ╔═╡ fb0c6c7e-f5df-11ea-38d0-2d98c9dc232f

# ╔═╡ ebd72fb8-f5e0-11ea-0630-573337dff753
md"## SVD"

# ╔═╡ f00d1eaa-f5e0-11ea-21df-d9cf6f7af9b9
md"Grab the rank-1 approx of an image.

First multiplication table"

# ╔═╡ b6478e1a-f5f6-11ea-3b92-6d4f067285f4
url = "https://arbordayblog.org/wp-content/uploads/2018/06/oak-tree-sunset-iStock-477164218.jpg"

# ╔═╡ d4a049a2-f5f8-11ea-2f34-4bc0e3a5954a
download(url, "tree.jpg")

# ╔═╡ f2c11f88-f5f8-11ea-3e02-c1d4fa22031e
begin
	image = load("tree.jpg")
	image = image[1:5:end, 1:5:end]
end

# ╔═╡ f7e38aaa-f5f8-11ea-002f-09dd1fa21181
reds = [Float64(c.r) for c in image]

# ╔═╡ 29062f7a-f5f9-11ea-2682-1374e7694e32
picture = Float64.(channelview(image))

# ╔═╡ 6156fd1e-f5f9-11ea-06a9-211c7ab813a4
pr, pg, pb = eachslice(picture, dims=1)

# ╔═╡ a9766e68-f5f9-11ea-0019-6f9d02050521
[RGB.(pr, 0, 0) RGB.(0, pg, 0) RGB.(0, 0, pb)]

# ╔═╡ fee66076-f5f9-11ea-2316-abc57b62a57c
RGB.(image

# ╔═╡ 6532b388-f5f9-11ea-2ae2-f9b12e441bb3
pr

# ╔═╡ 0c0ee362-f5f9-11ea-0f75-2d2810c88d65
begin
	Ur, Σr, Vr = svd(pr)
	Ug, Σg, Vg = svd(pg)
	Ub, Σb, Vb = svd(pb)
end

# ╔═╡ 3c28c4c2-f5fa-11ea-1947-9dfe91ea1535
RGB.(sum(outer(Ur[:,i], Vr[:,i]) .* Σr[i] for i in 1:5), 0, 0)

# ╔═╡ f56f40e4-f5fa-11ea-3a99-156565445c2e
@bind n Slider(1:100, show_value=true)

# ╔═╡ 7ba6e6a6-f5fa-11ea-2bcd-616d5a3c898b
RGB.(sum(outer(Ur[:,i], Vr[:,i]) .* Σr[i] for i in 1:n),
	sum(outer(Ug[:,i], Vg[:,i]) .* Σg[i] for i in 1:n),
	sum(outer(Ub[:,i], Vb[:,i]) .* Σb[i] for i in 1:n))

# ╔═╡ 8a22387e-f5fb-11ea-249b-435af5c0a6b6

# ╔═╡ 8df84fcc-f5d5-11ea-312f-bf2a3b3ce2ce
md"## Appendix"

# ╔═╡ 91980bcc-f5d5-11ea-211f-e9a08ff0fb19
function with_terminal(f)
	local spam_out, spam_err
	@color_output false begin
		spam_out = @capture_out begin
			spam_err = @capture_err begin
				f()
			end
		end
	end
	spam_out, spam_err

	HTML("""

$(Markdown.htmlesc(spam_out))

	""")
end

# ╔═╡ 466901ea-f5d5-11ea-1db5-abf82c96eabf
with_terminal() do
	dump(M)
end

# ╔═╡ b38c4aae-f5d5-11ea-39b6-7b0c7d529019
with_terminal() do
	dump(Diagonal(M))
end

# ╔═╡ 8b60629e-f5d6-11ea-27c8-d934460d3a57
with_terminal() do
	dump(M3)
end

# ╔═╡ cde79f38-f5d6-11ea-3297-0b5b240f7b9e
with_terminal() do
	dump(sparse(M4))
end

# ╔═╡ 4f8684ea-f5fb-11ea-07be-11d8046f35df
with_terminal() do
	@time @inbounds RGB.(sum(outer(Ur[:,i], Vr[:,i]) .* Σr[i] for i in 1:20),
	sum(outer(Ug[:,i], Vg[:,i]) .* Σg[i] for i in 1:20),
	sum(outer(Ub[:,i], Vb[:,i]) .* Σb[i] for i in 1:20));
end

# ╔═╡ Cell order:
# ╟─b0ba5b8c-f5d1-11ea-1304-3f0e47f935fe
# ╟─ffa95430-f5d1-11ea-3cb7-5bb8d8f13701
# ╟─261c4df2-f5d2-11ea-2c72-7d4b09c46098
# ╠═ae24c8b2-f60b-11ea-2c7a-03857d1217b2
# ╠═b5177f70-f60b-11ea-14a9-f5a574cc5185
# ╠═bc14fc1a-f60b-11ea-207a-91b967f28076
# ╠═ef8f44b2-f5fc-11ea-1e4d-bd873cd39d6c
# ╠═fd9211c0-f5fc-11ea-1745-7f2dae88af9e
# ╠═3a1ed5b8-f5fd-11ea-2ecd-b9d08349651f
# ╠═676f6e3a-f5fd-11ea-3552-e7e7c6382276
# ╠═77ae1146-f5d2-11ea-1226-27d15c90a8df
# ╠═3592648a-f5fd-11ea-2ed0-a1d4a4c14d07
# ╠═7ffa02f6-f5d2-11ea-1b41-15c21f92e9ee
# ╟─21a15906-f5d3-11ea-3a71-f53eabc31acc
# ╟─3c5f0ba8-f5d3-11ea-22d3-4f7c513144bc
# ╠═4d2db24e-f5d4-11ea-10d4-b126f643abee
# ╠═4624cd26-f5d3-11ea-1cf8-7d6555eae5fa
# ╠═397ac764-f5fe-11ea-20cc-8d7cab19d410
# ╠═82c7046c-f5d3-11ea-04e2-ef7c0f4db5da
# ╠═60214c1a-f5fe-11ea-2d08-c59715bcedd0
# ╠═7a98f708-f5d3-11ea-2301-cde71ca469ff
# ╠═cf77d83c-f5fd-11ea-136e-8951de64e05e
# ╠═5813e1b2-f5ff-11ea-2849-a1def74fc065
# ╠═e51d4ef6-f5fd-11ea-28f1-a1616e17f715
# ╠═982590d4-f5ff-11ea-3802-73292c75ad6c
# ╠═c3f70e18-f5ff-11ea-35aa-31b22ca506b8
# ╠═c200deb0-f5fd-11ea-32e1-d96596b16ebc
# ╠═a5989a76-f5d3-11ea-3fe5-5f9959199fe8
# ╠═b0bd6c6a-f5d3-11ea-00ee-b7c155f23481
# ╠═bd94f82c-f5d3-11ea-2a1e-77ddaadbfeb9
# ╠═8d2c6910-f5d4-11ea-1928-1baf09815687
# ╟─81c35324-f5d4-11ea-2338-9f982d38732c
# ╟─2cfda0dc-f5d5-11ea-16c4-b5a33b90e37f
# ╠═0feb5674-f5d5-11ea-0714-7379a7d381a3
# ╠═150432d4-f5d5-11ea-32b2-19a2a91d9637
# ╠═21328d1c-f5d5-11ea-288e-4171ad35326d
# ╠═466901ea-f5d5-11ea-1db5-abf82c96eabf
# ╠═b38c4aae-f5d5-11ea-39b6-7b0c7d529019
# ╟─d8f36278-f5d5-11ea-3338-8573ce40e65e
# ╟─e90c55fc-f5d5-11ea-10f1-470ff772985d
# ╠═19775c3c-f5d6-11ea-15c2-89618e654a1e
# ╠═20125236-f5d6-11ea-3877-6b332497be62
# ╠═232d1dcc-f5d6-11ea-2710-658c75b9c7a4
# ╠═2b2feb9c-f5d6-11ea-191c-df20360b12a9
# ╠═50c74f6c-f5d6-11ea-29f1-5de9997d5d9f
# ╠═5de72b7c-f5d6-11ea-1b6f-35b830b5fb34
# ╠═8b60629e-f5d6-11ea-27c8-d934460d3a57
# ╠═2fd7e52e-f5d7-11ea-3b5a-1f338e2451e0
# ╠═cde79f38-f5d6-11ea-3297-0b5b240f7b9e
# ╠═aa09c008-f5d8-11ea-1bdc-b51ee6eb2478
# ╟─d941cd66-f5d8-11ea-26ff-47ba7779ab20
# ╟─62a6ec62-f5d9-11ea-071e-ed33c5dea0cd
# ╠═67274c3c-f5d9-11ea-3475-c9d228e3bd5a
# ╟─b6c7a918-f600-11ea-18ff-6521507358c6
# ╟─765c6552-f5d9-11ea-29d3-bfe7b4b04612
# ╠═0f6ecba6-f5da-11ea-22c2-2929e562f413
# ╠═126fb3ea-f5da-11ea-2f7d-0b3259a296ce
# ╠═5f79e8f4-f5da-11ea-2b55-ef344b8a3ba2
# ╟─9b9e2c2a-f5da-11ea-369b-b513b196515b
# ╠═e68b98ea-f5da-11ea-1a9d-db45e4f80241
# ╠═f20ccac4-f5da-11ea-0e69-413b5e49f423
# ╠═0bc792e8-f5db-11ea-0b7a-1502ddc8008e
# ╠═12a2e96c-f5db-11ea-1c3e-494ae7446886
# ╠═22487ce2-f5db-11ea-32e9-6f70ab2c0353
# ╠═389ae62e-f5db-11ea-1557-c3adbbee0e5c
# ╟─0c2b6408-f5d9-11ea-2b7f-7fece2eecc1f
# ╟─542a9556-f5db-11ea-0375-99f52416f6e4
# ╠═165788b2-f601-11ea-3e69-cdbbb6558e54
# ╠═22941bb8-f601-11ea-1d6e-0d955297bc2e
# ╠═2f75df7e-f601-11ea-2fc2-aff4f335af33
# ╠═53e6b612-f601-11ea-05a9-5395e69b3c41
# ╠═5a493052-f601-11ea-2f5f-f940412905f2
# ╠═8633afb2-f601-11ea-206b-e9c4b9621c2a
# ╠═9918c4fa-f601-11ea-3bf1-3506dcb437f7
# ╟─3e919766-f601-11ea-0485-05f45484bf8d
# ╠═1052993e-f601-11ea-2c55-0d67e31b670e
# ╠═68190822-f5db-11ea-117f-d10a161208c3
# ╠═71b44874-f5db-11ea-1f67-47bad9295e03
# ╠═6c51eddc-f5db-11ea-1235-332cdbb072fa
# ╠═86fb49ee-f5db-11ea-3bfa-c95c3b8775a3
# ╠═173cfab4-f5d9-11ea-0c7c-bf8b0888f6e7
# ╠═8bab3e36-f5db-11ea-187a-f31fa8cf357d
# ╠═e64291e8-f5db-11ea-0cab-8567b781408f
# ╠═0111a124-f5dc-11ea-0904-fdd88d7acac4
# ╠═06406158-f5dc-11ea-02b1-2519a0176993
# ╠═12051f06-f5dc-11ea-2f6c-0fdc50eeff01
# ╠═195854e6-f5dc-11ea-114b-2333f87173f7
# ╠═768116f6-f5dc-11ea-1cfe-b3016c574725
# ╠═9050426e-f5dc-11ea-373c-65456732bd34
# ╠═965bd07e-f5dc-11ea-2e85-d34996cf2fae
# ╠═a9daf4cc-f5dc-11ea-270b-2566f89f168c
# ╠═dc0c8b72-f5dc-11ea-3e6f-0f43cbf58f56
# ╠═0ed78e76-f5dd-11ea-2ad8-a35a69c0ef9a
# ╠═1648a0fa-f5dd-11ea-0292-495207e83de9
# ╠═1e9cefea-f5dd-11ea-3e5f-a189fd41c42e
# ╠═24664f20-f5dd-11ea-2a69-cd3e0ebd5c39
# ╠═43a4920c-f5dd-11ea-1ab1-0b3d673c0f1e
# ╠═5d767290-f5dd-11ea-2189-81198fd216ce
# ╠═6aae805e-f5dd-11ea-108c-733daae313dc
# ╠═9a023cf8-f5dd-11ea-3016-f95d433e6df0
# ╠═b4c82246-f5dd-11ea-068f-2f63a5a382e2
# ╠═d1f87b22-f5dd-11ea-3bc3-471d5b3a5202
# ╠═d790281e-f5dd-11ea-0d1c-f57da5018a6b
# ╠═d1578d4c-f601-11ea-2983-27dc131d39b8
# ╠═d9556f32-f601-11ea-3dd8-1bc876b7b719
# ╠═e36e4ec2-f5dd-11ea-34ea-1bcf5fd7c16d
# ╠═52e857ca-f5de-11ea-14bb-bdc0ac24ab90
# ╠═98f08990-f5de-11ea-1f56-1f2d73649773
# ╠═21bbb60a-f5df-11ea-2c1b-dd716a657df8
# ╠═a5d637ea-f5de-11ea-3b70-877e876bc9c9
# ╠═2668e100-f5df-11ea-12b0-073a578a5edb
# ╠═e8d727f2-f5de-11ea-1456-f72602e81e0d
# ╠═4c80c786-f5df-11ea-31ec-318439349648
# ╠═8d2bae22-f5df-11ea-10d3-859f4c3aa6c7
# ╠═e23debc8-f5df-11ea-2c1e-b58a64f9acd3
# ╠═70b8918e-f5e0-11ea-3c86-6fa72df5a28e
# ╠═a0934122-f5e0-11ea-1f3b-ab0021ac6906
# ╠═b9381b8a-f5e0-11ea-1f84-e39325203038
# ╠═4cf96558-f5e0-11ea-19be-db4c59a41120
# ╠═11de523c-f5e0-11ea-2f3d-c981c1b6a1fe
# ╠═fb0c6c7e-f5df-11ea-38d0-2d98c9dc232f
# ╠═ebd72fb8-f5e0-11ea-0630-573337dff753
# ╠═f00d1eaa-f5e0-11ea-21df-d9cf6f7af9b9
# ╠═b6478e1a-f5f6-11ea-3b92-6d4f067285f4
# ╠═d4a049a2-f5f8-11ea-2f34-4bc0e3a5954a
# ╠═f2c11f88-f5f8-11ea-3e02-c1d4fa22031e
# ╠═f7e38aaa-f5f8-11ea-002f-09dd1fa21181
# ╠═29062f7a-f5f9-11ea-2682-1374e7694e32
# ╠═6156fd1e-f5f9-11ea-06a9-211c7ab813a4
# ╠═a9766e68-f5f9-11ea-0019-6f9d02050521
# ╠═fee66076-f5f9-11ea-2316-abc57b62a57c
# ╠═6532b388-f5f9-11ea-2ae2-f9b12e441bb3
# ╠═0c0ee362-f5f9-11ea-0f75-2d2810c88d65
# ╠═3c28c4c2-f5fa-11ea-1947-9dfe91ea1535
# ╠═f56f40e4-f5fa-11ea-3a99-156565445c2e
# ╠═7ba6e6a6-f5fa-11ea-2bcd-616d5a3c898b
# ╠═8a22387e-f5fb-11ea-249b-435af5c0a6b6
# ╠═4f8684ea-f5fb-11ea-07be-11d8046f35df
# ╠═8df84fcc-f5d5-11ea-312f-bf2a3b3ce2ce
# ╠═91980bcc-f5d5-11ea-211f-e9a08ff0fb19

max_stars_repo_path: lecture_notebooks/taking advantage of structure.jl

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
abstract type Community_Agent <: Agent end

@with_kw mutable struct DER_Agent <: Community_Agent
    name::String = "anonymus"
    location::String = "under_AGG_1"
    period::Period = Period()
    demand::Any = Demand(period = period) #Type Any was the only thing working on passing in multiple elements, to be reslved
    dev_ID::Any = Demand(period = period)
    dev_RT::Any = Demand(period = period)
    preferences::Preferences = Preferences()
    battery::Any = Battery()
    pv = PV()
    OPEX_DA::Float64 = 0.0
    OPEX_ID::Float64 = 0.0
    OPEX_RT::Float64 = 0.0
    #HERE A FUNCTION SHOULD BE ADDED TO CHECK THE SUM OF THE PROBABILITIES PER FORECAST
end


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataFrames
using CSV

function products_report(solution; marginal_costs = true)::DataFrame
    df = DataFrame()
    df."product name" = String[]
    df."location name" = String[]
    df."latitude (deg)" = Float64[]
    df."longitude (deg)" = Float64[]
    df."year" = Int[]
    df."amount (tonne)" = Float64[]
    df."marginal cost (\$/tonne)" = Float64[]
    T = length(solution["Energy"]["Plants (GJ)"])
    for (prod_name, prod_dict) in solution["Products"]
        for (location_name, location_dict) in prod_dict
            for year = 1:T
                marginal_cost = location_dict["Marginal cost (\$/tonne)"][year]
                latitude = round(location_dict["Latitude (deg)"], digits = 6)
                longitude = round(location_dict["Longitude (deg)"], digits = 6)
                amount = location_dict["Amount (tonne)"][year]
                push!(
                    df,
                    [
                        prod_name,
                        location_name,
                        latitude,
                        longitude,
                        year,
                        amount,
                        marginal_cost,
                    ],
                )
            end
        end
    end
    return df
end

write_products_report(solution, filename) = CSV.write(filename, products_report(solution))

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
# https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html

import Base: +, -, *, /

"""
Quaternion

``q = w + x⋅i + y⋅j + z⋅k``
"""
@export struct Qt <: Number
   w::Cfloat
   x::Cfloat
   y::Cfloat
   z::Cfloat
   Qt(w::Real = 0, x::Real = 0, y::Real = 0, z::Real = 0) = new(w, x, y, z)
end

Qt(u::Vec{3}, θ::Real) = Qt(cos(θ/2), (u * sin(θ/2))...)
Qt(v::Vec{3}) = Qt(0, v...)
Qt(v::Vec{4}) = Qt(v...)
Qt(v::Vector{<:Real}) = Qt(Vec(v...))

const i = Qt(0, 1)
const j = Qt(0, 0, 1)
const k = Qt(0, 0, 0, 1)

Base.getindex(q::Qt, i::Integer) = getfield(q, i)
Base.iterate(q::Qt, i = 1) = i > 4 ? nothing : (q[i], i + 1)

Vec{3}(q::Qt) = Vec(q.x, q.y, q.z)
Vec(q::Qt) = Vec(q...)

Base.promote_rule(::Type{Qt}, ::Type{<:Real}) = Qt
Base.convert(::Type{Qt}, n::Real) = Qt(n)

-(q::Qt) = Qt(-1 * Vec(q))
-(q::Qt, p::Qt) = Qt(Vec(q) - Vec(p))
+(q::Qt, p::Qt) = Qt(Vec(q) + Vec(p))

*(q::Qt, n::Real) = Qt(n * Vec(q))
*(n::Real, q::Qt) = Qt(n * Vec(q))

/(q::Qt, n::Real) = Qt(Vec(q) / n)
/(::Real, ::Qt)   = error("Cannot divide real numbers by Quaternions")
/(q::Qt, p::Qt) = q * inv(p)

LinearAlgebra.norm(q::Qt) = norm(Vec(q))

Base.conj(q::Qt) = Qt(q.w) - Qt(Vec{3}(q))
Base.adjoint(q::Qt) = conj(q)
Base.inv(q::Qt) = q' / norm(q)

Mat{4}((a, b, c, d)::Qt) = Mat4([
   a  -b  -c  -d
   b   a  -d   c
   c   d   a  -b
   d  -c   b   a
])

*((w₁, x₁, y₁, z₁)::Qt, (w₂, x₂, y₂, z₂)::Qt) = Qt(
   (w₁ * w₂ - x₁ * x₂ - y₁ * y₂ - z₁ * z₂),
   (w₁ * x₂ + x₁ * w₂ + y₁ * z₂ - z₁ * y₂),
   (w₁ * y₂ + y₁ * w₂ + z₁ * x₂ - x₁ * z₂),
   (w₁ * z₂ + z₁ * w₂ + x₁ * y₂ - y₁ * x₂),
)

@export begin
   const X = Vec(1, 0, 0)
   const Y = Vec(0, 1, 0)
   const Z = Vec(0, 0, 1)
end

"""
      rotation(q::Qt) -> Mat4

Creates a `4x4` rotation matrix from the quaternion `q`.
"""
@export rotation((w, x, y, z)::Qt) = Mat4([
   (1 - 2y^2 - 2z^2) (2x*y - 2w*z)     (2x*z + 2w*y)      0
   (2x*y + 2w*z)     (1 - 2x^2 - 2z^2) (2y*z - 2w*x)      0
   (2x*z - 2w*y)     (2y*z + 2w*x)     (1 - 2x^2 - 2y^2)  0
   0                 0                 0                  1
])

function *(q::Qt, v::Vec{N}) where {N}
   @assert 4 >= N >= 1 "
      Quaternions can only act as rotations on 1-4 dimensional vectors.
      Got vector with $N dimensions
   "
   return rotation(q)[1:N, 1:N] * v
end

function Base.show(io::IO, (w, x, y, z)::Qt)
   f(x) = abs(x) == 1 ? "- $(abs(x))" : "+ $(x)"
   write(io, "($w $(f(x))i $(f(y))j $(f(z))k)")
end

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays

include("memorylayouttests.jl")
include("multests.jl")
include("ldivtests.jl")

@testset "concat" begin
    @testset "Vcat" begin
        A = Vcat(Vector(1:10), Vector(1:20))
        @test @inferred(length(A)) == 30
        @test @inferred(A[5]) == A[15] == 5
        @test_throws BoundsError A[31]
        @test reverse(A) == Vcat(Vector(reverse(1:20)), Vector(reverse(1:10)))
        b = Array{Int}(undef, 31)
        @test_throws DimensionMismatch copyto!(b, A)
        b = Array{Int}(undef, 30)
        @test @allocated(copyto!(b, A)) == 0
        @test b == vcat(A.arrays...)

        A = Vcat(1:10, 1:20)
        @test @inferred(length(A)) == 30
        @test @inferred(A[5]) == A[15] == 5
        @test_throws BoundsError A[31]
        @test reverse(A) == Vcat(reverse(1:20), reverse(1:10))
        b = Array{Int}(undef, 31)
        @test_throws DimensionMismatch copyto!(b, A)
        b = Array{Int}(undef, 30)
        @test @allocated(copyto!(b, A)) == 0
        @test b == vcat(A.arrays...)

        A = Vcat(randn(2,10), randn(4,10))
        @test @inferred(length(A)) == 60
        @test @inferred(size(A)) == (6,10)
        @test_throws BoundsError A[61]
        @test_throws BoundsError A[7,1]
        b = Array{Float64}(undef, 7,10)
        @test_throws DimensionMismatch copyto!(b, A)
        b = Array{Float64}(undef, 6,10)
        @test @allocated(copyto!(b, A)) == 0
        @test b == vcat(A.arrays...)

        @test_throws ArgumentError Vcat()
    end
    @testset "Hcat" begin
        A = Hcat(1:10, 2:11)
        @test_throws BoundsError A[1,3]
        @test @inferred(size(A)) == (10,2)
        @test @inferred(A[5]) == @inferred(A[5,1]) == 5
        @test @inferred(A[11]) == @inferred(A[1,2]) == 2
        b = Array{Int}(undef, 11, 2)
        @test_throws DimensionMismatch copyto!(b, A)
        b = Array{Int}(undef, 10, 2)
        @test @allocated(copyto!(b, A)) == 0
        @test b == hcat(A.arrays...)

        A = Hcat(Vector(1:10), Vector(2:11))
        b = Array{Int}(undef, 10, 2)
        copyto!(b, A)
        @test b == hcat(A.arrays...)
        @test @allocated(copyto!(b, A)) == 0

        A = Hcat(1, zeros(1,5))
        @test A == hcat(1, zeros(1,5))

        A = Hcat(Vector(1:10), randn(10, 2))
        b = Array{Float64}(undef, 10, 3)
        copyto!(b, A)
        @test b == hcat(A.arrays...)
        @test @allocated(copyto!(b, A)) == 0
    end

    @testset "Special pads" begin
        A = Vcat([1,2,3], Zeros(7))
        B = Vcat([1,2], Zeros(8))

        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.arrays[1] isa Vector{Float64}
        @test C.arrays[2] isa Zeros{Float64}
        @test C == Vector(A) + Vector(B)

        B = Vcat([1,2], Ones(8))

        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.arrays[1] isa Vector{Float64}
        @test C.arrays[2] isa Ones{Float64}
        @test C == Vector(A) + Vector(B)

        B = Vcat([1,2], randn(8))

        C = @inferred(A+B)
        @test C isa BroadcastArray{Float64}
        @test C == Vector(A) + Vector(B)

        B = Vcat(SVector(1,2), Ones(8))
        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.arrays[1] isa Vector{Float64}
        @test C.arrays[2] isa Ones{Float64}
        @test C == Vector(A) + Vector(B)

        A = Vcat(SVector(3,4), Zeros(8))
        B = Vcat(SVector(1,2), Ones(8))
        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.arrays[1] isa SVector{2,Int}
        @test C.arrays[2] isa Ones{Float64}
        @test C == Vector(A) + Vector(B)
    end
end

@testset "Kron"  begin
    A = [1,2,3]
    B = [4,5,6,7]

    @test Array(Kron(A)) == A
    K = Kron(A,B)
    @test [K[k] for k=1:length(K)] == Array(K) == kron(A,B)

    A = randn(3)
    K = Kron(A,B)
    @test K isa Kron{Float64}
    @test all(isa.(K.arrays, Vector{Float64}))
    @test [K[k] for k=1:length(K)] == Array(K) == Array(Kron{Float64}(A,B)) == kron(A,B)

    # C = [7,8,9,10,11]
    # K = Kron(A,B,C)
    # @time [K[k] for k=1:length(K)] == Array(Kron(A,B)) == kron(A,B)

    A = randn(3,2)
    B = randn(4,6)
    K = Kron(A,B)
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == kron(A,B)
    K = Kron(A,B')
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B')) == kron(A,B')
    K = Kron(A',B)
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B)) == kron(A',B)
    K = Kron(A',B')
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B')) == kron(A',B')

    A = randn(3,2)
    B = randn(4,6)
    K = Kron(A,B)
    C = similar(K)
    copyto!(C, K)
    @test C == kron(A,B)
    @test_throws DimensionMismatch randn(2,2) .= K

    A = rand(Int,3,2)
    K = Kron(A,B)
    @test K isa Kron{Float64}
    @test all(isa.(K.arrays, Matrix{Float64}))
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(K) == Array(Kron{Float64}(A,B)) == kron(A,B)
end

@testset "BroadcastArray" begin
    A = randn(6,6)
    B = BroadcastArray(exp, A)
    @test Matrix(B) == exp.(A)

    C = BroadcastArray(+, A, 2)
    @test C == A .+ 2
    D = BroadcastArray(+, A, C)
    @test D == A + C

    x = Vcat([3,4], [1,1,1,1,1], 1:3)
    @test x .+ (1:10) isa Vcat
    @test (1:10) .+ x isa Vcat
    @test x + (1:10) isa Vcat
    @test (1:10) + x isa Vcat
    @test x .+ (1:10) == (1:10) .+ x == (1:10) + x == x + (1:10) == Vector(x) + (1:10)

    @test exp.(x) isa Vcat
    @test exp.(x) == exp.(Vector(x))
    @test x .+ 2 isa Vcat
    @test (x .+ 2).arrays[end] ≡ x.arrays[end] .+ 2 ≡ 3:5
    @test x .* 2 isa Vcat
    @test 2 .+ x isa Vcat
    @test 2 .* x isa Vcat
end

@testset "Cache" begin
    A = 1:10
    C = cache(A)
    @test size(C) == (10,)
    @test axes(C) == (Base.OneTo(10),)
    @test all(Vector(C) .=== Vector(A))

    A = reshape(1:10^2, 10,10)
    C = cache(A)
    @test size(C) == (10,10)
    @test axes(C) == (Base.OneTo(10),Base.OneTo(10))
    @test all(Array(C) .=== Array(A))

    A = reshape(1:10^3, 10,10,10)
    C = cache(A)
    @test size(C) == (10,10,10)
    @test axes(C) == (Base.OneTo(10),Base.OneTo(10),Base.OneTo(10))
    @test all(Array(C) .=== Array(A))

    A = reshape(1:10^3, 10,10,10)
    C = cache(A)
    LazyArrays.resizedata!(C,5,5,5)
    LazyArrays.resizedata!(C,8,8,8)
    @test all(C.data .=== Array(A)[1:8,1:8,1:8])
end

@testset "Diff and Cumsum" begin
    x = Vcat([3,4], [1,1,1,1,1], 3)
    y = @inferred(cumsum(x))
    @test y isa Vcat
    @test y == cumsum(Vector(x)) == Cumsum(x)
    @test @inferred(diff(x)) == diff(Vector(x)) == Diff(x)

    @test sum(x) == sum(Vector(x)) == last(y)
    @test cumsum(Vcat(4)) === Vcat(4)

    A = randn(3,4)
    @test Cumsum(A; dims=1) == cumsum(A; dims=1)
    @test Cumsum(A; dims=2) == cumsum(A; dims=2)
    @test Diff(A; dims=1) == diff(A; dims=1)
    @test Diff(A; dims=2) == diff(A; dims=2)
end

@testset "broadcast Vcat" begin
    x = Vcat(1:2, [1,1,1,1,1], 3)
    y = 1:8
    f = (x,y) -> cos(x*y)
    @test f.(x,y) isa Vcat
    @test @inferred(broadcast(f,x,y)) == f.(Vector(x), Vector(y))

    @test (x .+ y) isa Vcat
    @test (x .+ y).arrays[1] isa AbstractRange
    @test (x .+ y).arrays[end] isa Int

    z = Vcat(1:2, [1,1,1,1,1], 3)
    (x .+ z)  isa BroadcastArray
    (x + z) isa BroadcastArray
    @test Vector( x .+ z) == Vector( x + z) == Vector(x) + Vector(z)

    # Lazy mixed with Static treats as Lazy
    s = SVector(1,2,3,4,5,6,7,8)
    @test f.(x , s) isa Vcat
    @test f.(x , s) == f.(Vector(x), Vector(s))

    # these are special cased
    @test Vcat(1, Ones(5))  + Vcat(2, Fill(2.0,5)) ≡ Vcat(3, Fill(3.0,5))
    @test Vcat(SVector(1,2,3), Ones(5))  + Vcat(SVector(4,5,6), Fill(2.0,5)) ≡
        Vcat(SVector(5,7,9), Fill(3.0,5))
end

@testset "maximum/minimum Vcat" begin
    x = Vcat(1:2, [1,1,1,1,1], 3)
    @test maximum(x) == 3
    @test minimum(x) == 1
end


text:
# Verify that pulling data directly out of the SparseMatrixCSC form gives
# exactly what the previous hand-coded algorithm gave.
@testset "sparse matrix" begin
  function sparse_matrix_reference(A::AbstractSparseMatrix)
    m, n = size(A)
    @assert m==n

    col_start = Array{Int}(undef, n)
    col_len = Array{Int}(undef, n)
    row = Array{Int}(undef, 0)
    data = Array{Float64}(undef, 0)
    for j in 1:n
      if j==1
        col_start[j] = 1
      else
        col_start[j] = col_start[j-1] + col_len[j-1]
      end

      col_len[j] = 0
      for i in 1:n
        if A[i,j] != 0.0
          col_len[j] += 1
          push!(row, i)
          push!(data, A[i,j])
        end
      end
    end
    return col_start, col_len, row, data
  end

  sparse_matrix_reference(A::Matrix) = sparse_matrix_reference(sparse(A))

  function sparse_matrix_csc(A::SparseMatrixCSC)
    m, n = size(A)
    @assert m==n
    col_start = A.colptr[1:end-1]
    col_len = diff(A.colptr)
    row = rowvals(A)
    data = nonzeros(A)
    return col_start, col_len, row, data
  end

  sparse_matrix_csc(A::Matrix) = sparse_matrix_csc(convert(SparseMatrixCSC, A))

  Random.seed!(42)
  for i in 1:100
    n = rand(5:20)
    M = sprandn(n, n, 0.1)
    @test sparse_matrix_csc(M) == sparse_matrix_reference(M)

    Mf = Matrix(M)
    @test sparse_matrix_csc(Mf) == sparse_matrix_reference(Mf)
    @test sparse_matrix_csc(Mf) == sparse_matrix_csc(M)
  end
end


text:
function _clean_high_level_definition(definition::String)
    # Required Julia 1.5 or higher with how triple quoted strings are dealt with.
    definition = replace(definition, " " => "")
    definition = replace(definition, "\n" => "")

    return definition
end

@testset "_get_aws_sdk_js_files" begin
    apply(Patches._github_tree_patch) do
        files = _get_aws_sdk_js_files("foobar", GitHub.OAuth2("foobar"))

        @test length(files) == 1
        @test files[1] == OrderedDict("path" => "test-2020-01-01.normal.json")
    end
end

@testset "_filter_latest_service_version" begin
    migration_hub_v1 = Dict("path" => "AWSMigrationHub-2017-05-31.normal.json")
    migration_hub_v2 = Dict("path" => "AWSMigrationHub-2020-01-01.normal.json")

    access_analyzer_v1 = Dict("path" => "accessanalyzer-2019-11-01.normal.json")
    access_analyzer_v2 = Dict("path" => "accessanalyzer-2020-01-01.normal.json")

    @testset "empty" begin
        @test isempty(_filter_latest_service_version([]))
    end

    @testset "single service - single version" begin
        result = _filter_latest_service_version([migration_hub_v1])
        @test result == [migration_hub_v1]
    end

    @testset "single service - multiple versions" begin
        result = _filter_latest_service_version([migration_hub_v1, migration_hub_v2])
        @test result == [migration_hub_v2]
    end

    @testset "multiple service - single version" begin
        result = _filter_latest_service_version([migration_hub_v1, access_analyzer_v1])
        @test result == [access_analyzer_v1, migration_hub_v1]
    end

    @testset "multiple services - multiple versions" begin
        result = _filter_latest_service_version([migration_hub_v1, migration_hub_v2, access_analyzer_v1, access_analyzer_v2])
        @test result == [access_analyzer_v2, migration_hub_v2]
    end
end

@testset "_get_service_and_version" begin
    @testset "empty string" begin
        filename = ""
        @test_throws InvalidFileName _get_service_and_version(filename)
    end

    @testset "invalid filename" begin
        filename = "This is an invalid file name."
        @test_throws InvalidFileName _get_service_and_version(filename)
    end

    @testset "valid filename" begin
        filename = "AWSMigrationHub-2017-05-31.normal.json"
        service, version = _get_service_and_version(filename)

        @test service == "AWSMigrationHub"
        @test version == "2017-05-31"
    end
end

@testset "_generate_low_level_definitions" begin
    services = JSON.parsefile(joinpath(@__DIR__, "resources/services.json"))

    @testset "rest-xml" begin
        expected = "const s3 = AWS.RestXMLService(\"s3\", \"2006-03-01\")"
        response = _generate_low_level_definition(services["s3"])

        @test response == expected
    end

    @testset "rest-json" begin
        expected = "const glacier = AWS.RestJSONService(\"glacier\", \"2012-06-01\", LittleDict(\"x-amz-glacier-version\" => \"2012-06-01\"))"
        response = _generate_low_level_definition(services["glacier"])

        @test response == expected
    end

    @testset "ec2 / query" begin
        expected = "const ec2 = AWS.QueryService(\"ec2\", \"2016-11-15\")"
        response = _generate_low_level_definition(services["ec2"])

        @test response == expected
    end

    @testset "json" begin
        expected = "const budgets = AWS.JSONService(\"budgets\", \"2016-10-20\", \"1.1\", \"AWSBudgetServiceGateway\")"
        response = _generate_low_level_definition(services["budgets"])

        @test response == expected
    end

    @testset "signingName preference" begin
        expected = "const serviceid = AWS.RestXMLService(\"signingName\", \"2021-04-09\")"
        response = _generate_low_level_definition(services["signingName"])

        @test response == expected
    end

    @testset "invalid protocol" begin
        @test_throws ProtocolNotDefined _generate_low_level_definition(services["invalid"])
    end
end

@testset "_generate_low_level_definition" begin
    service = Dict(
        "serviceId" => "sample_service",
        "protocol" => "invalid-protocol",
        "endpointPrefix" => "endpoint",
        "apiVersion" => "api_version",
        "jsonVersion" => "json_version",
        "targetPrefix" => "target"
    )

    @testset "Invalid Protocol" begin
        @test_throws ProtocolNotDefined _generate_low_level_definition(service)
    end

    @testset "rest-xml" begin
        service["protocol"] = "rest-xml"
        expected_result = "const sample_service = AWS.RestXMLService(\"endpoint\", \"api_version\")"
        result = _generate_low_level_definition(service)

        @test result == expected_result
    end

    @testset "rest-json" begin
        service["protocol"] = "rest-json"
        expected_result = "const sample_service = AWS.RestJSONService(\"endpoint\", \"api_version\")"
        result = _generate_low_level_definition(service)

        @test result == expected_result
    end

    @testset "json" begin
        service["protocol"] = "json"
        expected_result = "const sample_service = AWS.JSONService(\"endpoint\", \"api_version\", \"json_version\", \"target\")"
        result = _generate_low_level_definition(service)

        @test result == expected_result
    end

    @testset "query" begin
        service["protocol"] = "query"
        expected_result = "const sample_service = AWS.QueryService(\"endpoint\", \"api_version\")"
        result = _generate_low_level_definition(service)

        @test result == expected_result
    end

    @testset "ec2" begin
        service["protocol"] = "ec2"
        expected_result = "const sample_service = AWS.QueryService(\"endpoint\", \"api_version\")"
        result = _generate_low_level_definition(service)

        @test result == expected_result
    end
end

@testset "_clean_documentation" begin
    documentation = "
To remove one or more tags, use the RemoveTagsFromCertificate action. \$ \\ To view all of the tags that have been applied to the certificate, use the ListTagsForCertificate action."
    expected_result = "To remove one or more tags, use the RemoveTagsFromCertificate action.   To view all of the tags that have been applied to the certificate, use the ListTagsForCertificate action."
    result = _clean_documentation(documentation)

    @test result == expected_result
end

@testset "_clean_uri" begin
    @testset "no parameters" begin
        uri = "/v1/configurations/"
        expected = "/v1/configurations/"

        @test _clean_uri(uri) == expected
    end

    @testset "single parameter" begin
        uri = "/v1/configurations/{parameter-one}"
        expected = "/v1/configurations/\$(parameter_one)"

        @test _clean_uri(uri) == expected
    end

    @testset "multiple parameters" begin
        uri = "/v1/configurations/{parameter_one}/{parameter_two}"
        expected = "/v1/configurations/\$(parameter_one)/\$(parameter_two)"

        @test _clean_uri(uri) == expected
    end

    @testset "hyphen not in parameter" begin
        uri = "/v1/configuration-parameters/{parameter-one}"
        expected = "/v1/configuration-parameters/\$(parameter_one)"

        @test _clean_uri(uri) == expected
    end

    @testset "remove plus signs" begin
        uri = "/v1/configuration-parameters/{parameter-one+}"
        expected = "/v1/configuration-parameters/\$(parameter_one)"

        @test _clean_uri(uri) == expected
    end
end

@testset "_format_name" begin
    @testset "single captial" begin
        function_name = "Testfunctionname"
        expected = "testfunctionname"

        @test _format_name(function_name) == expected
    end

    @testset "all capitals" begin
        function_name = "TESTFUNCTIONNAME"
        expected = "testfunctionname"

        @test _format_name(function_name) == expected
    end

    @testset "subsequent capitals" begin
        function_name = "TestFUNCTIONName"
        expected = "test_functionname"

        @test _format_name(function_name) == expected
    end

    @testset "normal" begin
        function_name = "TestFunctionName"
        expected = "test_function_name"

        @test _format_name(function_name) == expected
    end
end

@testset "_get_function_parameters" begin
    shapes = JSON.parsefile(joinpath(@__DIR__, "resources/shapes.json"))

    @testset "required params" begin
        input = "RequiredParams"

        required_params, optional_params = _get_function_parameters(input, shapes)

        @test required_params == Dict("RequiredParam" => LittleDict("location" => "", "documentation" => "Required param"))
        @test isempty(optional_params)
    end

    @testset "optional params" begin
        input = "OptionalParams"

        required_params, optional_params = _get_function_parameters(input, shapes)

        @test isempty(required_params)
        @test optional_params == Dict(
            "OptionalParam1" => Dict("documentation" => "Optional param 1", "idempotent" => false),
            "OptionalParam2" => Dict("documentation" => "Optional param 2", "idempotent" => false)
        )
    end

    @testset "required and optional params" begin
        input = "RequiredAndOptionalParams"

        required_params, optional_params = _get_function_parameters(input, shapes)

        @test required_params == Dict(
            "RequiredParam1" => LittleDict("location" => "", "documentation" => "Required param 1"),
            "RequiredParam2" => LittleDict("location" => "", "documentation" => "Required param 2")
        )

        @test optional_params == Dict(
            "OptionalParam" => Dict("documentation" => "Optional param", "idempotent" => false)
        )
    end

    @testset "no params" begin
        input = "NoParams"

        required_params, optional_params = _get_function_parameters(input, shapes)

        @test isempty(required_params)
        @test isempty(optional_params)
    end
end

@testset "_generate_high_level_definitions" begin
    service_name = "sample_service"
    protocol = "rest-xml"
    operations = JSON.parsefile(joinpath(@__DIR__, "resources/operations.json"))
    shapes = JSON.parsefile(joinpath(@__DIR__, "resources/shapes.json"))

    expected_result = """
    \"\"\"
        sample_operation(required_param1, required_param2)
        sample_operation(required_param1, required_param2, params::Dict{String,<:Any})

    The documentation for this operation.

    # Arguments
    - `required_param1`: Required param 1
    - `required_param2`: Required param 2

    # Optional Parameters
    Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
    - `"OptionalParam"`: Optional param
    \"\"\"
    sample_operation(RequiredParam1, RequiredParam2; aws_config::AbstractAWSConfig=global_aws_config()) = sample_service("POST", "/", Dict{String, Any}("RequiredParam1"=>RequiredParam1, "RequiredParam2"=>RequiredParam2); aws_config=aws_config)
    sample_operation(RequiredParam1, RequiredParam2, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sample_service("POST", "/", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RequiredParam1"=>RequiredParam1, "RequiredParam2"=>RequiredParam2), params)); aws_config=aws_config)
    """

    result = _generate_high_level_definitions(
        service_name,
        protocol,
        operations,
        shapes
    )

    @test size(result)[1] == 1

    expected_result = _clean_high_level_definition(expected_result)
    result = _clean_high_level_definition(result[1])

    @test result == expected_result
end

@testset "_generate_high_level_definition" begin
    service_name = "service_name"
    name = "FunctionName"
    method = "GET"
    request_uri = "request_uri"
    documentation = "Documentation for $name."

    @testset "locationless and non-idempotent" begin
        required_params = Dict("RequiredParam" => Dict("location" => "", "documentation" => "This parameter is required."))
        optional_params = Dict("OptionalParam" => Dict("idempotent" => false, "documentation" => "This parameter is optional."))

        @testset "rest protocol" begin
            protocol = "rest-xml"
            expected_result = """
            \"\"\"
                function_name(required_param)
                function_name(required_param, params::Dict{String,<:Any})

            Documentation for FunctionName.

            # Arguments
            - `required_param`: This parameter is required.

            # Optional Parameters
            Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
            - `"OptionalParam"`: This parameter is optional.
            \"\"\"
            function_name(RequiredParam; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("GET", "request_uri", Dict{String, Any}("RequiredParam"=>RequiredParam); aws_config=aws_config)
            function_name(RequiredParam, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("GET", "request_uri", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RequiredParam"=>RequiredParam), params)); aws_config=aws_config)
            """

            result = _generate_high_level_definition(
                service_name,
                protocol,
                name,
                method,
                request_uri,
                required_params,
                optional_params,
                documentation
            )

            expected_result = _clean_high_level_definition(expected_result)
            result = _clean_high_level_definition(result)

            @test result == expected_result
        end

        @testset "json, query, ec2 protocol" begin
            protocol = "ec2"
            expected_result = """
            \"\"\"
                function_name(required_param)
                function_name(required_param, params::Dict{String,<:Any})

            Documentation for FunctionName.

            # Arguments
            - `required_param`: This parameter is required.

            # Optional Parameters
            Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
            - `"OptionalParam"`: This parameter is optional.
            \"\"\"
            function_name(RequiredParam; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("FunctionName", Dict{String, Any}("RequiredParam"=>RequiredParam); aws_config=aws_config)
            function_name(RequiredParam, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("FunctionName", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RequiredParam"=>RequiredParam), params)); aws_config=aws_config)
            """

            result = _generate_high_level_definition(
                service_name,
                protocol,
                name,
                method,
                request_uri,
                required_params,
                optional_params,
                documentation
            )

            expected_result = _clean_high_level_definition(expected_result)
            result = _clean_high_level_definition(result)

            @test result == expected_result
        end
    end

    @testset "header location and idempotent" begin
        required_params = Dict("RequiredParam" => Dict("location" => "header", "documentation" => "This parameter is required."))
        optional_params = Dict("OptionalParam" => Dict("idempotent" => true, "documentation" => "This parameter is optional."))

        @testset "rest protocol" begin
            protocol = "rest-xml"
            expected_result = """
            \"\"\"
                function_name(required_param)
                function_name(required_param, params::Dict{String,<:Any})

            Documentation for FunctionName.

            # Arguments
            - `required_param`: This parameter   is required.

            # Optional Parameters
            Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
            - `"OptionalParam"`: This parameter i  s optional.
            \"\"\"
            function_name(RequiredParam; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("GET", "request_uri", Dict{String, Any}("OptionalParam"=>string(uuid4()), "headers"=>Dict{String, Any}("RequiredParam"=>RequiredParam)); aws_config=aws_config)
            function_name(RequiredParam, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("GET", "request_uri", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("OptionalParam"=>string(uuid4()), "headers"=>Dict{String, Any}("RequiredParam"=>RequiredParam)), params)); aws_config=aws_config)
            """
            result = _generate_high_level_definition(
                service_name,
                protocol,
                name,
                method,
                request_uri,
                required_params,
                optional_params,
                documentation
            )

            expected_result = _clean_high_level_definition(expected_result)
            result = _clean_high_level_definition(result)

            @test result == expected_result
        end

        @testset "json, query, ec2 protocol" begin
            protocol = "ec2"
            expected_result = """
            \"\"\"
                function_name(required_param)
                function_name(required_param, params::Dict{String,<:Any})

            Documentation for FunctionName.

            # Arguments
            - `required_param`: This parameter   is required.

            # Optional Parameters
            Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
            - `"OptionalParam"`: This parameter i  s optional.
            \"\"\"
            function_name(RequiredParam; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("FunctionName", Dict{String, Any}("RequiredParam"=>RequiredParam, "OptionalParam"=>string(uuid4())); aws_config=aws_config)
            function_name(RequiredParam, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_name("FunctionName", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RequiredParam"=>RequiredParam, "OptionalParam"=>string(uuid4())), params)); aws_config=aws_config)
            """

            result = _generate_high_level_definition(
                service_name,
                protocol,
                name,
                method,
                request_uri,
                required_params,
                optional_params,
                documentation
            )

            expected_result = _clean_high_level_definition(expected_result)
            result = _clean_high_level_definition(result)

            @test result == expected_result
        end
    end
end

@testset "string wrapping functionality" begin
    @testset "_validindex" begin
        str = "jμΛIα"  # μ, Λ, α have 2 codeunits.
        @test _validindex(str, 1) == 1 # j
        @test _validindex(str, 2) == 2 # μ
        @test _validindex(str, 3) == 2 # μ still
        @test _validindex(str, 4) == 4 # Λ
        @test_throws BoundsError _validindex(str, 0)
        @test_throws BoundsError _validindex(str, 10)
    end

    @testset "_splitline" begin
        str = "This is a short sentence."

        @testset "limit < 1" begin
            @test_throws DomainError _splitline(str, 0)
            @test_throws DomainError _splitline(str, -1)
        end

        @testset "limit == 1" begin
            result = _splitline(str, 1)
            @test result isa Tuple{String,String}
            line1, line2 = result
            @test line1 == string(first(str)) == "T"
            @test line2 == str[2:end] == "his is a short sentence."
        end

        @testset "limit >= ncodeunits" begin
            for limit in (ncodeunits(str), ncodeunits(str) + 1)
                result = _splitline(str, limit)
                @test result isa Tuple{String,String}
                line1, line2 = result
                @test line1 == str
                @test line2 == ""
            end
        end

        @testset "split on whitespace when possible" begin
            abc = "Aa Bb Cc"
            @test _splitline(abc, 1) == ("A", "a Bb Cc")  # No preceding whitespace to split on
            @test _splitline(abc, 2) == ("Aa", " Bb Cc")
            @test _splitline(abc, 3) == ("Aa ", "Bb Cc")
            @test _splitline(abc, 4) == ("Aa ", "Bb Cc")  # 4 == `B`, split on preceding whitespace
            @test _splitline(abc, 5) == ("Aa ", "Bb Cc")  # 5 == 'b', split on preceding whitespace
            @test _splitline(abc, 6) == ("Aa Bb ", "Cc")
            @test _splitline(abc, ncodeunits(abc) - 1) == ("Aa Bb ", "Cc")
        end

        @testset "does not try to split mid-character" begin
            str = "jμΛIα"  # 'μ' starts at str[2], 'Λ' starts at str[4]
            @test _splitline(str, 2) == ("jμ", "ΛIα")
            @test _splitline(str, 3) == ("jμ", "ΛIα") # should not try to split mid-'μ'
            @test _splitline(str, 4) == ("jμΛ", "Iα")
        end

        @testset "does not split on punctuation" begin
            str = "\"arn:aws:health:us-west-1::event/EBS/AWS\""
            result = _splitline(str, ncodeunits(str) - 1)
            # don't split escaped closing quote `\"` into `\` and `"`
            @test result == ("\"arn:aws:health:us-west-1::event/EBS/AWS", "\"")
        end
    end

    @testset "_wraplines" begin
        str = "This sentence contains exactly `η = 50` codeunits"

        @testset "limit < 1" begin
            @test_throws DomainError _wraplines(str, 0)
            @test_throws DomainError _wraplines(str, -1)
        end

        @testset "limit == 1" begin
            wrapped = _wraplines(str, 1)
            @test wrapped isa String
            @test startswith(wrapped, "T\nh\ni\ns\n\ns\ne")
        end

        @testset "limit >= ncodeunits" begin
            for limit in (50, 99)
                wrapped = _wraplines(str, limit)
                @test wrapped isa String
                @test wrapped == str
            end
        end

        @testset "1 < limit < ncodeunits" begin
            @test _wraplines(str, 20) == """
                This sentence
                contains exactly
                `η = 50` codeunits"""
            @test _wraplines(str, 25) == """
                This sentence contains
                exactly `η = 50`
                codeunits"""
            @test _wraplines(str, 30) == """
                This sentence contains
                exactly `η = 50` codeunits"""
        end

        @testset "trailing whitespace is stripped" begin
            str = "16charactersthen    fourspaces "
            @test _wraplines(str, 16) == "16charactersthen\n    fourspaces"
            @test _wraplines(str, 17) == "16charactersthen\n   fourspaces"
            @test _wraplines(str, 18) == "16charactersthen\n  fourspaces"
        end

        @testset "has default `limit=92` argument" begin
            str = string(
                "Lorem ipsum dolor sit amet, consectetur adipiscing elit, ",
                "sed do eiusmod tempor incididunt ut labore et dolore magna aliqua."
            )
            @test _wraplines(str) == _wraplines(str, 92)
        end

        @testset "optional `delim` keyword" begin
            str = string(
                "- Lorem ipsum dolor sit amet, consectetur adipiscing elit, ",
                "sed do eiusmod tempor incididunt ut labore et dolore magna aliqua."
            )
            @test _wraplines(str, 50; delim="\n  ") == """
                - Lorem ipsum dolor sit amet, consectetur
                  adipiscing elit, sed do eiusmod tempor incididunt
                  ut labore et dolore magna aliqua."""
        end
    end
end


text:
module TestMeta
using Test

using Base.Meta: ParseError

using Catlab.Meta

strip_all(expr) = strip_lines(expr, recurse=true)
parse_fun(expr) = parse_function(strip_all(expr))

# Function generation
@test (generate_function(JuliaFunction(:(f(x,y)))) ==
       strip_all(:(function f(x,y) end)))
@test (generate_function(JuliaFunction(:(f(x::Int,y::Int)), :Int)) ==
       strip_all(:(function f(x::Int,y::Int)::Int end)))
@test (generate_function(JuliaFunction(:(f(x)), :Bool, :(isnothing(x)))) ==
       strip_all(:(function f(x)::Bool isnothing(x) end)))

fun_with_docstring_expr = quote
  """Is nothing"""
  function f(x)::Bool
    isnothing(x)
  end
end
@test (strip_all(generate_function(
        JuliaFunction(:(f(x)), :Bool, :(isnothing(x)), "Is nothing"))) ==
       strip_all(fun_with_docstring_expr).args[1])

# Function parsing
@test_throws ParseError parse_fun(:(f(x,y)))
@test (parse_fun(:(function f(x,y) x end)) ==
       JuliaFunction(:(f(x,y)), nothing, quote x end))

@test parse_fun((quote
  """ My docstring
  """
  function f(x,y) x end
end).args[1]) == JuliaFunction(:(f(x,y)), nothing, quote x end, " My docstring\n")

@test (parse_fun(:(function f(x::Int,y::Int)::Int x end)) ==
       JuliaFunction(:(f(x::Int,y::Int)), :Int, quote x end))

@test (parse_fun(:(f(x,y) = x)) ==
       JuliaFunction(:(f(x,y)), nothing, quote x end))

sig = JuliaFunctionSig(:f, [:Int,:Int])
@test parse_function_sig(:(f(x::Int,y::Int))) == sig
@test parse_function_sig(:(f(::Int,::Int))) == sig
@test parse_function_sig(:(f(x,y))) == JuliaFunctionSig(:f, [:Any,:Any])

# Type transformations
bindings = Dict((:r => :R, :s => :S, :t => :T))
@test replace_symbols(bindings, :(foo(x::r,y::s)::t)) == :(foo(x::R,y::S)::T)
@test replace_symbols(bindings, :(foo(xs::Vararg{r}))) == :(foo(xs::Vararg{R}))

end


text:
"""
    macro @agent(name, varargs...)

Basic Macro to create an instance of Agent. It generates an Agent type with the characteristics specified in its arguments.

The only compulsory input is the name of the agent.
 - **dims**: Dimensions.

Then you can add an arbitrary set of additional parameters and rules for the agent.

# Parameter Types

 - **:Local**
Local parameter with Float signature. It can represent the individual mass of each agent, its radius...
It is declared individually for each agent so it can be updated with UpdateLocal types.
Has to be declared as:

    localParameter::Local
    [localParameter1,localParameter2,localParameter3,...]::Local

 - **:Identity**
Parameter with Int signature that can be used to identify individually each agent in the model, its an agent type...
It is declared individually for each agent so it can be updated with UpdateLocal types.
Has to be declared as:

    idName::Identity
    [idName1,idName2,idName3,...]::Identity

 - **:Global**
Global parameter with Float signature. It can represent the temperature felt by all the agents...
It is declared globally for all agents and can be updated with UpdateGlobal.
Has to be declared as:

    globalParameter::Global
    [globalParameter1,globalParameter2,globalParameter3,...]::Global

 - **:GlobalArray**
Global array parameter with Float signature. It can represent the force or interaction matrix between different type of agents...
It is declared globally for all agents and can be updated with UpdateGlobal.
Has to be declared as:

    globalArrayParameter::GlobalArray
    [globalArrayParameter1,globalArrayParameter2,...]::GlobalArray

- **:Medium**
    Medium parameter with Float signature. These paremeters are the parameters of a medium in which the agents are embeded.
    Has to be declared as:

        mediumParameter::Medium
        [mediumParameter1,mediumParameter2,...]::Medium

# Update rules

 - **Equation**
Type for declaring dynamical equations for the variables.
Has to be declared as:

    Equation =
        begin
            ... code including the differential equations ...
        end

 - **UpdateGlobal**
Type for declaring global update rules for the global parameters and global arrays.
Has to be declared as:

    UpdateGlobal =
        begin
            ... code including the update rules ...
        end

 - **UpdateLocal**
Type for declaring local update rules for the local parameters and ids.
Has to be declared as:

    UpdateLocal =
        begin
            ... code including the update rules ...
        end

 - **UpdateInteraction**
Type for declaring interaction update rules for the interaction parameters.
The variables declared in these updates will be updated during the integration steps.
Has to be declared as:

    UpdateInteraction =
        begin
            ... code including the update rules ...
        end

 - **UpdateLocalInteraction**
Type for declaring interaction update rules for the local parameters and ids.
The variables declared in these updates will be updated just once per time step instead of once per integration point.
For integrators that evaluate just once as the Euler integrator, UpdateInteraction and UpdateLocalInteraction are equivalent.
Has to be declared as:

    UpdateLocalInteraction =
    begin
            ... code including the update rules ...
        end

 - **UpdateMedium**
Type for declaring the equations of motion that regulate the medium.
Has to be declared as:

    UpdateMedium =
        begin
                ... code including the update rules ...
            end

 - **UpdateMediumInteraction**
Type for declaring the interaction between the agents and the medium.
The agents are always modelled as delta sources of interaction with the medium located at their coordinates (x,y,z).
Has to be declared as:

UpdateMediumInteraction =
    begin
            ... code including the update rules ...
        end

# Custom made blocks of code

Any function which returns a piece of code with the structue as the above declared is valid for declaration. This allows repetitive algorithms to be reused and allow a fluent customization.

```
function localWithIncrements(variable,step)
    return begin
        \$variable::Local

        UpdateLocal = \$variable += \$step
    end
end

@agent(
    Agent, #Name of the agent

    l1::Local,

    localWithIncrement(:b)
)

```
"""
macro agent(dims, varargs...)

    m = Agent()
    m.dims = dims

    if dims == 0
        nothing
    elseif dims == 1
        push!(m.declaredSymbols["Local"],:x)
    elseif dims == 2
        push!(m.declaredSymbols["Local"],:x,:y)
    elseif dims == 3
        push!(m.declaredSymbols["Local"],:x,:y,:z)
    else
        error("Dims has to be an integer between 0 and 3.")
    end
    push!(m.declaredSymbols["Identity"],:id)

    #Add all contributions
    for ii in varargs
        if typeof(ii) != Expr
            error(ii, " should be and expression containing at least a name and a type. Example variable::Variable.")
        end

        if ii.head == :call
            t = eval(ii)
            if typeof(t) != Expr
                error("Function ", ii, " should return a piece of code to be added to the agent.")
            elseif t.head != :block
                t = quote $t end
            end

            t = t.args
        else
            t = [ii]
        end

        for i in t
            if typeof(i) == LineNumberNode
                nothing
            elseif i.head == :(::)
                name = i.args[1]
                type = i.args[2]

                if !(type in VALID_TYPES)
                    error(i, " should be and expression containing at least a name and a type. Example variable::Variable.")
                elseif typeof(name) == Expr
                    if !(name.head == :vect)
                        error(name, " should be a symbol or an array of symbols. Check how to declare ", type, " types.")
                    else
                        names = Symbol[]
                        for j in name.args
                            if typeof(j) != Symbol
                                error(j, " in ", i, " should be a symbol.  Check how to declare ", type, " types.")
                            end
                            push!(names,j)
                        end
                        name = names
                    end
                end

                checkDeclared_(m,name)
                if typeof(name) == Array{Symbol,1}
                    append!(m.declaredSymbols[string(type)],name)
                else
                    push!(m.declaredSymbols[string(type)],name)
                end

            elseif i.head == :(=)
                if i.args[1] in VALID_UPDATES
                    if !(string(i.args[1]) in keys(m.declaredUpdates))
                        if i.args[2].head == :block
                            m.declaredUpdates[string(i.args[1])] = i.args[2]
                        else
                            m.declaredUpdates[string(i.args[1])] = quote $(i.args[2]) end
                        end
                    elseif i.args[2].head == :block
                        append!(m.declaredUpdates[string(i.args[1])].args,i.args[2].args)
                    end
                else
                    error(i.args[1], " is not a valid type.")
                end
            else
                error(i, " is not an understood rule or variable declaration. Error in ", ii)
            end
        end
    end

    return m
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
@testset "965.univalued-binary-tree.jl" begin
    @test is_unival_tree(TreeNode(1, TreeNode(1, TreeNode(1), TreeNode(1)), TreeNode(1, nothing, TreeNode(1)))) == true
    @test is_unival_tree(TreeNode(1, TreeNode(1, TreeNode(5), TreeNode(1)), TreeNode(1))) == false
end


text:
module Collocation

export
    CollocationMatrix,
    collocation_points,
    collocation_points!,
    collocation_matrix,
    collocation_matrix!

using ..BSplines
using ..DifferentialOps
using ..Recombinations:
    num_constraints,
    RecombinedBSplineBasis  # for Documenter

using BandedMatrices
using SparseArrays

include("matrix.jl")
include("points.jl")

"""
    collocation_matrix(
        B::AbstractBSplineBasis,
        x::AbstractVector,
        [deriv::Derivative = Derivative(0)],
        [MatrixType = CollocationMatrix{Float64}];
        clip_threshold = eps(eltype(MatrixType)),
    )

Return collocation matrix mapping B-spline coefficients to spline values at the
collocation points `x`.

# Extended help

The matrix elements are given by the B-splines evaluated at the collocation
points:

```math
C_{ij} = b_j(x_i) \\quad \\text{for} \\quad
i ∈ [1, N_x] \\text{ and } j ∈ [1, N_b],
```

where `Nx = length(x)` is the number of collocation points, and
`Nb = length(B)` is the number of B-splines in `B`.

To obtain a matrix associated to the B-spline derivatives, set the `deriv`
argument to the order of the derivative.

Given the B-spline coefficients ``\\{u_j, 1 ≤ j ≤ N_b\\}``, one can recover the
values (or derivatives) of the spline at the collocation points as `v = C * u`.
Conversely, if one knows the values ``v_i`` at the collocation points,
the coefficients ``u`` of the spline passing by the collocation points may be
obtained by inversion of the linear system `u = C \\ v`.

The `clip_threshold` argument allows one to ignore spurious, negligible values
obtained when evaluating B-splines. These values are typically unwanted, as they
artificially increase the number of elements (and sometimes the bandwidth) of
the matrix.
They may appear when a collocation point is located on a knot.
By default, `clip_threshold` is set to the machine epsilon associated to the
matrix element type (see `eps` in the Julia documentation).
Set it to zero to disable this behaviour.

## Matrix types

The `MatrixType` optional argument allows to select the type of returned matrix.

Due to the compact support of B-splines, the collocation matrix is
[banded](https://en.wikipedia.org/wiki/Band_matrix) if the collocation points
are properly distributed. Therefore, it makes sense to store it in a
`BandedMatrix` (from the
[BandedMatrices](https://github.com/JuliaMatrices/BandedMatrices.jl) package),
as this will lead to memory savings and especially to time savings if
the matrix needs to be inverted.

### Supported matrix types

- `CollocationMatrix{T}`: thin wrapper over `BandedMatrix{T}`, with efficient LU
  factorisations without pivoting (see [`CollocationMatrix`](@ref) for details).
  This option performs much better than sparse matrices for inversion of linear
  systems.
  On the other hand, for matrix-vector or matrix-matrix
  multiplications, `SparseMatrixCSC` may perform better, especially when using OpenBLAS (see
  [BandedMatrices issue](https://github.com/JuliaMatrices/BandedMatrices.jl/issues/110)).
  May fail with an error for non-square matrix shapes, or if the distribution of
  collocation points is not adapted. In these cases, the effective
  bandwidth of the matrix may be larger than the expected bandwidth.

- `SparseMatrixCSC{T}`: regular sparse matrix; correctly handles any matrix
  shape.

- `Matrix{T}`: a regular dense matrix. Generally performs worse than the
  alternatives, especially for large problems.

See also [`collocation_matrix!`](@ref).
"""
function collocation_matrix(
        B::AbstractBSplineBasis, x::AbstractVector,
        deriv::Derivative = Derivative(0),
        ::Type{MatrixType} = CollocationMatrix{Float64};
        kwargs...,
    ) where {MatrixType}
    Nx = length(x)
    Nb = length(B)
    C = allocate_collocation_matrix(MatrixType, (Nx, Nb), order(B); kwargs...)
    collocation_matrix!(C, B, x, deriv; kwargs...)
end

collocation_matrix(B, x, ::Type{M}; kwargs...) where {M} =
    collocation_matrix(B, x, Derivative(0), M; kwargs...)

allocate_collocation_matrix(::Type{M}, dims, k) where {M <: AbstractMatrix} =
    M(undef, dims...)

allocate_collocation_matrix(::Type{SparseMatrixCSC{T}}, dims, k) where {T} =
    spzeros(T, dims...)

function allocate_collocation_matrix(
        ::Type{M}, dims, k,
    ) where {M <: CollocationMatrix}
    # Number of upper and lower diagonal bands.
    #
    # The matrices **almost** have the following number of upper/lower bands (as
    # cited sometimes in the literature):
    #
    #  - for even k: Nb = (k - 2) / 2 (total = k - 1 bands)
    #  - for odd  k: Nb = (k - 1) / 2 (total = k bands)
    #
    # However, near the boundaries U/L bands can become much larger.
    # Specifically for the second and second-to-last collocation points (j = 2
    # and N - 1). For instance, the point j = 2 sees B-splines in 1:k, leading
    # to an upper band of size k - 2.
    #
    # TODO is there a way to reduce the number of bands??
    bands = collocation_bandwidths(k)
    T = eltype(M)
    data = BandedMatrix{T}(undef, dims, bands)
    CollocationMatrix(data) :: M
end

collocation_bandwidths(k) = (k - 2, k - 2)

"""
    collocation_matrix!(
        C::AbstractMatrix{T}, B::AbstractBSplineBasis, x::AbstractVector,
        [deriv::Derivative = Derivative(0)]; clip_threshold = eps(T))

Fill preallocated collocation matrix.

See [`collocation_matrix`](@ref) for details.
"""
function collocation_matrix!(
        C::AbstractMatrix{T}, B::AbstractBSplineBasis, x::AbstractVector,
        deriv::Derivative = Derivative(0);
        clip_threshold = eps(T),
    ) where {T <: AbstractFloat}
    if axes(C, 1) != axes(x, 1) || size(C, 2) != length(B)
        throw(ArgumentError("wrong dimensions of collocation matrix"))
    end

    fill!(C, 0)
    b_lo, b_hi = max_bandwidths(C)

    @inbounds for i ∈ axes(C, 1)
        jlast, bs = evaluate_all(B, x[i], deriv, T)  # = (b_{jlast + 1 - k}(x), …, b_{jlast}(x))

        for (δj, bj) ∈ pairs(bs)
            j = jlast + 1 - δj

            # Skip very small values (and zeros).
            # This is important for SparseMatrixCSC, which also stores explicit
            # zeros.
            abs(bj) <= clip_threshold && continue

            if (i > j && i - j > b_lo) || (i < j && j - i > b_hi)
                # This will cause problems if C is a BandedMatrix, and (i, j) is
                # outside the allowed bands. This may be the case if the collocation
                # points are not properly distributed.
                @warn "Non-zero value outside of matrix bands: b[$j](x[$i]) = $bj"
            end

            C[i, j] = bj
        end
    end

    C
end

# Maximum number of bandwidths allowed in a matrix container.
max_bandwidths(A::BandedMatrix) = bandwidths(A)
max_bandwidths(A::CollocationMatrix) = max_bandwidths(parent(A))
max_bandwidths(A::AbstractMatrix) = size(A) .- 1

end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
"""
    AdjustBrightness(δ = 0.2)
    AdjustBrightness(distribution)

Adjust the brightness of an image by a factor chosen uniformly
from `f ∈ [1-δ, 1+δ]` by multiplying each color channel by `f`.

You can also pass any `Distributions.Sampleable` from which the
factor is selected.

## Example

{cell=AdjustBrightness}
```julia
using DataAugmentation, TestImages

item = Image(testimage("lighthouse"))
tfm = AdjustBrightness(0.2)
titems = [apply(tfm, item) for _ in 1:8]
showgrid(titems; ncol = 4, npad = 16)
```
"""
struct AdjustBrightness{S<:Sampleable} <: Transform
    dist::S
end

AdjustBrightness(f::Real) = AdjustBrightness(Uniform(max(0, 1 - f), 1 + f))

getrandstate(tfm::AdjustBrightness) = rand(tfm.dist)

function apply(tfm::AdjustBrightness, item::Image; randstate = getrandstate(tfm))
    factor = randstate
    return setdata(item, adjustbrightness(itemdata(item), factor))
end

function apply!(buf, tfm::AdjustBrightness, item::Image; randstate = getrandstate(tfm))
    factor = randstate
    adjustbrightness!(itemdata(buf), itemdata(item), factor)
    return buf
end

function adjustbrightness(img, factor)
    return adjustbrightness!(copy(img), factor)
end

adjustbrightness!(img, factor) = adjustbrightness!(img, img, factor)

# TODO: add methods for non-RGB/Gray images
function adjustbrightness!(dst::AbstractArray{U}, img::AbstractArray{T}, factor) where {T, U}
    map!(dst, img) do x
        convert(U, clamp01(x * factor))
    end
end

##
"""
    AdjustContrast(factor = 0.2)
    AdjustContrast(distribution)

Adjust the contrast of an image by a factor chosen uniformly
from `f ∈ [1-δ, 1+δ]`.

Pixels `c` are transformed `c + μ*(1-f)` where `μ` is the mean color
of the image.

You can also pass any `Distributions.Sampleable` from which the
factor is selected.

## Example

{cell=AdjustBrightness}
```julia
using DataAugmentation, TestImages

item = Image(testimage("lighthouse"))
tfm = AdjustContrast(0.2)
titems = [apply(tfm, item) for _ in 1:8]
showgrid(titems; ncol = 4, npad = 16)
```
"""
struct AdjustContrast{S<:Sampleable} <: Transform
    dist::S
end

AdjustContrast(f::Real) = AdjustContrast(Uniform(max(0, 1 - f), 1 + f))

getrandstate(tfm::AdjustContrast) = rand(tfm.dist)

function apply(tfm::AdjustContrast, item::Image; randstate = getrandstate(tfm))
    factor = randstate
    return setdata(item, adjustcontrast(itemdata(item), factor))
end

function apply!(buf, tfm::AdjustContrast, item::Image; randstate = getrandstate(tfm))
    factor = randstate
    adjustcontrast!(itemdata(buf), itemdata(item), factor)
    return buf
end

function adjustcontrast(img, factor)
    return adjustcontrast!(copy(img), factor)
end

# TODO: add methods for non-RGB/Gray images
function adjustcontrast!(dst::AbstractArray{U}, img::AbstractArray{T}, factor) where {T, U}
    μ = mean(img)
    map!(dst, img) do x
        convert(U, clamp01(x + μ * (1 - factor)))
    end
end

adjustcontrast!(img, factor) = adjustcontrast!(img, img, factor)


text:
function AESCTR(plaintext, iv::Array{UInt8, 1}, key::AbstractAESKey, cache::AbstractAESCache)
    len = length(plaintext)
    ciphertext = similar(Array{UInt8, 1}, len)
    AESCTR!(ciphertext, plaintext, iv, key, cache)
end

function AESCTR!(ciphertext, plaintext, iv::Array{UInt8, 1}, key::AbstractAESKey, cache::AbstractAESCache)
    len = length(plaintext)
    pad = 16 - (len % 16)
    iters = Int((len + pad) / 16)
    counter = cache.modecache
    counter .= iv
    ciphertextblock = AESBlock(ciphertext)
    plaintextblock = AESBlock(plaintext)
    for i in 1:(iters-1)
        AESEncryptBlock!(ciphertextblock, counter, key, cache)
        @. ciphertextblock ⊻= plaintextblock
        increment!(counter)
        increment!(ciphertextblock)
        increment!(plaintextblock)
    end
    # for the last iteration, we dont need counter no more
    # we can use the inplace property of AES
    AESEncryptBlock!(counter, counter, key, cache)
    offset = 16(iters-1)
    for i in 1:(16-pad)
        ciphertext[offset+i] = UInt8(plaintext[offset+i]) ⊻ counter[i]
    end
    ciphertext
end

# Encryption is same as decryption in CTR mode for a give Key and IV
const AESCTR_D = AESCTR
const AESCTR_D! = AESCTR!


text:
using CellFishing: CellFishing, HammingIndexes
using Test
using Random: Random, shuffle
using SparseArrays: SparseMatrixCSC
using Distributions: Poisson

resetseed!() = Random.seed!(1234)

@testset "BitVectors" begin
    for T in [CellFishing.BitVec64,
              CellFishing.BitVec128,
              CellFishing.BitVec256,
              CellFishing.BitVec512]
        @test T <: CellFishing.BitVec
        @test convert(T, 0) isa T
        @test T(0) isa T
        @test zero(T) === T(0)
        @test T(false) === T(0)
        @test T(true) === T(1)
        @test (T(0b0011) & T(0b0101)) === T(0b0001)
        @test (T(0b0011) | T(0b0101)) === T(0b0111)
        @test (T(0b0011) ⊻ T(0b0101)) === T(0b0110)
        @test T(0b0010) << 1 === T(0b0100)
        @test T(0b0010) >> 1 === T(0b0001)
        @test T(0b0101)[1] === true
        @test T(0b0101)[2] === false
        @test T(0b0101) % UInt32 === UInt32(0b0101)
        @test T(0b0101) % UInt64 === UInt64(0b0101)
        @test count_ones(T(0b0101)) === 2
        @test rand(T) isa T
    end
end

@testset "HammingIndexes" begin
    linear = HammingIndexes.LinearSearch()
    mindex = HammingIndexes.MultiIndexSearch()
    resetseed!()
    data = rand(UInt64, 1000)
    db = HammingIndexes.HammingIndex(data)
    for r in 0:16
        ok = true
        for i in 1:length(data)
            x = HammingIndexes.findwithin(linear, r, data[i], db)
            y = HammingIndexes.findwithin(mindex, r, data[i], db)
            ok &= x == y
        end
        @test ok
    end
    for k in 0:16
        ok = true
        for i in 1:length(data)
            x = HammingIndexes.findknn(linear, k, data[i], db)
            y = HammingIndexes.findknn(mindex, k, data[i], db)
            # indexes may be different
            ok &= x.distances == y.distances
        end
        @test ok
    end

    resetseed!()
    data = rand(UInt64, 1_000_000)
    db = HammingIndexes.HammingIndex(data)
    queries = data[rand(1:length(data), 100)]
    for r in [16, 18, 20]
        ok = true
        for q in queries
            x = HammingIndexes.findwithin(linear, r, q, db)
            y = HammingIndexes.findwithin(mindex, r, q, db)
            ok &= x == y
        end
        @test ok
    end
    for k in 0:16
        ok = true
        for q in queries
            x = HammingIndexes.findknn(linear, k, q, db)
            y = HammingIndexes.findknn(mindex, k, q, db)
            # indexes may be different
            ok &= x.distances == y.distances
        end
        @test ok
    end
end

@testset "CMatrix" begin
    m, n = 100, 200
    resetseed!()
    counts = rand(0:100, m, n)
    # with the default parameters
    M = CellFishing.CMatrix(counts)
    @test M isa AbstractMatrix{Int}
    @test all(M .== counts)
    for j in 1:n
        @test M[:,j] == counts[:,j]
    end
    # with specific parameters
    M = CellFishing.CMatrix(
        counts,
        level=9,
        shuffle=false,
        compressor="lz4hc",
        blocksize=m,
        nthreads=2)
    @test M isa AbstractMatrix{Int}
    @test all(M .== counts)
    for j in 1:n
        @test M[:,j] == counts[:,j]
    end
end

@testset "Features" begin
    m = 100
    names = string.("feature:", 1:m)
    features = CellFishing.Features(names)
    @test CellFishing.nfeatures(features) == m
    @test CellFishing.nselected(features) == 0

    # add features
    features = CellFishing.addfeatures(features, ["feature:10", "feature:32", "feature:99"])
    @test CellFishing.nfeatures(features) == m
    @test CellFishing.nselected(features) == 3
    @test CellFishing.selectedfeatures(features) == ["feature:10", "feature:32", "feature:99"]
    @test_throws ArgumentError CellFishing.addfeatures(features, ["foo", "bar"])

    # drop features
    features = CellFishing.dropfeatures(features, ["feature:32"])
    @test CellFishing.nfeatures(features) == m
    @test CellFishing.nselected(features) == 2
    @test CellFishing.selectedfeatures(features) == ["feature:10", "feature:99"]
    @test_throws ArgumentError CellFishing.dropfeatures(features, ["foo", "bar"])
end

@testset "CellIndex" begin
    resetseed!()
    m, n = 100, 200
    counts = rand(0:1000, m, n)
    featurenames = string.("feature:", 1:m)
    features = CellFishing.selectfeatures(counts, featurenames, n_min_features=m)
    index = CellFishing.CellIndex(counts, features, metadata=string.("cell:", 1:n))
    @test CellFishing.nbits(index) == 128
    @test CellFishing.ncells(index) == n
    @test occursin("CellFishing.CellIndex(<#cells=200, hash=128×4>)", sprint(show, index))
    #@test index.featurenames[1] == "feature:1"
    @test index.metadata[1] == "cell:1"
    @test index.counts === nothing
    @test occursin("CellIndex", sprint(show, index))
    mktempdir() do tmpdir
        tmpfile = joinpath(tmpdir, "index")
        CellFishing.save(tmpfile, index)
        @test CellFishing.load(tmpfile) isa CellFishing.CellIndex
        @test CellFishing.load(tmpfile).counts === nothing
    end
    # accept Union{T,Missing} values (but no missings)
    index = CellFishing.CellIndex(Matrix{Union{Int,Missing}}(counts), features, metadata=string.("cell:", 1:n))
    @test index isa CellFishing.CellIndex
    # keep raw counts
    index = CellFishing.CellIndex(counts, features, keep_counts=true)
    @test index.counts == counts[features.selected,:]
    mktempdir() do tmpdir
        tmpfile = joinpath(tmpdir, "index")
        CellFishing.save(tmpfile, index)
        @test CellFishing.load(tmpfile) isa CellFishing.CellIndex
        @test CellFishing.load(tmpfile).counts !== nothing
    end

    resetseed!()
    m, n = 100, 200
    counts = rand(0:1000, m, n)
    featurenames = string.("feature:", 1:m)
    n_all = n_ok = 0
    for MatType in [Matrix, SparseMatrixCSC]
        for n_bits in [64, 128, 256, 512],
            superbit in [1, 8, 17],
            n_dims in [64, m-1],
            randomize in false:true,
            normalize in false:true,
            standardize in false:true,
            index in false:true
            features = CellFishing.selectfeatures(MatType(counts), featurenames, n_min_features=m)
            idx = CellFishing.CellIndex(
                MatType(counts), features,
                n_bits=n_bits, superbit=superbit,
                n_dims=n_dims, randomize=randomize,
                normalize=normalize, standardize=standardize,
                index=index)
            perm = shuffle(1:m)
            neighbors = CellFishing.findneighbors(1, CellFishing.ExpressionMatrix(MatType(counts)[perm,:], featurenames[perm]), idx)
            totaldist::Int = 0
            for i in 1:10
                j = neighbors.indexes[1,i]
                totaldist += neighbors.hammingdistances[1,i]
                n_all += 1
                n_ok += i == j
            end
            @test totaldist ≤ 1
            S = CellFishing.similarities(neighbors, idx)
            @test all(-1 .≤ S .≤ 1)
            neighbors = CellFishing.findneighbors(3, [1, 2, 3, 5, 7], idx)
            @test neighbors.indexes[1,:] == [1, 2, 3, 5, 7]
            @test all(neighbors.hammingdistances[1,:] .== 0)
        end
    end
    @test n_ok / n_all ≈ 1.0
end

@testset "DEG" begin
    resetseed!()
    m, n = 100, 3000
    λ = rand(m) * 100
    counts = vcat([rand(Poisson(λ[i]), n)' for i in 1:m]...)
    featurenames = string.("feature:", 1:m)
    features = CellFishing.selectfeatures(counts, featurenames, n_min_features=m)
    index = CellFishing.CellIndex(counts, features, keep_counts=true)
    for _ in 1:10
        n = 1000
        counts = vcat([rand(Poisson(λ[i]), n)' for i in 1:m]...)
        neighbors = CellFishing.findneighbors(1, counts, featurenames, index)
        degs = CellFishing.finddegs(counts, featurenames, neighbors.indexes[1,:], index)
        #@show sum(degs.negatives .< -3), sum(degs.positives .< -3)
        @test 10 < sum(degs.negatives .< -3) < 70
        @test 10 < sum(degs.positives .< -3) < 70
    end
end