text:
"""
    PyObject(x)

Convert `x` to a Python object.

This is the default type returned by most API functions.

It supports attribute access (`x.attr`), indexing (`x[i,j]`, `delete!(x, i, j)`),
calling `x(a, b, kw=c)`, `length`, iteration (yielding `PyObject`s), comparisons to other
`PyObject`s (`x == y`), and arithmetic with other `PyObject`s and `Number`s (`x + y`, `x * 3`).

Note that comparisons between `PyObject`s (`==`, `≠`, `≤`, etc.) return `PyObject`,
except `isequal` and `isless` which return `Bool`. Use [`pyeq`](@ref) and friends to
return `Bool`.

A number of special properties are also defined for convenience to convert the object to
another type. To avoid clashes with object attributes, they all have the prefix `jl!`.
- `x.jl!(T)` is `pyconvert(T, x)`
- `x.jl!i` is `pyconvert(Int, x)`
- `x.jl!b` is `pytruth(x)`
- `x.jl!s` is `pystr(String, x)`.
- `x.jl!r` is `pyrepr(String, x)`.
- `x.jl!f` is `pyconvert(Float64, x)`
- `x.jl!c` is `pyconvert(Complex{Float64}, x)`
- `x.jl!iter(T=PyObject)` is `PyIterable{T}(x)`
- `x.jl!list(T=PyObject)` is `PyList{T}(x)`
- `x.jl!set(T=PyObject)` is `PySet{T}(x)`
- `x.jl!dict(K=PyObject, V=PyObject)` is `PyDict{K,V}(x)`
- `x.jl!io(...)` is `PyIO(x; ...)`
- `x.jl!pandasdf(...)` is `PyPandasDataFrame(x; ...)`
- `x.jl!buffer(...)` is `PyBuffer(x, ...)`
- `x.jl!array(...)` is `PyArray{...}(x)`
- `x.jl!vector(...)` is `PyVector{...}(x)`
- `x.jl!matrix(...)` is `PyMatrix{...}(x)`
"""
mutable struct PyObject
    ref::PyRef
    make::Any
    PyObject(::Val{:nocopy}, o::PyRef) = new(o)
    PyObject(o) = new(PyRef(o))
    PyObject(::Val{:lazy}, mk) = new(PyRef(), mk)
end
export PyObject

ispyreftype(::Type{PyObject}) = true
pyptr(o::PyObject) = begin
    ref = getfield(o, :ref)
    ptr = ref.ptr
    if isnull(ptr)
        val = try
            getfield(o, :make)()
        catch err
            C.PyErr_SetString(C.PyExc_Exception(), "Error retrieving object value: $err")
            return ptr
        end
        ptr = ref.ptr = C.PyObject_From(val)
    end
    ptr
end
Base.unsafe_convert(::Type{CPyPtr}, o::PyObject) = checknull(pyptr(o))
pynewobject(p::Ptr, check::Bool = false) =
    (check && isnull(p)) ? pythrow() : PyObject(Val(:nocopy), pynewref(p))
pyborrowedobject(p::Ptr, check::Bool = false) =
    (check && isnull(p)) ? pythrow() : PyObject(Val(:nocopy), pyborrowedref(p))
pylazyobject(mk) = PyObject(Val(:lazy), mk)

C.PyObject_TryConvert__initial(o, ::Type{PyObject}) = C.putresult(pyborrowedobject(o))

Base.convert(::Type{PyObject}, x::PyObject) = x
Base.convert(::Type{Any}, x::PyObject) = x
Base.convert(::Type{T}, x::PyObject) where {T} = x isa T ? x : pyconvert(T, x)
Base.convert(::Type{PyObject}, x) = PyObject(x)

### Cache some common values

const _pynone = pylazyobject(() -> pynone(PyRef))
pynone(::Type{PyObject}) = _pynone

const _pytrue = pylazyobject(() -> pybool(PyRef, true))
const _pyfalse = pylazyobject(() -> pybool(PyRef, false))
pybool(::Type{PyObject}, x::Bool) = x ? _pytrue : _pyfalse

### IO

Base.string(o::PyObject) = pystr(String, o)

Base.print(io::IO, o::PyObject) = print(io, string(o))

Base.show(io::IO, o::PyObject) = begin
    s = pyrepr(String, o)
    if get(io, :typeinfo, Any) == typeof(o)
        print(io, s)
    elseif startswith(s, "<") && endswith(s, ">")
        print(io, "<py ", s[2:end])
    else
        print(io, "<py ", s, ">")
    end
end

function Base.show(io::IO, ::MIME"text/plain", o::PyObject)
    prefix = get(io, :typeinfo, Any) != typeof(o)
    h, w = displaysize(io)
    h -= 3
    x = try
        pystr(String, pypprintmodule().pformat(o, width = w))
    catch
        pyrepr(String, o)
    end
    if get(io, :limit, true)
        if '\n' ∈ x
            # multiple lines
            # each one is truncated to one screen width
            if prefix
                print(io, "py:")
                h -= 1
            end
            xs = split(x, '\n')
            printlines(xs, nl = true) =
                for (i, x) in enumerate(xs)
                    (nl || i > 1) && print(io, '\n')
                    if length(x) ≤ w
                        print(io, x)
                    else
                        print(io, x[1:nextind(x, 0, w - 1)], '…')
                    end
                end
            if length(xs) ≤ h
                # all lines fit on screen
                printlines(xs, prefix)
            else
                # too many lines, skip the middle ones
                h -= 1
                h2 = cld(h, 2)
                h3 = (length(xs) + 1) - (h - h2)
                printlines(xs[1:h2], prefix)
                linelen = min(
                    checkbounds(Bool, xs, h2) ? length(xs[h2]) : 0,
                    checkbounds(Bool, xs, h3) ? length(xs[h3]) : 0,
                    w,
                )
                msg = "... [skipping $(h3-h2-1) lines] ..."
                pad = fld(linelen - length(msg), 2)
                print(io, "\n", pad > 0 ? " "^pad : "", msg)
                printlines(xs[h3:end])
            end
        elseif length(x) ≤ (prefix ? w - 4 : w)
            # one short line
            print(io, prefix ? "py: " : "", x)
            return
        else
            # one long line
            if prefix
                println(io, "py:")
                h -= 1
            end
            a = h * w
            if length(x) ≤ a
                # whole string fits on screen
                print(io, x)
            else
                # too long, skip the middle
                h -= 1
                h2 = cld(h, 2)
                i2 = nextind(x, 0, h2 * w)
                i3 = prevind(x, ncodeunits(x) + 1, (h - h2) * w)
                println(io, x[1:i2])
                println(
                    io,
                    " ... [skipping $(length(x[nextind(x,i2):prevind(x,i3)])) characters] ...",
                )
                print(io, x[i3:end])
            end
        end
    else
        print(io, prefix ? "py: " : "", x)
    end
end

Base.show(io::IO, mime::MIME, o::PyObject) = _py_mime_show(io, mime, o)
Base.showable(mime::MIME, o::PyObject) = _py_mime_showable(mime, o)

### PROPERTIES

Base.getproperty(o::PyObject, k::Symbol) =
    if k == :jl!
        ((::Type{T}) where {T}) -> pyconvert(T, o)
    elseif k == :jl!i
        pyconvert(Int, o)
    elseif k == :jl!b
        pytruth(o)
    elseif k == :jl!s
        pystr(String, o)
    elseif k == :jl!r
        pyrepr(String, o)
    elseif k == :jl!f
        pyconvert(Cdouble, o)
    elseif k == :jl!c
        pyconvert(Complex{Cdouble}, o)
    elseif k == :jl!iter
        ((::Type{T} = PyObject) where {T}) -> PyIterable{T}(o)
    elseif k == :jl!list
        ((::Type{T} = PyObject) where {T}) -> PyList{T}(o)
    elseif k == :jl!set
        ((::Type{T} = PyObject) where {T}) -> PySet{T}(o)
    elseif k == :jl!dict
        ((::Type{K} = PyObject, ::Type{V} = PyObject) where {K,V}) -> PyDict{K,V}(o)
    elseif k == :jl!io
        (; opts...) -> PyIO(o; opts...)
    elseif k == :jl!pandasdf
        (; opts...) -> PyPandasDataFrame(o; opts...)
    elseif k == :jl!buffer
        (args...) -> PyBuffer(o, args...)
    elseif k == :jl!array
        f_arr() = PyArray(o)
        f_arr(::Type{T}) where {T} = PyArray{T}(o)
        f_arr(::Type{T}, N::Integer) where {T} = PyArray{T,N}(o)
        f_arr(::Type{T}, N::Integer, ::Type{R}) where {T,R} = PyArray{T,N,R}(o)
        f_arr(::Type{T}, N::Integer, ::Type{R}, M::Bool) where {T,R} = PyArray{T,N,R,M}(o)
        f_arr(::Type{T}, N::Integer, ::Type{R}, M::Bool, L::Bool) where {T,R} =
            PyArray{T,N,R,M,L}(o)
        f_arr
    elseif k == :jl!vector
        f_vec() = PyVector(o)
        f_vec(::Type{T}) where {T} = PyVector{T}(o)
        f_vec(::Type{T}, ::Type{R}) where {T,R} = PyVector{T,R}(o)
        f_vec(::Type{T}, ::Type{R}, M::Bool) where {T,R} = PyVector{T,R,M}(o)
        f_vec(::Type{T}, ::Type{R}, M::Bool, L::Bool) where {T,R} = PyVector{T,R,M,L}(o)
        f_vec
    elseif k == :jl!matrix
        f_mat() = PyMatrix(o)
        f_mat(::Type{T}) where {T} = PyMatrix{T}(o)
        f_mat(::Type{T}, ::Type{R}) where {T,R} = PyMatrix{T,R}(o)
        f_mat(::Type{T}, ::Type{R}, M::Bool) where {T,R} = PyMatrix{T,R,M}(o)
        f_mat(::Type{T}, ::Type{R}, M::Bool, L::Bool) where {T,R} = PyMatrix{T,R,M,L}(o)
        f_mat
    else
        pygetattr(PyObject, o, k)
    end

Base.setproperty!(o::PyObject, k::Symbol, v) = pysetattr(o, k, v)

if hasproperty(Base, :hasproperty)
    Base.hasproperty(o::PyObject, k::Symbol) = pyhasattr(o, k)
else
    Compat.hasproperty(o::PyObject, k::Symbol) = pyhasattr(o, k)
end

function Base.propertynames(o::PyObject)
    # this follows the logic of rlcompleter.py
    @py ```
    def members(o):
        def classmembers(c):
            r = dir(c)
            if hasattr(c, '__bases__'):
                for b in c.__bases__:
                    r += classmembers(b)
            return r
        words = set(dir(o))
        words.discard('__builtins__')
        if hasattr(o, '__class__'):
            words.add('__class__')
            words.update(classmembers(o.__class__))
        return words
    $(words::Vector{Symbol}) = members($o)
    ```
    words
end

### CALL

(f::PyObject)(args...; kwargs...) = pycall(PyObject, f, args...; kwargs...)

### ITERATE & INDEX

Base.IteratorSize(::Type{PyObject}) = Base.SizeUnknown()
Base.eltype(::Type{PyObject}) = PyObject

Base.getindex(o::PyObject, k) = pygetitem(PyObject, o, k)
Base.getindex(o::PyObject, k...) = pygetitem(PyObject, o, k)
Base.setindex!(o::PyObject, v, k) = pysetitem(o, k, v)
Base.setindex!(o::PyObject, v, k...) = pysetitem(o, k, v)
Base.delete!(o::PyObject, k) = pydelitem(o, k)
Base.delete!(o::PyObject, k...) = pydelitem(o, k)

Base.length(o::PyObject) = Int(pylen(o))

Base.iterate(o::PyObject, it::PyRef = pyiter(PyRef, o)) = begin
    vo = C.PyIter_Next(it)
    if !isnull(vo)
        (pynewobject(vo), it)
    elseif C.PyErr_IsSet()
        pythrow()
    else
        nothing
    end
end

Base.in(x, o::PyObject) = pycontains(o, x)
Base.hash(o::PyObject) = reinterpret(UInt, Int(pyhash(o)))

### COMPARISON

Base.:(==)(x::PyObject, y::PyObject) = pyeq(PyObject, x, y)
Base.:(!=)(x::PyObject, y::PyObject) = pyne(PyObject, x, y)
Base.:(<=)(x::PyObject, y::PyObject) = pyle(PyObject, x, y)
Base.:(<)(x::PyObject, y::PyObject) = pylt(PyObject, x, y)
Base.:(>=)(x::PyObject, y::PyObject) = pyge(PyObject, x, y)
Base.:(>)(x::PyObject, y::PyObject) = pygt(PyObject, x, y)
Base.isequal(x::PyObject, y::PyObject) = pyeq(Bool, x, y)
Base.isless(x::PyObject, y::PyObject) = pylt(Bool, x, y)

### ARITHMETIC

Base.zero(::Type{PyObject}) = pyint(0)
Base.one(::Type{PyObject}) = pyint(1)

# unary
Base.:(-)(o::PyObject) = pyneg(PyObject, o)
Base.:(+)(o::PyObject) = pypos(PyObject, o)
Base.abs(o::PyObject) = pyabs(PyObject, o)
Base.:(~)(o::PyObject) = pyinv(PyObject, o)

# binary
Base.:(+)(o1::PyObject, o2::PyObject) = pyadd(PyObject, o1, o2)
Base.:(-)(o1::PyObject, o2::PyObject) = pysub(PyObject, o1, o2)
Base.:(*)(o1::PyObject, o2::PyObject) = pymul(PyObject, o1, o2)
Base.:(/)(o1::PyObject, o2::PyObject) = pytruediv(PyObject, o1, o2)
Base.:fld(o1::PyObject, o2::PyObject) = pyfloordiv(PyObject, o1, o2)
Base.:mod(o1::PyObject, o2::PyObject) = pymod(PyObject, o1, o2)
Base.:(^)(o1::PyObject, o2::PyObject) = pypow(PyObject, o1, o2)
Base.:(<<)(o1::PyObject, o2::PyObject) = pylshift(PyObject, o1, o2)
Base.:(>>)(o1::PyObject, o2::PyObject) = pyrshift(PyObject, o1, o2)
Base.:(&)(o1::PyObject, o2::PyObject) = pyand(PyObject, o1, o2)
Base.:xor(o1::PyObject, o2::PyObject) = pyxor(PyObject, o1, o2)
Base.:(|)(o1::PyObject, o2::PyObject) = pyor(PyObject, o1, o2)

Base.:(+)(o1::PyObject, o2::Number) = pyadd(PyObject, o1, o2)
Base.:(-)(o1::PyObject, o2::Number) = pysub(PyObject, o1, o2)
Base.:(*)(o1::PyObject, o2::Number) = pymul(PyObject, o1, o2)
Base.:(/)(o1::PyObject, o2::Number) = pytruediv(PyObject, o1, o2)
Base.:fld(o1::PyObject, o2::Number) = pyfloordiv(PyObject, o1, o2)
Base.:mod(o1::PyObject, o2::Number) = pymod(PyObject, o1, o2)
Base.:(^)(o1::PyObject, o2::Number) = pypow(PyObject, o1, o2)
Base.:(<<)(o1::PyObject, o2::Number) = pylshift(PyObject, o1, o2)
Base.:(>>)(o1::PyObject, o2::Number) = pyrshift(PyObject, o1, o2)
Base.:(&)(o1::PyObject, o2::Number) = pyand(PyObject, o1, o2)
Base.:xor(o1::PyObject, o2::Number) = pyxor(PyObject, o1, o2)
Base.:(|)(o1::PyObject, o2::Number) = pyor(PyObject, o1, o2)

Base.:(+)(o1::Number, o2::PyObject) = pyadd(PyObject, o1, o2) # Defining +(::Any, ::PyObject) like this hangs Julia v1.5.2-v1.5.3 (at least) during precompilation
Base.:(-)(o1::Number, o2::PyObject) = pysub(PyObject, o1, o2)
Base.:(*)(o1::Number, o2::PyObject) = pymul(PyObject, o1, o2)
Base.:(/)(o1::Number, o2::PyObject) = pytruediv(PyObject, o1, o2)
Base.:fld(o1::Number, o2::PyObject) = pyfloordiv(PyObject, o1, o2)
Base.:mod(o1::Number, o2::PyObject) = pymod(PyObject, o1, o2)
# Base.:(^)(o1::Number, o2::PyObject) = pypow(PyObject, o1, o2)
# Base.:(<<)(o1::Number, o2::PyObject) = pylshift(PyObject, o1, o2)
# Base.:(>>)(o1::Number, o2::PyObject) = pyrshift(PyObject, o1, o2)
Base.:(&)(o1::Number, o2::PyObject) = pyand(PyObject, o1, o2)
Base.:xor(o1::Number, o2::PyObject) = pyxor(PyObject, o1, o2)
Base.:(|)(o1::Number, o2::PyObject) = pyor(PyObject, o1, o2)

# ternary
Base.powermod(o1::PyObject, o2::Union{PyObject,Number}, o3::Union{PyObject,Number}) =
    pypow(PyObject, o1, o2, o3)

### DOCUMENTATION

# function Base.Docs.getdoc(o::PyObject)
#     function tryget(f, g=identity)
#         a = try
#             f()
#         catch
#             return nothing
#         end
#         pyisnone(a) ? nothing : g(a)
#     end
#     # function name(o)
#     #     n = tryget(()->o.__name__)
#     #     n === nothing && return nothing
#     #     m = tryget(()->o.__module__)
#     #     (m === nothing || string(m) == "builtins") ? string(n) : string(m, ".", n)
#     # end
#     getname(o) = tryget(()->o.__name__, string)
#     getdoc(o) = tryget(()->o.__doc__, string ∘ ins.cleandoc)

#     docs = []

#     # Short description
#     ins = pyimport("inspect")
#     desc, name = if ins.ismodule(o).jl!b
#         (pyhasattr(o, "__path__") ? "package" : "module"), gname(o)
#     elseif ins.isgetsetdescriptor(o).jl!b
#         ("getset descriptor", "$(getname(o.__objclass__)).$(o.__name__)")
#     elseif ins.ismemberdescriptor(o).jl!b
#         ("member descriptor", "$(getname(o.__objclass__)).$(o.__name__)")
#     elseif ins.isclass(o).jl!b
#         ("class", getname(o))
#     elseif ins.isfunction(o).jl!b || ins.isbuiltin(o).jl!b
#         ("function", getname(o))
#     elseif ins.ismethod(o).jl!b
#         ("method", getname(o))
#     else
#         ("object of type", getname(pytype(o)))
#     end
#     push!(docs, Markdown.Paragraph(["Python ", desc, " ", Markdown.Code(name), "."]))

#     if ins.isroutine(o).jl!b || ins.isclass(o).jl!b
#         try
#             push!(docs, Markdown.Code("python", "$(o.__name__)$(ins.signature(o))"))
#         catch
#         end
#     end

#     # Maybe document the class instead
#     doc = getdoc(o)
#     if doc === nothing && !ins.ismodule(o).jl!b && !ins.isclass(o).jl!b && !ins.isroutine(o).jl!b && !ins.isdatadescriptor(o).jl!b
#         o = pyhasattr(o, "__origin__") ? o.__origin__ : pytype(o)
#         doc = getdoc(o)
#     end
#     doc === nothing || push!(docs, Markdown.Paragraph([Markdown.Text(doc)]))

#     # Done
#     Markdown.MD(docs)
# end
# Base.Docs.Binding(o::PyObject, k::Symbol) = getproperty(o, k)


text:
using RvSpectraKitLearn
using MultivariateStats

# Likely need to update this to point to directory containing the spectra you wish to analyze
#include(joinpath(Pkg.dir("RvSpectraKitLearn"),"path_to_spectra.jl"))
path_to_spectra = joinpath(Pkg.dir(),"..","..","Code","XavierSpotData")
cd(path_to_spectra)

function get_filenames_multi_lats(; num_max_files::Integer = 0)
  files = readdir()
  datafiles = String[]
  strengths = Float64[]
  inclinations = Float64[]
  latitudes = Float64[]
  resolutions = Float64[]
  samplings = Float64[]
  phases = Float64[]
  for f in files
    m = match(r"(\w+)_(\d+)pc_incl(\d+)_lat(\d+)_reso(\d+)k_(-?\d+\.?\d*)pxsampl_(-?\d+\.?\d*)\.csv",f)
    if m==nothing continue end
    effect_name = m.captures[1]
    resolution = parse(m.captures[5])
    if resolution> 149
    push!(strengths,parse(m.captures[2]))
	  push!(inclinations,parse(m.captures[3]))
    push!(latitudes,parse(m.captures[4]))
    push!(resolutions,parse(m.captures[5]))
    push!(samplings,parse(m.captures[6]))
    push!(phases,parse(m.captures[7]))
    push!(datafiles,f)
    end
    if num_max_files>0 && length(datafiles)==num_max_files break end
  end
  @assert(length(datafiles)>0)
  return (datafiles,phases)
  #idx = sortperm(phases)
  #return (datafiles[idx],phases[idx])
end

(datafiles, phases) = get_filenames_multi_lats(); # num_max_files=4);
(lambda, obs_noisefree) = read_filelist(datafiles);

datafiles

sampling = 3.0
planet_amplitude = 10.0
rvs_true = 0.0*cos(2pi*phases)

snr = 200.0
obs = make_noisy_spectra(obs_noisefree,snr,sampling=sampling)

pca_out = fit_pca_default(obs);
pca_scores = transform(pca_out,obs)   # compute scores using Julia's PCA routine for comparison's sake
pca_eford_out = fit_pca_eford(obs);   # compute first few scores using Eric's itterative PCA routine

doppler_comp_simple = calc_doppler_component_simple(lambda,obs);
genpca_simple_out = fit_gen_pca_rv(obs,doppler_comp_simple)

rvs_out_simple = est_rvs_from_pc(obs,genpca_simple_out[1],genpca_simple_out[2][:,1])

genpca_simple_out_1 = genpca_simple_out[3]
genpca_simple_out_2 = genpca_simple_out[3]
plot(rvs_out_simple,genpca_simple_out[3][:,2],"r.")
plot(rvs_out_simple,genpca_simple_out[3][:,3],"b.")
plot(rvs_out_simple,genpca_simple_out[3][:,4],"g.")
 xlabel("Apparent Velocity (m/s)")
 ylabel("Next PC Coefficients")

plot(rvs_out_simple,genpca_simple_out[3][:,2],"r.")
lat10_mask = find(x->ismatch(r"lat10",x),datafiles)
lat40_mask = find(x->ismatch(r"lat40",x),datafiles)
lat70_mask = find(x->ismatch(r"lat70",x),datafiles)
plot(rvs_out_simple[lat70_mask],genpca_simple_out[3][lat70_mask,2],"b.")
plot(rvs_out_simple[lat40_mask],genpca_simple_out[3][lat40_mask,2],"g.")
plot(rvs_out_simple[lat10_mask],genpca_simple_out[3][lat10_mask,2],"m.")
 xlabel("Apparent Velocity (m/s)")
 ylabel("Next PC Coefficient")

abs(genpca_simple_out_1.-genpca_simple_out_2)

doppler_comp_gp = calc_doppler_component_gp(lambda,obs);
genpca_gp_out = fit_gen_pca_rv(obs,doppler_comp_gp)
rvs_out_gp = est_rvs_from_pc(obs,genpca_gp_out[1],genpca_gp_out[2][:,1])
rms_simple = sqrt(sumabs2(rvs_out_simple.-rvs_true)/length(rvs_true))
rms_gp = sqrt(sumabs2(rvs_out_gp.-rvs_true)/length(rvs_true))
println("snr = ", snr, " rms = ", rms_simple, "  ", rms_gp)

plot(rvs_out_simple,genpca_simple_out[3][:,3],"r.")
 plot(rvs_out_gp,genpca_gp_out[3][:,3],"b.")

using PyPlot
plot(phases,rvs_true,"b-")
 plot(phases,rvs_out_simple,"g.")
 xlabel("Phase")
 ylabel("Apparent Velocity (m/s)")

plot(phases,rvs_out_gp,"g.")
 plot(phases,rvs_out_simple.-rvs_true,"r.")
 plot(phases,rvs_out_gp.-rvs_true,"g.")

ext: jl

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
type SwitchChoiceModel <: ChoiceModel
	casecms::Vector{Tuple{Vector{UInt}, ChoiceModel}}
	defaultcm::ChoiceModel
	function SwitchChoiceModel(casecms::Vector, defaultcm::ChoiceModel)
		new(casecms, defaultcm)
	end
end

function godelnumber(cm::SwitchChoiceModel, cc::ChoiceContext)
	for (cpids, casecm) in cm.casecms
		if cc.cpid in cpids
			return godelnumber(casecm, cc)
		end
	end
	godelnumber(cm.defaultcm, cc)
end

function resetstate!(cm::SwitchChoiceModel)
	foreach(c->resetstate!(c[2]), cm.casecms)
	resetstate!(cm.defaultcm)
end

getparams(cm::SwitchChoiceModel) = [vcat(map(c->getparams(c[2]), cm.casecms)...); getparams(cm.defaultcm)]

function setparams!(cm::SwitchChoiceModel, params)
	if length(params) != numparams(cm)
		error("expected $(numparams(cm)) model parameter(s), but got $(length(params))")
	end
	startidx = 1
	for (cpids, casecm) in cm.casecms
		endidx = startidx + numparams(casecm) - 1
		if endidx >= startidx
			setparams!(casecm, params[startidx:endidx])
		end
		startidx = endidx + 1
	end
	endidx = startidx + numparams(cm.defaultcm) - 1
	if endidx >= startidx
		setparams!(cm.defaultcm, params[startidx:endidx])
	end
end

paramranges(cm::SwitchChoiceModel) = [vcat(map(c->paramranges(c[2]), cm.casecms)...); paramranges(cm.defaultcm)]

show(io::IO, cm::SwitchChoiceModel) = print(io, "Switch choice model with cases: ", map(c->c[2], cm.casecms)..., "; default: ", cm.defaultcm)


text:
function overlap_ellipse(	xAB::Float64,
						yAB::Float64,
						zAB::Float64,
						a11::Float64,
						a12::Float64,
						a13::Float64,
						a21::Float64,
						a22::Float64,
						a23::Float64,
						a31::Float64,
						a32::Float64,
						a33::Float64,
						b11::Float64,
						b12::Float64,
						b13::Float64,
						b21::Float64,
						b22::Float64,
						b23::Float64,
						b31::Float64,
						b32::Float64,
						b33::Float64)

	number_of_iterations_overlap_criterion::Int64 = 20

	a::Float64 = 0.0
	b::Float64 = 0.0
	c::Float64 = 0.0
	x::Float64 = 0.0
	fb::Float64 = 0.0
	fx::Float64 = 0.0
	overlapfun::Float64 = 0.0
	iter::Int64 = 0

	const4t::Float64 = (a12 * a21 - a11 * a22 + a11 * b22 - a12 * b21 - a21 * b12 + a22 * b11 - b11 * b22 + b12 * b21) * zAB^2 + (yAB * (a11 * a23 - a13 * a21 + a11 * a32 - a12 * a31 - a11 * b23 + a13 * b21 + a21 * b13 - a23 * b11 - a11 * b32 + a12 * b31 + a31 * b12 - a32 * b11 + b11 * b23 - b13 * b21 + b11 * b32 - b12 * b31) - xAB * (a12 * a23 - a13 * a22 + a21 * a32 - a22 * a31 - a12 * b23 + a13 * b22 + a22 * b13 - a23 * b12 - a21 * b32 + a22 * b31 + a31 * b22 - a32 * b21 + b12 * b23 - b13 * b22 + b21 * b32 - b22 * b31)) * zAB + (a23 * a32 - a22 * a33 + a22 * b33 - a23 * b32 - a32 * b23 + a33 * b22 - b22 * b33 + b23 * b32) * xAB^2 + (a12 * a33 - a13 * a32 + a21 * a33 - a23 * a31 - a12 * b33 + a13 * b32 + a32 * b13 - a33 * b12 - a21 * b33 + a23 * b31 + a31 * b23 - a33 * b21 + b12 * b33 - b13 * b32 + b21 * b33 - b23 * b31) * xAB * yAB + (a13 * a31 - a11 * a33 + a11 * b33 - a13 * b31 - a31 * b13 + a33 * b11 - b11 * b33 + b13 * b31) * yAB^2
	const3t::Float64 = (3.0 * a11 * a22 - 3.0 * a12 * a21 - 2.0 * a11 * b22 + 2.0 * a12 * b21 + 2.0 * a21 * b12 - 2.0 * a22 * b11 + b11 * b22 - b12 * b21) * zAB^2 + (xAB * (3.0 * a12 * a23 - 3.0 * a13 * a22 + 3.0 * a21 * a32 - 3.0 * a22 * a31 - 2.0 * a12 * b23 + 2.0 * a13 * b22 + 2.0 * a22 * b13 - 2.0 * a23 * b12 - 2.0 * a21 * b32 + 2.0 * a22 * b31 + 2.0 * a31 * b22 - 2.0 * a32 * b21 + b12 * b23 - b13 * b22 + b21 * b32 - b22 * b31) - yAB * (3.0 * a11 * a23 - 3.0 * a13 * a21 + 3.0 * a11 * a32 - 3.0 * a12 * a31 - 2.0 * a11 * b23 + 2.0 * a13 * b21 + 2.0 * a21 * b13 - 2.0 * a23 * b11 - 2.0 * a11 * b32 + 2.0 * a12 * b31 + 2.0 * a31 * b12 - 2.0 * a32 * b11 + b11 * b23 - b13 * b21 + b11 * b32 - b12 * b31)) * zAB + (3.0 * a22 * a33 - 3.0 * a23 * a32 - 2.0 * a22 * b33 + 2.0 * a23 * b32 + 2.0 * a32 * b23 - 2.0 * a33 * b22 + b22 * b33 - b23 * b32) * xAB^2 + (3.0 * a13 * a32 - 3.0 * a12 * a33 - 3.0 * a21 * a33 + 3.0 * a23 * a31 + 2.0 * a12 * b33 - 2.0 * a13 * b32 - 2.0 * a32 * b13 + 2.0 * a33 * b12 + 2.0 * a21 * b33 - 2.0 * a23 * b31 - 2.0 * a31 * b23 + 2.0 * a33 * b21 - b12 * b33 + b13 * b32 - b21 * b33 + b23 * b31) * xAB * yAB + (3.0 * a11 * a33 - 3.0 * a13 * a31 - 2.0 * a11 * b33 + 2.0 * a13 * b31 + 2.0 * a31 * b13 - 2.0 * a33 * b11 + b11 * b33 - b13 * b31) * yAB^2
	const2t::Float64 = (3.0 * a12 * a21 - 3.0 * a11 * a22 + a11 * b22 - a12 * b21 - a21 * b12 + a22 * b11) * zAB^2 + (yAB * (3.0 * a11 * a23 - 3.0 * a13 * a21 + 3.0 * a11 * a32 - 3.0 * a12 * a31 - a11 * b23 + a13 * b21 + a21 * b13 - a23 * b11 - a11 * b32 + a12 * b31 + a31 * b12 - a32 * b11) - xAB * (3.0 * a12 * a23 - 3.0 * a13 * a22 + 3.0 * a21 * a32 - 3.0 * a22 * a31 - a12 * b23 + a13 * b22 + a22 * b13 - a23 * b12 - a21 * b32 + a22 * b31 + a31 * b22 - a32 * b21)) * zAB + (3.0 * a23 * a32 - 3.0 * a22 * a33 + a22 * b33 - a23 * b32 - a32 * b23 + a33 * b22) * xAB^2 + (3.0 * a12 * a33 - 3.0 * a13 * a32 + 3.0 * a21 * a33 - 3.0 * a23 * a31 - a12 * b33 + a13 * b32 + a32 * b13 - a33 * b12 - a21 * b33 + a23 * b31 + a31 * b23 - a33 * b21) * xAB * yAB + (3.0 * a13 * a31 - 3.0 * a11 * a33 + a11 * b33 - a13 * b31 - a31 * b13 + a33 * b11) * yAB^2
	const1t::Float64 = (a11 * a22 - a12 * a21) * zAB^2 + (xAB * (a12 * a23 - a13 * a22 + a21 * a32 - a22 * a31) - yAB * (a11 * a23 - a13 * a21 + a11 * a32 - a12 * a31)) * zAB + (a22 * a33 - a23 * a32) * xAB^2 + (a13 * a32 - a12 * a33 - a21 * a33 + a23 * a31) * xAB * yAB + (a11 * a33 - a13 * a31) * yAB^2
	const3n::Float64 = a11 * a23 * a32 - a11 * a22 * a33 + a12 * a21 * a33 - a12 * a23 * a31 - a13 * a21 * a32 + a13 * a22 * a31 + a11 * a22 * b33 - a11 * a23 * b32 - a11 * a32 * b23 + a11 * a33 * b22 - a12 * a21 * b33 + a12 * a23 * b31 + a12 * a31 * b23 - a12 * a33 * b21 + a13 * a21 * b32 - a13 * a22 * b31 - a13 * a31 * b22 + a13 * a32 * b21 + a21 * a32 * b13 - a21 * a33 * b12 - a22 * a31 * b13 + a22 * a33 * b11 + a23 * a31 * b12 - a23 * a32 * b11 - a11 * b22 * b33 + a11 * b23 * b32 + a12 * b21 * b33 - a12 * b23 * b31 - a13 * b21 * b32 + a13 * b22 * b31 + a21 * b12 * b33 - a21 * b13 * b32 - a22 * b11 * b33 + a22 * b13 * b31 + a23 * b11 * b32 - a23 * b12 * b31 - a31 * b12 * b23 + a31 * b13 * b22 + a32 * b11 * b23 - a32 * b13 * b21 - a33 * b11 * b22 + a33 * b12 * b21 + b11 * b22 * b33 - b11 * b23 * b32 - b12 * b21 * b33 + b12 * b23 * b31 + b13 * b21 * b32 - b13 * b22 * b31
	const2n::Float64 = 3.0 * a11 * a22 * a33 - 3.0 * a11 * a23 * a32 - 3.0 * a12 * a21 * a33 + 3.0 * a12 * a23 * a31 + 3.0 * a13 * a21 * a32 - 3.0 * a13 * a22 * a31 - 2.0 * a11 * a22 * b33 + 2.0 * a11 * a23 * b32 + 2.0 * a11 * a32 * b23 - 2.0 * a11 * a33 * b22 + 2.0 * a12 * a21 * b33 - 2.0 * a12 * a23 * b31 - 2.0 * a12 * a31 * b23 + 2.0 * a12 * a33 * b21 - 2.0 * a13 * a21 * b32 + 2.0 * a13 * a22 * b31 + 2.0 * a13 * a31 * b22 - 2.0 * a13 * a32 * b21 - 2.0 * a21 * a32 * b13 + 2.0 * a21 * a33 * b12 + 2.0 * a22 * a31 * b13 - 2.0 * a22 * a33 * b11 - 2.0 * a23 * a31 * b12 + 2.0 * a23 * a32 * b11 + a11 * b22 * b33 - a11 * b23 * b32 - a12 * b21 * b33 + a12 * b23 * b31 + a13 * b21 * b32 - a13 * b22 * b31 - a21 * b12 * b33 + a21 * b13 * b32 + a22 * b11 * b33 - a22 * b13 * b31 - a23 * b11 * b32 + a23 * b12 * b31 + a31 * b12 * b23 - a31 * b13 * b22 - a32 * b11 * b23 + a32 * b13 * b21 + a33 * b11 * b22 - a33 * b12 * b21
	const1n::Float64 = 3.0 * a11 * a23 * a32 - 3.0 * a11 * a22 * a33 + 3.0 * a12 * a21 * a33 - 3.0 * a12 * a23 * a31 - 3.0 * a13 * a21 * a32 + 3.0 * a13 * a22 * a31 + a11 * a22 * b33 - a11 * a23 * b32 - a11 * a32 * b23 + a11 * a33 * b22 - a12 * a21 * b33 + a12 * a23 * b31 + a12 * a31 * b23 - a12 * a33 * b21 + a13 * a21 * b32 - a13 * a22 * b31 - a13 * a31 * b22 + a13 * a32 * b21 + a21 * a32 * b13 - a21 * a33 * b12 - a22 * a31 * b13 + a22 * a33 * b11 + a23 * a31 * b12 - a23 * a32 * b11

	b = 0.5

	fb = (((const4t * b + const3t) * b + const2t) * b + const1t) * b/(((const3n * b + const2n) * b + const1n) * b)
	if fb >= 1.0
		overlapfun = 1.0
	else
		iter = 1
		a = 0.0
		c = 1.0
		while fb < 1 && iter < number_of_iterations_overlap_criterion
			# Pick point in upper half.
			x = 0.5 * (b + c)
			fx = (((const4t * x + const3t) * x + const2t) * x + const1t) * x/(((const3n * x + const2n) * x + const1n) * x)
			if fx < fb
				c = x
			else
				b = x
				fb = fx
			end
			# Pick point in lower half.
			x = 0.5 * (a + b)
			fx = (((const4t * x + const3t) * x + const2t) * x + const1t) * x/(((const3n * x + const2n) * x + const1n) * x)

			if fx < fb
				a = x
			else
				b = x
				fb = fx
			end

			iter = iter + 1
		end
		overlapfun = fb
	end

	return overlapfun
end

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
# .................................. save to file for FITS array ...................................................

function _fits_save(FITS)

    o = IOBuffer()

    for i ∈ eachindex(FITS)
        a = _write_header(FITS[i])
        Base.write(o,Array{UInt8,1}(a.data))
        b = _write_data(FITS[i])
        b.size > 0 && Base.write(o,Array{UInt8,1}(b.data))
    end

    return _fits_write_IO(o, FITS[1].filename)                    # same filename in all HDUs

end

function _fits_write_IO(o::IO, filename::String)

    s = Base.open(filename,"w")
        Base.write(s,o.data)
        Base.close(s)

end

# .................................. write-io (header sector) .......................................................

function _write_header(FITS_HDU)

    o = IOBuffer()

    Base.seekstart(o)

    records = FITS_HDU.header.records
      nrecs = Base.length(records)
     blanks = [Base.repeat(' ',80) for i=1:(36 - nrecs % 36)]      # complement last header block with blank records
    records = Base.append!(records,blanks)
     nbytes = Base.write(o,Array{UInt8,1}(join(records)))

    return o

end

# .................................. write-io for FITS_HDU (data sector) ..........................................

function _write_data(FITS_HDU)

    hdutype = FITS_HDU.dataobject.hdutype

    hdutype == "PRIMARY"  && return _write_IMAGE_data(FITS_HDU)
    hdutype == "IMAGE"    && return _write_IMAGE_data(FITS_HDU)
    hdutype == "TABLE"    && return _write_TABLE_data(FITS_HDU)
    hdutype == "BINTABLE" && return _write_BINTABLE_data(FITS_HDU)

    return error("FitsError: '$hdutype': not a 'FITS standard extension'")

end

function _write_IMAGE_data(FITS_HDU)

    o = IOBuffer()

    Base.seekstart(o)

     data = FITS_HDU.dataobject.data
    ndata = Base.length(data)

    ndata == 0 && return o

         E = Base.eltype(data)
         E <: Real || error("FitsError: incorrect DataType (Real type mandatory for image HDUs)")
    nbytes = E.size
     bzero = _fits_bzero(E)
      data = Base.vec(data)
      data = data .- E(bzero)                                           # change between Int and UInt (if applicable)
      data = hton.(data)                                                # change from 'host' to 'network' ordering

    [Base.write(o,data[i]) for i=1:ndata]                               # write data
    [Base.write(o,E(0)) for i=1:((2880÷nbytes)-ndata % (2880÷nbytes))]  # complement with type E zero elements

    return o

end

function _write_TABLE_data(FITS_HDU)

    o = IOBuffer()

    Base.seekstart(o)

    records = join(FITS_HDU.dataobject.data)                      # Array of ASCII records
      nrecs = Base.length(records)                                # number of ASCII records
      lrecs = Base.length(records[1])                             # length of ASCII records
      nchar = 2880 - (nrecs * lrecs) % 2880                       # number of blanks to complement last data block
     blanks = Base.repeat(' ',nchar)                              # complement last data block with blanks
     nbytes = Base.write(o,Array{UInt8,1}(records * blanks))

    return o

end


text:
export EI_CLASS, EI_DATA, EI_VERSION, EI_OSABI, ELFMAG1, ELFMAG2, ELFMAG3, ELFMAG4, elf64_st_type, elf64_r_info, elf64_r_sym

include("./types.jl")

macro enum_export(T, syms...)
    if T isa Symbol
        Ttype = typeof(T)
        esc(quote
            @enum $T $(syms...)
            export $T
            $([:(export $sym) for sym in syms]...)
        end)
    elseif (T.head === :(::) && T isa Expr && length(T.args) == 2 && T.args[1] isa Symbol)
        if eltype(syms) == Expr && length(syms) == 1
            vals = filter(x -> x isa Expr, syms[1].args)
            esc(quote
                @enum $(T.args[1])::$(T.args[2]) $(syms...)
                export $(T.args[1])
                $([:(export $(sym.args[1])) for sym in vals]...)
            end)
        else
            esc(quote
                @enum $(T.args[1])::$(T.args[2]) $(syms...)
                export $(T.args[1])
                $([:(export $sym) for sym in syms]...)
            end)
        end
    end
end

const EI_CLASS = 5
@enum_export EHClass::UInt8 begin
    ELFCLASSNONE = 0
    ELFCLASS32 = 1
    ELFCLASS64 = 2
end

const EI_DATA = 6
@enum_export EHEndian::UInt8 begin
    ELFDATANONE = 0
    ELFDATA2MSB = 1
    ELFDATA2LSB = 2
end

const EI_VERSION = 7
@enum_export EHVersion::UInt8 begin
    EV_NONE = 0
    EV_CURRENT = 1
end

const EI_OSABI = 8
@enum_export EHOsAbi::UInt8 begin
    ELFOSABI_SYSV = 0
    ELFOSABI_HPUX = 1
    ELFOSABI_NETBSD = 2
    ELFOSABI_LINUX = 3
    ELFOSABI_SOLARIS = 6
    ELFOSABI_AIX = 7
    ELFOSABI_IRIX = 8
    ELFOSABI_FREEBSD = 9
    ELFOSABI_TRU64 = 10
    ELFOSABI_MODESTO = 11
    ELFOSABI_OPENBSD = 12
    ELFOSABI_ARM_AEABI = 64
    ELFOSABI_ARM = 97
    ELFOSABI_STANDALONE = 255
end

@enum_export EHType::Elf64_Half begin
    ET_NONE = 0
    ET_REL = 1
    ET_EXEC = 2
    ET_DYN = 3
    ET_CORE = 4
    ET_NUM = 5
    ET_LOOS = 0xfe00
    ET_HIOS = 0xfeff
    ET_LOPROC = 0xff00
    ET_HIPROC = 0xffff
end

# for e_machine in Ehdr
@enum_export EHMachine::Elf64_Half begin
    EM_NONE = 0
    EM_M32 = 1
    EM_SPARC = 2
    EM_386 = 3
    EM_68K = 4
    EM_88K = 5
    EM_IAMCU = 6
    EM_860 = 7
    EM_MIPS = 8
    EM_S370 = 9
    EM_MIPS_RS3_LE = 10
    EM_PARISC = 15
    EM_VPP500 = 17
    EM_SPARC32PLUS = 18
    EM_960 = 19
    EM_PPC = 20
    EM_PPC64 = 21
    EM_S390 = 22
    EM_SPU = 23
    EM_V800 = 36
    EM_FR20 = 37
    EM_RH32 = 38
    EM_RCE = 39
    EM_ARM = 40
    EM_FAKE_ALPHA = 41
    EM_SH = 42
    EM_SPARCV9 = 43
    EM_TRICORE = 44
    EM_ARC = 45
    EM_H8_300 = 46
    EM_H8_300H = 47
    EM_H8S = 48
    EM_H8_500 = 49
    EM_IA_64 = 50
    EM_MIPS_X = 51
    EM_COLDFIRE = 52
    EM_68HC12 = 53
    EM_MMA = 54
    EM_PCP = 55
    EM_NCPU = 56
    EM_NDR1 = 57
    EM_STARCORE = 58
    EM_ME16 = 59
    EM_ST100 = 60
    EM_TINYJ = 61
    EM_X86_64 = 62
    EM_PDSP = 63
    EM_PDP10 = 64
    EM_PDP11 = 65
    EM_FX66 = 66
    EM_ST9PLUS = 67
    EM_ST7 = 68
    EM_68HC16 = 69
    EM_68HC11 = 70
    EM_68HC08 = 71
    EM_68HC05 = 72
    EM_SVX = 73
    EM_ST19 = 74
    EM_VAX = 75
    EM_CRIS = 76
    EM_JAVELIN = 77
    EM_FIREPATH = 78
    EM_ZSP = 79
    EM_MMIX = 80
    EM_HUANY = 81
    EM_PRISM = 82
    EM_AVR = 83
    EM_FR30 = 84
    EM_D10V = 85
    EM_D30V = 86
    EM_V850 = 87
    EM_M32R = 88
    EM_MN10300 = 89
    EM_MN10200 = 90
    EM_PJ = 91
    EM_OPENRISC = 92
    EM_ARC_COMPACT = 93
    EM_XTENSA = 94
    EM_VIDEOCORE = 95
    EM_TMM_GPP = 96
    EM_NS32K = 97
    EM_TPC = 98
    EM_SNP1K = 99
    EM_ST200 = 100
    EM_IP2K = 101
    EM_MAX = 102
    EM_CR = 103
    EM_F2MC16 = 104
    EM_MSP430 = 105
    EM_BLACKFIN = 106
    EM_SE_C33 = 107
    EM_SEP = 108
    EM_ARCA = 109
    EM_UNICORE = 110
    EM_EXCESS = 111
    EM_DXP = 112
    EM_ALTERA_NIOS2 = 113
    EM_CRX = 114
    EM_XGATE = 115
    EM_C166 = 116
    EM_M16C = 117
    EM_DSPIC30F = 118
    EM_CE = 119
    EM_M32C = 120
    EM_TSK3000 = 131
    EM_RS08 = 132
    EM_SHARC = 133
    EM_ECOG2 = 134
    EM_SCORE7 = 135
    EM_DSP24 = 136
    EM_VIDEOCORE3 = 137
    EM_LATTICEMICO32 = 138
    EM_SE_C17 = 139
    EM_TI_C6000 = 140
    EM_TI_C2000 = 141
    EM_TI_C5500 = 142
    EM_TI_ARP32 = 143
    EM_TI_PRU = 144
    EM_MMDSP_PLUS = 160
    EM_CYPRESS_M8C = 161
    EM_R32C = 162
    EM_TRIMEDIA = 163
    EM_QDSP6 = 164
    EM_8051 = 165
    EM_STXP7X = 166
    EM_NDS32 = 167
    EM_ECOG1X = 168
    EM_MAXQ30 = 169
    EM_XIMO16 = 170
    EM_MANIK = 171
    EM_CRAYNV2 = 172
    EM_RX = 173
    EM_METAG = 174
    EM_MCST_ELBRUS = 175
    EM_ECOG16 = 176
    EM_CR16 = 177
    EM_ETPU = 178
    EM_SLE9X = 179
    EM_L10M = 180
    EM_K10M = 181
    EM_AARCH64 = 183
    EM_AVR32 = 185
    EM_STM8 = 186
    EM_TILE64 = 187
    EM_TILEPRO = 188
    EM_MICROBLAZE = 189
    EM_CUDA = 190
    EM_TILEGX = 191
    EM_CLOUDSHIELD = 192
    EM_COREA_1ST = 193
    EM_COREA_2ND = 194
    EM_ARC_COMPACT2 = 195
    EM_OPEN8 = 196
    EM_RL78 = 197
    EM_VIDEOCORE5 = 198
    EM_78KOR = 199
    EM_56800EX = 200
    EM_BA1 = 201
    EM_BA2 = 202
    EM_XCORE = 203
    EM_MCHP_PIC = 204
    EM_KM32 = 210
    EM_KMX32 = 211
    EM_EMX16 = 212
    EM_EMX8 = 213
    EM_KVARC = 214
    EM_CDP = 215
    EM_COGE = 216
    EM_COOL = 217
    EM_NORC = 218
    EM_CSR_KALIMBA = 219
    EM_Z80 = 220
    EM_VISIUM = 221
    EM_FT32 = 222
    EM_MOXIE = 223
    EM_AMDGPU = 224
    EM_RISCV = 243
    EM_BPF = 247
    EM_CSKY = 252
end

const EI_NIDENT = 16
const EIMAG1 = 1
const ELFMAG1 = 0x7f
const EIMAG2 = 2
const ELFMAG2 = UInt8('E')
const EIMAG3 = 3
const ELFMAG3 = UInt8('L')
const EIMAG4 = 4
const ELFMAG4 = UInt8('F')

# for sh_name in Shdr
@enum_export SHName::UInt32 begin
    SHN_UNDEF = 0
    SHN_LORESERVE = 0xff00
    # SHN_LOPROC = 0xff00
    # SHN_BEFORE = 0xff00
    SHN_AFTER = 0xff01
    SHN_HIPROC = 0xff1f
    SHN_LOOS = 0xff20
    SHN_HIOS = 0xff3f
    SHN_ABS = 0xfff1
    SHN_COMMON = 0xfff2
    SHN_XINDEX = 0xffff
    # SHN_HIRESERVE = 0xffff
end

# for sh_type in Shdr
@enum_export SHType::UInt32 begin
    SHT_NULL = 0
    SHT_PROGBITS = 1
    SHT_SYMTAB = 2
    SHT_STRTAB = 3
    SHT_RELA = 4
    SHT_HASH = 5
    SHT_DYNAMIC = 6
    SHT_NOTE = 7
    SHT_NOBITS = 8
    SHT_REL = 9
    SHT_SHLIB = 10
    SHT_DYNSYM = 11
    SHT_INIT_ARRAY = 14
    SHT_FINI_ARRAY = 15
    SHT_PREINIT_ARRAY = 16
    SHT_GROUP = 17
    SHT_SYMTAB_SHNDX = 18
    SHT_NUM = 19
    SHT_LOOS = 0x60000000
    SHT_GNU_ATTRIBUTES = 0x6ffffff5
    SHT_GNU_HASH = 0x6ffffff6
    SHT_GNU_LIBLIST = 0x6ffffff7
    SHT_CHECKSUM = 0x6ffffff8
    SHT_LOSUNW = 0x6ffffffa
    # SHT_SUNW_move = 0x6ffffffa
    SHT_SUNW_COMDAT = 0x6ffffffb
    SHT_SUNW_syminfo = 0x6ffffffc
    SHT_GNU_verdef = 0x6ffffffd
    SHT_GNU_verneed = 0x6ffffffe
    SHT_GNU_versym = 0x6fffffff
    # SHT_HISUNW = 0x6fffffff
    # SHT_HIOS = 0x6fffffff
    SHT_LOPROC = 0x70000000
    SHT_HIPROC = 0x7fffffff
    SHT_LOUSER = 0x80000000
    SHT_HIUSER = 0x8fffffff
end

# for sh_flags in Shdr
@enum_export SHFlags::UInt32 begin
    SHF_WRITE = (1 << 0)
    SHF_ALLOC = (1 << 1)
    SHF_EXECINSTR = (1 << 2)
    SHF_MERGE = (1 << 4)
    SHF_STRINGS = (1 << 5)
    SHF_INFO_LINK = (1 << 6)
    SHF_LINK_ORDER = (1 << 7)
    SHF_OS_NONCONFORMING = (1 << 8)
    SHF_GROUP = (1 << 9)
    SHF_TLS = (1 << 10)
    SHF_COMPRESSED = (1 << 11)
    SHF_MASKOS = 0x0ff00000
    SHF_MASKPROC = 0xf0000000
    SHF_ORDERED = (1 << 30)
    SHF_EXCLUDE = (1 << 31)
end

# if sh_flags = SHF_COMPRESSED
@enum_export CHType::UInt32 begin
    ELFCOMPRESS_ZLIB = 1
    ELFCOMPRESS_LOOS = 0x60000000
    ELFCOMPRESS_HIOS = 0x6fffffff
    ELFCOMPRESS_LOPROC = 0x70000000
    ELFCOMPRESS_HIPROC = 0x7fffffff
end

# for si_boundto
@enum_export SymInfoBT::UInt16 begin
    SYMINFO_BT_SELF = 0xffff
    SYMINFO_BT_PARENT = 0xfffe
    SYMINFO_BT_LOWRESERVE = 0xff00
end

# for si_flags
@enum_export SymInfoFlags::UInt16 begin
    SYMINFO_FLG_DIRECT = 0x0001
    SYMINFO_FLG_PASSTHRU = 0x0002
    SYMINFO_FLG_COPY = 0x0004
    SYMINFO_FLG_LAZYLOAD = 0x0008
end

@enum_export SymInfoVersion::UInt8 begin
    SYMINFO_NONE
    SYMINFO_CURRENT
end

# for st_info
@enum_export STBind::UInt8 begin
    STB_LOCAL = 0 #  Local symbol
    STB_GLOBAL = 1 #  Global symbol
    STB_WEAK = 2 #  Weak symbol
    STB_NUM = 3 #  Number of defined types.
    # STB_LOOS = 10 #  Start of OS-specific
    STB_GNU_UNIQUE = 10 #  Unique symbol.
    STB_HIOS = 12 #  End of OS-specific
    STB_LOPROC = 13 #  Start of processor-specific
    STB_HIPROC = 15 #  End of processor-specific
end

# for st_info

@enum_export STType::UInt8 begin
    STT_NOTYPE = 0 #  Symbol type is unspecified
    STT_OBJECT = 1 #  Symbol is a data object
    STT_FUNC = 2 #  Symbol is a code object
    STT_SECTION = 3 #  Symbol associated with a section
    STT_FILE = 4 #  Symbol's name is file name
    STT_COMMON = 5 #  Symbol is a common data object
    STT_TLS = 6 #  Symbol is thread-local data object
    STT_NUM = 7 #  Number of defined types.
    # STT_LOOS = 10 #  Start of OS-specific
    STT_GNU_IFUNC = 10 #  Symbol is indirect code object
    STT_HIOS = 12 #  End of OS-specific
    STT_LOPROC = 13 #  Start of processor-specific
    STT_HIPROC = 15 #  End of processor-specific
end

# for st_other
@enum_export STVisibility::UInt8 begin
    STV_DEFAULT = 0 #  Default symbol visibility rules
    STV_INTERNAL = 1 #  Processor specific hidden class
    STV_HIDDEN = 2 #  Sym unavailable in other modules
    STV_PROTECTED = 3 #  Not preemptible, not exported
end

elf32_r_sym(val) = val >>> 8
elf32_r_type(val) = val & 0xff
elf32_r_info(sym, type) = (sym << 8) + (type & 0xff)

elf64_st_type(val) = STType(UInt8(val & 0xf))
elf64_r_sym(i) = i >>> 32
elf64_r_info(sym, type) = (Elf64_Xword(sym) << 32) + type

# for p_type
@enum_export PType::UInt32 begin
    PT_NULL = 0 #   Program header table entry unused
    PT_LOAD = 1 #   Loadable program segment
    PT_DYNAMIC = 2 #   Dynamic linking information
    PT_INTERP = 3 #   Program interpreter
    PT_NOTE = 4 #   Auxiliary information
    PT_SHLIB = 5 #   Reserved
    PT_PHDR = 6 #   Entry for header table itself
    PT_TLS = 7 #   Thread-local storage segment
    PT_NUM = 8 #   Number of defined types
    PT_LOOS = 0x60000000 #  Start of OS-specific
    PT_GNU_EH_FRAME = 0x6474e550 #  GCC .eh_frame_hdr segment
    PT_GNU_STACK = 0x6474e551 #  Indicates stack executability
    PT_GNU_RELRO = 0x6474e552 #  Read-only after relocation
    PT_LOSUNW = 0x6ffffffa #
    # PT_SUNWBSS = 0x6ffffffa #  Sun Specific segment
    PT_SUNWSTACK = 0x6ffffffb #  Stack segment
    PT_HISUNW = 0x6fffffff #
    # PT_HIOS = 0x6fffffff #  End of OS-specific
    PT_LOPROC = 0x70000000 #  Start of processor-specific
    PT_HIPROC = 0x7fffffff #  End of processor-specific
end

# for p_flags
@enum_export PFlags::UInt32 begin
    PF_X = 1
    PF_W = 2
    PF_R = 4
    PF_MASKOS = 0x0ff00000
    PF_MASKPROC = 0xf0000000
end

# values for note segmets
@enum_export NoteTypes::UInt32 begin
    NT_PRSTATUS = 1 # Contains copy of prstatus struct
    # NT_PRFPREG = 2 # Contains copy of fpregset struct.
    NT_FPREGSET = 2 # Contains copy of fpregset struct
    NT_PRPSINFO = 3 # Contains copy of prpsinfo struct
    NT_PRXREG = 4 # Contains copy of prxregset struct
    # NT_TASKSTRUCT = 4 # Contains copy of task structure
    NT_PLATFORM = 5 # String from sysinfo(SI_PLATFORM)
    NT_AUXV = 6 # Contains copy of auxv array
    NT_GWINDOWS = 7 # Contains copy of gwindows struct
    NT_ASRS = 8 # Contains copy of asrset struct
    NT_PSTATUS = 10 # Contains copy of pstatus struct
    NT_PSINFO = 13 # Contains copy of psinfo struct
    NT_PRCRED = 14 # Contains copy of prcred struct
    NT_UTSNAME = 15 # Contains copy of utsname struct
    NT_LWPSTATUS = 16 # Contains copy of lwpstatus struct
    NT_LWPSINFO = 17 # Contains copy of lwpinfo struct
    NT_PRFPXREG = 20 # Contains copy of fprxregset struct
    NT_SIGINFO = 0x53494749 # Contains copy of siginfo_t, size might increase
    NT_FILE = 0x46494c45 # Contains information about mapped files
    NT_PRXFPREG = 0x46e62b7f # Contains copy of user_fxsr_struct
    NT_PPC_VMX = 0x100 # PowerPC Altivec/VMX registers
    NT_PPC_SPE = 0x101 # PowerPC SPE/EVR registers
    NT_PPC_VSX = 0x102 # PowerPC VSX registers
    NT_PPC_TAR = 0x103 # Target Address Register
    NT_PPC_PPR = 0x104 # Program Priority Register
    NT_PPC_DSCR = 0x105 # Data Stream Control Register
    NT_PPC_EBB = 0x106 # Event Based Branch Registers
    NT_PPC_PMU = 0x107 # Performance Monitor Registers
    NT_PPC_TM_CGPR = 0x108 # TM checkpointed GPR Registers
    NT_PPC_TM_CFPR = 0x109 # TM checkpointed FPR Registers
    NT_PPC_TM_CVMX = 0x10a # TM checkpointed VMX Registers
    NT_PPC_TM_CVSX = 0x10b # TM checkpointed VSX Registers
    NT_PPC_TM_SPR = 0x10c # TM Special Purpose Registers
    NT_PPC_TM_CTAR = 0x10d # TM checkpointed Target Address Register
    NT_PPC_TM_CPPR = 0x10e # TM checkpointed Program Priority Register
    NT_PPC_TM_CDSCR = 0x10f # TM checkpointed Data Stream Control Register
    NT_PPC_PKEY = 0x110 # Memory Protection Keys registers.
    NT_386_TLS = 0x200 # i386 TLS slots (struct user_desc)
    NT_386_IOPERM = 0x201 # x86 io permission bitmap (1=deny)
    NT_X86_XSTATE = 0x202 # x86 extended state using xsave
    NT_S390_HIGH_GPRS = 0x300 # s390 upper register halves
    NT_S390_TIMER = 0x301 # s390 timer register
    NT_S390_TODCMP = 0x302 # s390 TOD clock comparator register
    NT_S390_TODPREG = 0x303 # s390 TOD programmable register
    NT_S390_CTRS = 0x304 # s390 control registers
    NT_S390_PREFIX = 0x305 # s390 prefix register
    NT_S390_LAST_BREAK = 0x306 # s390 breaking event address
    NT_S390_SYSTEM_CALL = 0x307 # s390 system call restart data
    NT_S390_TDB = 0x308 # s390 transaction diagnostic block
    NT_S390_VXRS_LOW = 0x309 # s390 vector registers 0-15 upper half.
    NT_S390_VXRS_HIGH = 0x30a # s390 vector registers 16-31.
    NT_S390_GS_CB = 0x30b # s390 guarded storage registers.
    NT_S390_GS_BC = 0x30c # s390 guarded storage broadcast control block.
    NT_S390_RI_CB = 0x30d # s390 runtime instrumentation.
    NT_ARM_VFP = 0x400 # ARM VFP/NEON registers
    NT_ARM_TLS = 0x401 # ARM TLS register
    NT_ARM_HW_BREAK = 0x402 # ARM hardware breakpoint registers
    NT_ARM_HW_WATCH = 0x403 # ARM hardware watchpoint registers
    NT_ARM_SYSTEM_CALL = 0x404 # ARM system call number
    NT_ARM_SVE = 0x405 # ARM Scalable Vector Extension registers
    NT_ARM_PAC_MASK = 0x406 # ARM pointer authentication code masks.
    NT_ARM_PACA_KEYS = 0x407 # ARM pointer authentication address keys.
    NT_ARM_PACG_KEYS = 0x408 # ARM pointer authentication generic key.
    NT_VMCOREDD = 0x700 # Vmcore Device Dump Note.
    NT_MIPS_DSP = 0x800 # MIPS DSP ASE registers.
    NT_MIPS_FP_MODE = 0x801 # MIPS floating-point mode.
    # NT_MIPS_MSA = 0x802 # MIPS SIMD registers.
end

const NT_VERSION = 1

@enum_export DTag::UInt32 begin
    DT_NULL = 0 # Marks end of dynamic section
    DT_NEEDED = 1 # Name of needed library
    DT_PLTRELSZ = 2 # Size in bytes of PLT relocs
    DT_PLTGOT = 3 # Processor defined value
    DT_HASH = 4 # Address of symbol hash table
    DT_STRTAB = 5 # Address of string table
    DT_SYMTAB = 6 # Address of symbol table
    DT_RELA = 7 # Address of Rela relocs
    DT_RELASZ = 8 # Total size of Rela relocs
    DT_RELAENT = 9 # Size of one Rela reloc
    DT_STRSZ = 10 # Size of string table
    DT_SYMENT = 11 # Size of one symbol table entry
    DT_INIT = 12 # Address of init function
    DT_FINI = 13 # Address of termination function
    DT_SONAME = 14 # Name of shared object
    DT_RPATH = 15 # Library search path (deprecated)
    DT_SYMBOLIC = 16 # Start symbol search here
    DT_REL = 17 # Address of Rel relocs
    DT_RELSZ = 18 # Total size of Rel relocs
    DT_RELENT = 19 # Size of one Rel reloc
    DT_PLTREL = 20 # Type of reloc in PLT
    DT_DEBUG = 21 # For debugging; unspecified
    DT_TEXTREL = 22 # Reloc might modify .text
    DT_JMPREL = 23 # Address of PLT relocs
    DT_BIND_NOW = 24 # Process relocations of object
    DT_INIT_ARRAY = 25 # Array with addresses of init fct
    DT_FINI_ARRAY = 26 # Array with addresses of fini fct
    DT_INIT_ARRAYSZ = 27 # Size in bytes of DT_INIT_ARRAY
    DT_FINI_ARRAYSZ = 28 # Size in bytes of DT_FINI_ARRAY
    DT_RUNPATH = 29 # Library search path
    DT_FLAGS = 30 # Flags for the object being loaded
    # DT_ENCODING = 32 # Start of encoded range
    DT_PREINIT_ARRAY = 32 # Array with addresses of preinit fct
    DT_PREINIT_ARRAYSZ = 33 # size in bytes of DT_PREINIT_ARRAY
    DT_SYMTAB_SHNDX = 34 # Address of SYMTAB_SHNDX section
    DT_NUM = 35 # Number used
    DT_LOOS = 0x6000000d # Start of OS-specific
    DT_HIOS = 0x6ffff000 # End of OS-specific
    DT_LOPROC = 0x70000000 # Start of processor-specific
    DT_HIPROC = 0x7fffffff # End of processor-specific
    # DT_PROCNUM = DT_MIPS_NUM # Most used by any processor
end

@enum_export DPtr::UInt32 begin
    DT_ADDRRNGLO = 0x6ffffe00
    DT_GNU_HASH = 0x6ffffef5 # GNU-style hash table.
    DT_TLSDESC_PLT = 0x6ffffef6
    DT_TLSDESC_GOT = 0x6ffffef7
    DT_GNU_CONFLICT = 0x6ffffef8 # Start of conflict section
    DT_GNU_LIBLIST = 0x6ffffef9 # Library list
    DT_CONFIG = 0x6ffffefa # Configuration information.
    DT_DEPAUDIT = 0x6ffffefb # Dependency auditing.
    DT_AUDIT = 0x6ffffefc # Object auditing.
    DT_PLTPAD = 0x6ffffefd # PLT padding.
    DT_MOVETAB = 0x6ffffefe # Move table.
    DT_SYMINFO = 0x6ffffeff # Syminfo table.
    # DT_ADDRRNGHI = 0x6ffffeff
    DT_ADDRNUM = 11
end

dt_addrtagidx(tag) = (DT_ADDRRNGHI - (tag)) # Reverse order!

@enum_export DVal::UInt32 begin
    DT_VALRNGLO = 0x6ffffd00
    DT_GNU_PRELINKED = 0x6ffffdf5 # Prelinking timestamp
    DT_GNU_CONFLICTSZ = 0x6ffffdf6 # Size of conflict section
    DT_GNU_LIBLISTSZ = 0x6ffffdf7 # Size of library list
    DT_CHECKSUM = 0x6ffffdf8
    DT_PLTPADSZ = 0x6ffffdf9
    DT_MOVEENT = 0x6ffffdfa
    DT_MOVESZ = 0x6ffffdfb
    DT_FEATURE_1 = 0x6ffffdfc # Feature selection (DTF_*).
    DT_POSFLAG_1 = 0x6ffffdfd # Flags for DT_* entries, effecting the following DT_* entry.
    DT_SYMINSZ = 0x6ffffdfe # Size of syminfo table (in bytes)
    DT_SYMINENT = 0x6ffffdff # Entry size of syminfo
    # DT_VALRNGHI = 0x6ffffdff
    DT_VALNUM = 12

end

dt_valtagidx(tag) = (DT_VALRNGHI - (tag)) # Reverse order!


text:
module LaserFileSplitter

using EzXML

function save_svg_in_chunks(filename::AbstractString, outfolder::AbstractString, max_num_chunks::Integer)
    svgdoc = readxml(filename);
    elms = collect(eachelement(root(svgdoc)))
    paths = filter(n->n.name == "path", elms)

    mkdir(outfolder)

    for num_chunks in 1:max_num_chunks

        chunk_folder = outfolder*"/$(num_chunks)_chunks"
        mkdir(chunk_folder)

        num_per_chunk = Int(ceil(length(paths) // num_chunks))
        chunked_paths = [paths[1 + (k-1)*num_per_chunk:k*num_per_chunk] for k in 1:num_chunks-1]
        push!(chunked_paths, paths[1 + (num_chunks-1)*num_per_chunk:end])

        @assert(sum(length.(chunked_paths))==length(paths))

        for k in 1:num_chunks
            # delete all other paths
            for i in 1:num_chunks
                if i != k
                    unlink!.(chunked_paths[i])
                end
            end

            # save chunk
            write(chunk_folder*"/$(num_chunks)_$(k).svg", svgdoc)

            # add back all other paths
            for i in 1:num_chunks
                if i != k
                    for path in chunked_paths[i]
                        link!(root(svgdoc), path)
                    end
                end
            end

        end
    end
    true
end

export save_svg_in_chunks

end

max_stars_repo_name: zoemcc/LaserFileSplitter

max_issues_repo_name: zoemcc/LaserFileSplitter

max_forks_repo_name: zoemcc/LaserFileSplitter


text:
"""
# AK135

## Reference

Kennett, B.L.N., Engdahl, E., Buland, R., 1995. Constraints on seismic velocities in
the Earth from travel-times. Geophys J Int 122, 108–124.
https://doi.org/10.1111/j.1365-246X.1995.tb03540.x
"""
const AK135 = LinearLayeredModel(6371., 136,
    # radius
    6371. .- reverse([ 0.0, 20.0, 20.0, 35.0, 35.0, 77.5, 120.0, 165.0, 210.0, 210.0,
    260.0, 310.0, 360.0, 410.0, 410.0, 460.0, 510.0, 560.0, 610.0, 660.0, 660.0, 710.0,
    760.0, 809.5, 859.0, 908.5, 958.0, 1007.5, 1057.0, 1106.5, 1156.0, 1205.5, 1255.0, 1304.5, 1354.0,
    1403.5, 1453.0, 1502.5, 1552.0, 1601.5, 1651.0, 1700.5, 1750.0, 1799.5, 1849.0,
    1898.5, 1948.0, 1997.5, 2047.0, 2096.5, 2146.0, 2195.5, 2245.0, 2294.5, 2344.0,
    2393.5, 2443.0, 2492.5, 2542.0, 2591.5, 2640.0, 2690.0, 2740.0, 2740.0, 2789.7,
    2839.3, 2891.5, 2891.5, 2939.3, 2989.7, 3040.0, 3090.3, 3140.7, 3191.0, 3241.3,
    3291.7, 3342.0, 3392.3, 3442.6, 3493.0, 3543.3, 3593.6, 3644.0, 3694.3, 3744.6,
    3795.0, 3845.3, 3895.6, 3945.9, 3996.3, 4046.6, 4096.9, 4147.3, 4197.6, 4247.9,
    4298.3, 4348.6, 4398.9, 4449.3, 4499.6, 4549.9, 4600.3, 4650.6, 4700.9, 4801.6,
    4851.9, 4902.2, 4952.6, 5002.9, 5053.2, 5103.6, 5153.5, 5153.5, 5204.6, 5255.3,
    5306.0, 5356.8, 5407.5, 5458.2, 5508.9, 5559.6, 5610.3, 5661.0, 5711.7, 5813.2,
    5863.9, 5914.6, 5965.3, 6016.0, 6066.7, 6117.4, 6168.1, 6218.9, 6269.6, 6320.3,
    6371.0]),
    # Vp
    reverse([ 5.8000, 5.8000, 6.5000, 6.5000, 8.0400, 8.0450, 8.0500, 8.1750, 8.3000,
    8.3000, 8.4825, 8.6650, 8.8475, 9.0300, 9.3600, 9.5280, 9.6960, 9.8640, 10.0320,
    10.2000, 10.7900, 10.9229, 11.0558, 11.1353, 11.2221, 11.3068, 11.3896, 11.4705,
    11.5495, 11.6269, 11.7026, 11.7766, 11.8491, 11.9200, 11.9895, 12.0577, 12.1245,
    12.1912, 12.2550, 12.3185, 12.3819, 12.4426, 12.5031, 12.5631, 12.6221, 12.6804,
    12.7382, 12.7956, 12.8526, 12.9096, 12.9668, 13.0222, 13.0783, 13.1336, 13.1894,
    13.2465, 13.3018, 13.3585, 13.4156, 13.4741, 13.5312, 13.5900, 13.6494, 13.6494,
    13.6530, 13.6566, 13.6602, 8.0000, 8.0382, 8.1283, 8.2213, 8.3122, 8.4001, 8.4861,
    8.5692, 8.6496, 8.7283, 8.8036, 8.8761, 8.9461, 9.0138, 9.0792, 9.1426, 9.2042,
    9.2634, 9.3205, 9.3760, 9.4297, 9.4814, 9.5306, 9.5777, 9.6232, 9.6673, 9.7100,
    9.7513, 9.7914, 9.8304, 9.8682, 9.9051, 9.9410, 9.9761, 10.0103, 10.0439, 10.0768,
    10.1415, 10.1739, 10.2049, 10.2329, 10.2565, 10.2745, 10.2854, 10.2890, 11.0427,
    11.0585, 11.0718, 11.0850, 11.0983, 11.1166, 11.1316, 11.1457, 11.1590, 11.1715,
    11.1832, 11.1941, 11.2134, 11.2219, 11.2295, 11.2364, 11.2424, 11.2477, 11.2521,
    11.2557, 11.2586, 11.2606, 11.2618, 11.2622]),
    # Vs
    reverse([3.4600, 3.4600, 3.8500, 3.8500, 4.4800, 4.4900, 4.5000, 4.5090, 4.5180,
    4.5230, 4.6090, 4.6960, 4.7830, 4.8700, 5.0800, 5.1860, 5.2920, 5.3980, 5.5040,
    5.6100, 5.9600, 6.0897, 6.2095, 6.2426, 6.2798, 6.3160, 6.3512, 6.3854, 6.4187,
    6.4510, 6.4828, 6.5138, 6.5439, 6.5727, 6.6008, 6.6285, 6.6555, 6.6815, 6.7073,
    6.7326, 6.7573, 6.7815, 6.8052, 6.8286, 6.8515, 6.8742, 6.8972, 6.9194, 6.9418,
    6.9627, 6.9855, 7.0063, 7.0281, 7.0500, 7.0720, 7.0931, 7.1144, 7.1369, 7.1586,
    7.1807, 7.2031, 7.2258, 7.2490, 7.2490, 7.2597, 7.2704, 7.2811, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 3.5043, 3.5187, 3.5314, 3.5435, 3.5551, 3.5661, 3.5765,
    3.5864, 3.5957, 3.6044, 3.6126, 3.6202, 3.6337, 3.6396, 3.6450, 3.6498, 3.6540,
    3.6577, 3.6608, 3.6633, 3.6653, 3.6667, 3.6675, 3.6678]),
    # ρ
    reverse([2.44900, 2.44900, 2.71420, 2.71420, 3.29760, 3.29940, 3.30130, 3.34870,
    3.39600, 3.39600, 3.46520, 3.53430, 3.60340, 3.67260, 3.79760, 3.86120, 3.92480,
    3.98850, 4.05210, 4.11580, 4.33930, 4.38960, 4.43990, 4.47010, 4.50290, 4.53500,
    4.56640, 4.59700, 4.62700, 4.65630, 4.68490, 4.71300, 4.74040, 4.76730, 4.79360,
    4.81950, 4.84480, 4.87000, 4.89420, 4.91820, 4.94230, 4.96530, 4.98820, 5.01090,
    5.03330, 5.05530, 5.07720, 5.09900, 5.12060, 5.14220, 5.16380, 5.18480, 5.20610,
    5.22700, 5.24810, 5.26980, 5.29070, 5.31220, 5.33380, 5.35600, 5.37760, 5.39990,
    5.42240, 5.42240, 5.42380, 5.42510, 5.42650, 9.91450, 9.99420, 10.07220, 10.14850,
    10.22330, 10.29640, 10.36790, 10.43780, 10.50620, 10.57310, 10.63850, 10.70230,
    10.76470, 10.82570, 10.88520, 10.94340, 11.00010, 11.05550, 11.10950, 11.16230,
    11.21370, 11.26390, 11.31270, 11.36040, 11.40690, 11.45210, 11.49620, 11.53910,
    11.58090, 11.62160, 11.66120, 11.69980, 11.73730, 11.77370, 11.80920, 11.84370,
    11.87720, 11.94140, 11.97220, 12.00010, 12.03110, 12.05930, 12.08670, 12.11330,
    12.13910, 12.70370, 12.72890, 12.75300, 12.77600, 12.79800, 12.81880, 12.83870,
    12.85740, 12.87510, 12.89170, 12.90720, 12.92170, 12.94740, 12.95860, 12.96880,
    12.97790, 12.98590, 12.99290, 12.99880, 13.00360, 13.00740, 13.01000, 13.01170,
    13.01220]),
    false, Float64[], Float64[], Float64[], Float64[], Float64[], false, Float64[], Float64[])


text:
function sparse_cube(datacube)
    radiance = []
    for i in 1:size(datacube)[3]
        push!(radiance,sparse(datacube[:,:,i]))
    end
    return VectorOfArray(radiance)
end


text:
@compile reference_copy """
__global__ void reference_copy(const float *input, float *output)
{
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    output[i] = input[i];
}
"""

@target ptx function kernel_copy(input::CuDeviceArray{Float32}, output::CuDeviceArray{Float32})
    i = blockId_x() + (threadId_x()-1) * numBlocks_x()
    output[i] = input[i]

    return nothing
end


text:
abstract type StepSize end

abstract type State{P} end

struct PimcSystem
    pot::Function
    force::Function
    mass::Float64
    beta::Float64

    function PimcSystem(pot::Function, force::Function, mass::Float64, temperature::Float64)
        beta = 1 / (kB * temperature)
        new(pot, force, mass, beta)
    end
end

bead_prev(P::Int, j::Int) = mod(j - 1, 1:P)
bead_next(P::Int, j::Int) = mod(j + 1, 1:P)


text:
using DataGraphs
using NamedGraphs
using Dictionaries
using Graphs

g = NamedDimGraph(grid((2, 2)); dims=(2, 2))
dg = NamedDimDataGraph{String,String}(g)

for v in vertices(dg)
  dg[v] = "V$v"
end
for e in edges(dg)
  dg[e] = "E$e"
end

dg_1c = dg[1, :]
dg_2c = dg[2, :]
dg_c1 = dg[:, 1]
dg_c2 = dg[:, 2]

@show nv(dg) == 4
@show ne(dg) == 4
@show nv(dg_1c) == 2
@show ne(dg_1c) == 1
@show nv(dg_2c) == 2
@show ne(dg_2c) == 1
@show nv(dg_c1) == 2
@show ne(dg_c1) == 1
@show nv(dg_c2) == 2
@show ne(dg_c2) == 1

dg_nn = dg[[(1, 1), (2, 2)]]

@show nv(dg_nn) == 2
@show ne(dg_nn) == 0


text:
using GeometricIntegrators.Solvers
using Test

n = 1
T = Float64

function F(x::Vector, b::Vector)
    b[:] = x.^2
end

function J(x::Vector, A::Matrix)
    A[:,:] = 2x
end

for Solver in (NewtonSolver, QuasiNewtonSolver, NLsolveNewton)
    x = ones(T, n)
    nl = Solver(x, F)
    solve!(nl)
    # println(nl.status.i, ", ", nl.status.rₐ,", ",  nl.status.rᵣ,", ",  nl.status.rₛ)
    for x in nl.x
        @test x ≈ 0 atol=1E-7
    end

    x = ones(T, n)
    nl = Solver(x, F; J! = J)
    solve!(nl)
    # println(nl.status.i, ", ", nl.status.rₐ,", ",  nl.status.rᵣ,", ",  nl.status.rₛ)
    for x in nl.x
        @test x ≈ 0 atol=1E-7
    end
end


text:
"""
`addfilter(t; readout = true)`

Create selectors (`categoricalselector`, `rangeselector`, `selector` are supported) and delete them for various
columns of table `t`. `readout` denotes whether the table will be displayed initially. Outputs the filtered table.
"""
function addfilter(t, r = 6; readout = true)
    t isa AbstractObservable || (t = Observable{Any}(t))

    wdg = Widget{:addfilter}(output = Observable{Any}(t[]))
    cols = map(Tables.columntable, t)

    wdg[:cols] = dropdown(map(collect∘propertynames, cols), placeholder = "Column to filter", value = nothing)

    selectoptions = OrderedDict(
        "categorical" => categoricalselector,
        "range" => rangeselector,
        "predicate" => selector
    )

    wdg[:selectortype] = dropdown(selectoptions, placeholder = "Selector type", value = nothing)
    wdg[:button] = button("Add selector")
    wdg[:filter] = button("Filter")

    wrap = Widgets.div(className = "column")
    container = Widgets.div(className = "columns is-multiline is-mobile")
    wdg[:selectors] = notifications([], wrap = wrap, container = container)

    @on begin
        &wdg[:button]
        push!(wdg[:selectors][], wdg[:selectortype][](cols[], wdg[:cols][], lazymap))
        wdg[:selectors][] = wdg[:selectors][]
    end

    on(wdg[:filter]) do x
        sels = (observe(i)[] for i in observe(wdg[:selectors])[])
        wdg.output[] = Tables.materializer(t[])(_filter(cols[], sels...))
    end

    @layout! wdg Widgets.div(
        Widgets.div(className = "level is-mobile", map(Widgets.div(className="level-item"), [:cols, :selectortype, :button, :filter])...),
        :selectors,
        toggled(head(_.output, r), readout = true),
        className = "interact-widget"
    )
    on(t) do x
        observe(wdg, :selectors)[] = []
        observe(wdg)[] = x
    end
    wdg
end

addfilter(t::Widgets.AbstractWidget; kwargs...) = addfilter(observe(t); kwargs...)


text:
using PatModules
@mainmodule example ("subpackage", "utils") begin

    function greetings()
        utils.do_nothing()
        subpackage.say_hi()
        # Can descend into subsubpackages
        subpackage.two.say_bye()
    end
end


text:
using StatsBase, Distributions, Distances
export SPF, DcSPF, SpfPosterior
export cardinality, get_normalising_const

struct SPF{T<:Union{Int, String}}
    mode::Path{T} # Mode
    γ::Real # Precision
    dist::Metric # Distance metric
    V::Vector{T} # Vertex Set
    K::Real # Maximum length path
end

SPF(Iᵐ, γ, d, 𝒱) = SPF(Iᵐ, γ, d, 𝒱, Inf)  # K defaults to ∞
SPF(Iᵐ, γ, d) = SPF(Iᵐ, γ, d, unique(Iᵐ), Inf)  # Vertex set dedaults to unique vals in mods

struct DcSPF
    mode::Path{Int} # Mode
    γ::Real         # Dispersion
    d::Metric
    p::DiscreteUnivariateDistribution  # Path length dist
    V::UnitRange
end
# K defaults to ∞
DcSPF(Iᵐ, γ, d, p) = SPF(Iᵐ, γ, d, p, unique(Iᵐ))  # Vertex set dedaults to unique vals in mods

struct SpfPosterior{T<:Union{Int, String}}
    data::Vector{Path{T}}
    I_prior::SPF{T}
    γ_prior::ContinuousUnivariateDistribution
    dist::Metric
    V::Vector{T}
    K::Real
    sample_size::Int
    function SpfPosterior(
        data::Vector{Path{S}},
        I_prior::SPF{S},
        γ_prior::ContinuousUnivariateDistribution
        ) where {S<:Union{Int, String}}

        dist = I_prior.dist
        V = I_prior.V
        K = I_prior.K
        sample_size = length(data)
        new{S}(data, I_prior, γ_prior, dist, V, K, sample_size)
    end
end

function log_eval(p::SpfPosterior, Z::Float64, Iᵐ::Path, γ::Float64)
    log_posterior = (
        -γ * sum_of_dists(p.data, Iᵐ, p.dist) # Unormalised log likleihood
        - p.sample_size * log(Z) # Normalising constant
        - p.I_prior.γ * p.dist(Iᵐ, p.I_prior.mode) # Iᵐ prior (unormalised)
        + logpdf(p.γ_prior, γ)

    )
    return log_posterior

end

function get_length_dist(p::SpfPosterior{T}, α::Float64) where T <:Union{Int, String}

    lprobs = counts(length.(p.data), 1:p.K) .+ α * length(p.data)
    lprobs = lprobs / sum(lprobs)

    return Categorical(lprobs)
end

function get_vertex_proposal_dist(p::SpfPosterior{T}) where {T<:Union{Int, String}}
    μ = vertex_counts(vcat(p.data...), p.V)
    μ /= sum(μ)
    return μ
end

max_stars_repo_path: src/Metric Models/SPF.jl

max_issues_repo_path: src/Metric Models/SPF.jl

max_forks_repo_path: src/Metric Models/SPF.jl


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

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

# ╔═╡ 972b2230-7634-11eb-028d-df7fc722ec70
html"""

Section 1.5

 Transformations II: Composability, Linearity and Nonlinearity

Lecture Video

"""

# ╔═╡ b7895bd2-7634-11eb-211e-ef876d23bd88
PlutoUI.TableOfContents(aside=true)

# ╔═╡ 230b0118-30b7-4035-ad31-520165a76fcc
md"""
#### Intializing packages

_When running this notebook for the first time, this could take up to 15 minutes. Hang in there!_
"""

# ╔═╡ 6b473b2d-4326-46b4-af38-07b61de287fc
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(name="ImageIO", version="0.5"),
		Pkg.PackageSpec(name="ImageShow", version="0.2"),
		Pkg.PackageSpec(name="FileIO", version="1.6"),
		Pkg.PackageSpec(name="PNGFiles", version="0.3.6"),
		Pkg.PackageSpec(name="Colors", version="0.12"),
		Pkg.PackageSpec(name="ColorVectorSpace", version="0.8"),

		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
		Pkg.PackageSpec(name="HypertextLiteral", version="0.5"),
		Pkg.PackageSpec(name="ForwardDiff", version="0.10"),
		Pkg.PackageSpec(name="NonlinearSolve", version="0.3"),
		Pkg.PackageSpec(name="StaticArrays"),
	])

	using PlutoUI
	using Colors, ColorVectorSpace, ImageShow, FileIO
	using PlutoUI
	using HypertextLiteral
	using LinearAlgebra
	using ForwardDiff
	using NonlinearSolve
	using StaticArrays
end

# ╔═╡ 230cba36-9d0a-4726-9e55-7df2c6743968
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ 96766502-7a06-11eb-00cc-29849773dbcf
# img_original = load(download(corgis_url));
img_original = load(download(longcorgi_url));
# img_original = load(download(theteam_url));

# ╔═╡ 890d30b9-2cd0-4d3a-99f6-f7d3d7858fda
corgis_url = "https://user-images.githubusercontent.com/6933510/108605549-fb28e180-73b4-11eb-8520-7e29db0cc965.png"

# ╔═╡ 85fba8fb-a9ea-444d-831b-ec6489b58b4f
longcorgi_url = "https://user-images.githubusercontent.com/6933510/110868198-713faa80-82c8-11eb-8264-d69df4509f49.png"

# ╔═╡ 06beabc3-2aa7-4e78-9bae-dc4b37251aa2
theteam_url = "https://news.mit.edu/sites/default/files/styles/news_article__image_gallery/public/images/202004/edelman%2520philip%2520sanders.png?itok=ZcYu9NFeg"

# ╔═╡ 26dd0e98-7a75-11eb-2196-5d7bda201b19
md"""
After you select your image, we suggest moving this line above just above the top of your browser.

---------------
"""

# ╔═╡ e0b657ce-7a03-11eb-1f9d-f32168cb5394
md"""
#  The fun stuff: playing with transforms
"""

# ╔═╡ 005ca75a-7622-11eb-2ba4-9f450e71df1f
let

range = -1.5:.1:1.5
md"""

This is a "scrubbable" matrix: click on the number and drag to change!

**A =**

``(``
 $(@bind a Scrubbable( range; default=1.0))
 $(@bind b Scrubbable( range; default=0.0))
``)``

``(``
$(@bind c Scrubbable(range; default=0.0 ))
$(@bind d Scrubbable(range; default=1.0))
``)``

"""
end

# ╔═╡ 23ade8ee-7a09-11eb-0e40-296c6b831d74
md"""
Grab a [linear](#a0afe3ae-76b9-11eb-2301-cde7260ddd7f) or [nonlinear](#a290d5e2-7a02-11eb-37db-41bf86b1f3b3) transform, or make up your own!
"""

# ╔═╡ 58a30e54-7a08-11eb-1c57-dfef0000255f
# T⁻¹ = id
#  T⁻¹ = rotate(α)
  T⁻¹ = shear(α)
#   T⁻¹ = lin(A) # uses the scrubbable
#   T⁻¹ = shear(α) ∘ shear(-α)
 # T⁻¹ = nonlin_shear(α)
 #   T⁻¹ =   inverse(nonlin_shear(α))
#    T⁻¹ =  nonlin_shear(-α)
#  T⁻¹ =  xy
# T⁻¹ = warp(α)
# T⁻¹ = ((x,y),)-> (x+α*y^2,y+α*x^2) # may be non-invertible

#T⁻¹ = ((x,y),)-> (x,y^2)
# T⁻¹  = flipy ∘ ((x,y),) ->  ( (β*x - α*y)/(β - y)  , -h*y/ (β - y)   )

# ╔═╡ 2efaa336-7630-11eb-0c17-a7d4a0141dac
md"""
zoom = $(@bind  z Scrubbable(.1:.1:3,  default=1))
"""

# ╔═╡ 7f28ac40-7914-11eb-1403-b7bec34aeb94
md"""
pan = [$(@bind panx Scrubbable(-1:.1:1, default=0)),
$(@bind pany Scrubbable(-1:.1:1, default=0)) ]
"""

# ╔═╡ ce55beee-7643-11eb-04bc-b517703facff
md"""
α= $(@bind α Slider(-30:.1:30, show_value=true, default=0))
β= $(@bind β Slider(-10:.1:10, show_value=true, default = 5))
h= $(@bind h Slider(.1:.1:10, show_value=true, default = 5))
"""

# ╔═╡ b76a5bd6-802f-11eb-0951-1f1092dee8de
1+1

# ╔═╡ 5d33f6ea-7e9c-11eb-2fb3-dbb7cb07c60c
md"""
pixels = $(@bind pixels Slider(1:1000, default=800, show_value=true))
"""

# ╔═╡ 45dccdec-7912-11eb-01b4-a97e30344f39
md"""
Show grid lines $(@bind show_grid CheckBox(default=true))
ngrid = $(@bind ngrid Slider(5:5:20, show_value=true, default = 10))
"""

# ╔═╡ d2fb356e-7f32-11eb-177d-4f47d6c9e59b
md"""
Circular Frame $(@bind circular CheckBox(default=true))
radius = $(@bind r Slider(.1:.1:1, show_value=true, default = 1))
"""

# ╔═╡ ca28189e-7e9a-11eb-21d6-bd819f3e0d3a
begin
		[
			begin

			 x,y = transform_ij_to_xy(i,j, pixels)

			X,Y = ( translate(-panx,-pany)  )([x,y])
			 X,Y = ( T⁻¹∘scale(1/z)∘translate(-panx,-pany) )([x,y])
			 i,j = transform_xy_to_ij(img,X,Y)
			 getpixel(img,i,j; circular=circular, r=r)
			end

			for i = 1:pixels, j = 1:pixels
		]
end

# ╔═╡ ccea7244-7f2f-11eb-1b7b-b9b8473a8c74
transform_xy_to_ij(img,0.0,0.0)

# ╔═╡ 55b5fc92-7a76-11eb-3fba-854c65eb87f9
md"""
Above: The original image is placed in a [-1,1] x [-1 1] box and transformed.
"""

# ╔═╡ 85686412-7a75-11eb-3d83-9f2f8a3c5509
A = [a b ; c d];

# ╔═╡ a7df7346-79f8-11eb-1de6-71f027c46643
md"""
## Pedagogical note: Why the Module 1 application = image processing

Image processing is a great way to learn Julia, some linear algebra, and some nonlinear mathematics.  We don't presume the audience will become professional image processors, but we do believe that the principles learned transcend so many applications... and everybody loves playing with their own images!
"""

# ╔═╡ 044e6128-79fe-11eb-18c1-395ae857dc73
md"""
# Last Lecture Leftovers
"""

# ╔═╡ 78d61e28-79f9-11eb-0605-e77d206cda84
md"""
## Interesting question about linear transformations
If a transformation takes lines into lines and preserves the origin, does it have to be linear?

Answer = **no**!

The example of a **perspective map** takes all lines into lines, but paralleograms generally do not become parallelograms.
"""

# ╔═╡ aad4d6e4-79f9-11eb-0342-b900a41cfbaf
md"""
[A nice interactive demo of perspective maps](https://www.khanacademy.org/humanities/renaissance-reformation/early-renaissance1/beginners-renaissance-florence/a/linear-perspective-interactive) from Khan academy.
"""

# ╔═╡ d42aec08-76ad-11eb-361a-a1f2c90fd4ec
Resource("https://cdn.kastatic.org/ka-perseus-images/1b351a3653c1a12f713ec24f443a95516f916136.jpg")

# ╔═╡ d9115c1a-7aa0-11eb-38e4-d977c5a6b75b
md"""
**Challenge exercise**: Rewrite this using Julia and Pluto!
"""

# ╔═╡ e965cf5e-79fd-11eb-201d-695b54d08e54
md"""
## Julia style (a little advanced): Reminder about defining vector valued functions
"""

# ╔═╡ 1e11c1ec-79fe-11eb-1867-9da72b3f3bc4
md"""
Many people find it hard to read

`f(v) = [ v[1]+v[2] , v[1]-v[2] ]  ` or
`  f = v ->  [ v[1]+v[2] , v[1]-v[2] ]  `

and instead prefer

`f((x,y)) = [ x+y , x-y ] ` or
` f = ((x,y),) -> [ x+y , x-y ] `.

All four of these will take a 2-vector to a 2-vector in the same way for the purposes of this lecture, i.e. `f( [1,2] )` can be defined by any of the four forms.

The forms with the `->` are anonymous functions.  (They are still
considered anonymous, even though we then name them `f`.)
"""

# ╔═╡ 28ef451c-7aa1-11eb-340c-ab3a1193a3c4
md"""
## Functions with parameters
The anonymous form comes in handy when one wants a function to depend on a **parameter**.
For example:

`f(α) = ((x,y),) -> [x + αy, x - αy]`

allows you to apply the `f(7)` function to the input vector `[1, 2]` by running
`f(7)([1, 2])` .
"""

# ╔═╡ a0afe3ae-76b9-11eb-2301-cde7260ddd7f
md"""
# Linear transformations: a collection
"""

# ╔═╡ fc2deb7c-7aa1-11eb-019f-d3e3c80b9ff1
md"""
Here are a few useful linear transformations:
"""

# ╔═╡ d364f91a-76b9-11eb-1807-75e733940d53
begin
	 id((x,y)) = SA[x, y]

	 scalex(α) = ((x, y),) -> SA[α*x,  y]
	 scaley(α) = ((x, y),) -> SA[x,   α*y]
	 scale(α)  = ((x, y),) -> SA[α*x, α*y]

	 swap((x,y))  = SA[y, x]
	 flipy((x,y)) = SA[x, -y]

	 rotate(θ) = ((x, y),) -> SA[cos(θ)*x + sin(θ)*y, -sin(θ)*x + cos(θ)*y]
	 shear(α)  = ((x, y),) -> SA[x + α*y, y]
end

# ╔═╡ 080d87e0-7aa2-11eb-18f5-2fb6a7a5bcb4
md"""
In fact we can write down the *most general* linear transformation in one of two ways:
"""

# ╔═╡ 15283aba-7aa2-11eb-389c-e9f215bd03e2
begin
	lin(a, b, c, d) = ((x, y),) -> ( a*x + b*y, c*x + d*y )

	lin(A) = v-> A * [v...]  # linear algebra version using matrix multiplication
end

# ╔═╡ 2612d2c2-7aa2-11eb-085a-1f27b6174995
md"""
The second version uses the matrix multiplication notation from linear algebra, which means exactly the same as the first version when

$$A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$$
"""

# ╔═╡ a290d5e2-7a02-11eb-37db-41bf86b1f3b3
md"""
# Nonlinear transformations: a collection
"""

# ╔═╡ b4cdd412-7a02-11eb-149a-df1888a0f465
begin
  translate(α,β)  = ((x, y),) -> SA[x+α, y+β]   # affine, but not linear

  nonlin_shear(α) = ((x, y),) -> SA[x, y + α*x^2]

  warp(α)    = ((x, y),) -> rotate(α*√(x^2+y^2))(SA[x, y])
  xy((r, θ)) = SA[ r*cos(θ), r*sin(θ) ]
  rθ(x)      = SA[norm(x), atan(x[2],x[1]) ]

  # exponentialish =  ((x,y),) -> [log(x+1.2), log(y+1.2)]
  # merc = ((x,y),) ->  [ log(x^2+y^2)/2 , atan(y,x) ] # (reim(log(complex(y,x)) ))
end

# ╔═╡ 704a87ec-7a1e-11eb-3964-e102357a4d1f
md"""
# Composition
"""

# ╔═╡ 4b0e8742-7a70-11eb-1e78-813f6ad005f4
let
	x = rand()

	( sin ∘ cos )(x) ≈ sin(cos(x))
end

# ╔═╡ 44792484-7a20-11eb-1c09-95b27b08bd34
md"""
## Composing functions in mathematics
[Wikipedia (math) ](https://en.wikipedia.org/wiki/Function_composition)

In math we talk about *composing* two functions to create a new function:
the function that takes $x$ to $\sin(\cos(x))$ is the **composition**
of the sine function and the cosine function.

We humans tend to blur the distinction between the sine function
and the value of $\sin(x)$ at some point $x$.  The sine function
is a mathematical object by itself.  It's a thing that can be evaluated
at as many $x$'s as you like.

If you look at the two sides of
` (sin ∘ cos )(x) ≈ sin(cos(x)) ` and see that they are exactly the same, it's time to ask yourself what's a key difference? On the left a function is built
` sin ∘ cos ` which is then evaluated at `x`.  On the right, that function
is never built.
"""

# ╔═╡ f650b788-7a70-11eb-0b20-779d2f18f111
md"""
## Composing functions in computer science
[wikipedia (cs)](https://en.wikipedia.org/wiki/Function_composition_(computer_science))

A key issue is a programming language is whether it's easy to name
the composition in that language.  In Julia one can create the function
`sin ∘ cos`  and one can readily check that ` (sin ∘ cos)(x) ` always yields the same value as `sin(cos(x))`.
"""

# ╔═╡ c852d398-7aa2-11eb-2ded-ab2e5236e9b2
md"""

## Composing functions in Julia
[Julia's  `∘` operator](https://docs.julialang.org/en/v1/manual/functions/#Function-composition-and-piping)
follows the [mathematical typography](https://en.wikipedia.org/wiki/Function_composition#Typography) convention, as was
shown in the `sin ∘ cos` example above. We can type this symbol as `\circ<TAB>`.

"""

# ╔═╡ 061076c2-7aa3-11eb-0d04-b7cbc60e6cb2
md"""
## Composition of software at a higher level

The trend these days is to have higher-order composition of functionalities.
A good example would be that an optimization can wrap a highly complicated
program which might include all kinds of solvers, and still run successfully.
This might require the ability of the outer software to have some awareness
of the inner software.  It can be quite magical when two very different pieces of software "compose", i.e. work together.  Julia's language construction encourages composability.  We will discuss this more in a future lecture.

"""

# ╔═╡ 014c14a6-7a72-11eb-119b-f5cfc82085ca
md"""
### Find your own examples

Take some of the Linear and Nonlinear Transformations (see the Table of Contents) and find some inverses by placing them in the `T=` section of "The fun stuff" at the top of this notebook.
"""

# ╔═╡ 89f0bc54-76bb-11eb-271b-3190b4d8cbc0
md"""
Linear transformations can be written in math using matrix multiplication notation as

$$\begin{pmatrix} a & b \\ c & d \end{pmatrix}
\begin{pmatrix} x \\ y \end{pmatrix}$$.
"""

# ╔═╡ f70f7ea8-76b9-11eb-3bd7-87d40a2861b1
md"""
By contrast, here are a few fun functions that cannot be written as matrix times
vector.  What characterizes the matrix ones from the non-matrix ones?
"""

# ╔═╡ bf28c388-76bd-11eb-08a7-af2671218017
md"""
This may be a little fancy, but we see that warp is a rotation, but the
rotation depends on the point where it is applied.
"""

# ╔═╡ 5655d2a6-76bd-11eb-3042-5b2dd3f6f44e
begin
	warp₂(α,x,y) = rotate(α*√(x^2+y^2))
	warp₂(α) = ((x,y),) -> warp₂(α,x,y)([x,y])
end

# ╔═╡ 56f1e4cc-7a03-11eb-187b-c5a917978eb9
warp3(α) = ((x,y),) -> rotate(α*√(x^2+y^2))([x,y])

# ╔═╡ 70dc4346-7a03-11eb-055e-111d2519a44c
warp3(1)([1,2])

# ╔═╡ 852592d6-76bd-11eb-1265-5f200e39113d
warp(1)([5,6])

# ╔═╡ 8e36f4a2-76bd-11eb-2fda-9d1424752812
warp₂(1.0)([5.0,6.0])

# ╔═╡ a8bf7128-7aa5-11eb-3ee9-953b0b5ccd01

# ╔═╡ ad700740-7a74-11eb-3369-15e5fd89194d
md"""
# Linear transformations: See a matrix, think beyond number arrays
"""

# ╔═╡ e051259a-7a74-11eb-12fc-99c5dc867fbd
md"""
Software writers and beginning linear algebra students see a matrix and see a lowly table of numbers.  We want you to see a *linear transformation* -- that's what professional mathematicians see.
"""

# ╔═╡ 1856ddae-7a78-11eb-3422-298e1103275b
md"""
What defines a linear transformation?  There are a few equivalent ways of giving a definition.
"""

# ╔═╡ 4b4fe818-7a78-11eb-2986-59e60063d346
md"""
**Linear transformation definitions:**
"""

# ╔═╡ 5d656494-7a78-11eb-12e8-d17856bd8c4d
md"""
- The intuitive definition:

   > The rectangles (gridlines) in the transformed image [above](#e0b657ce-7a03-11eb-1f9d-f32168cb5394) always become a lattice of congruent parallelograms.

- The easy operational (but devoid of intuition) definition:

   > A transformation is linear if it is defined by $v \mapsto A*v$ (matrix times vector) for some fixed matrix $A$.

- The scaling and adding definition:

   > 1. If you scale and then transform or if you transform and then scale, the result is always the same:
   >
   > $T(cv)=c \, T(v)$ ( $v$ is any vector, and $c$ any number.)
   >
   > 2. If you add and then transform or vice versa the result is the same:
   >
   > $T(v_1+v_2) = T(v_1) + T(v_2).$ ($v_1,v_2$ are any vectors.)

- The mathematician's definition:

   > (A consolidation of the above definition.) $T$ is linear if
   >
   > $T(c_1 v_1 + c_2 v_2) = c_1 T(v_1) + c_2 T(v_2)$ for all numbers $c_1,c_2$ and vectors $v_1,v_2$.  (This can be extended to beyond 2 terms.)

"""

# ╔═╡ b0e6d1ac-7a7d-11eb-0a9e-1310dcb5957f
md"""
### The matrix
"""

# ╔═╡ 7e4ad37c-7a84-11eb-1490-25090e133a7c
Resource("https://upload.wikimedia.org/wikipedia/en/c/c1/The_Matrix_Poster.jpg")

# ╔═╡ 96f47252-7a84-11eb-3d18-e3ba79dd20c2
md"""
No not that matrix!
"""

# ╔═╡ ae5b3a32-7a84-11eb-04c0-337a74105a58
md"""
The matrix for a linear transformation $T$ is easy to write down: The first
column is $T([1, 0])$ and the second is $T([0,1])$. That's it!
"""

# ╔═╡ c9f2b61e-7a84-11eb-3841-33739a226ff9
md"""
Once we have those, the linearity relation

$$T([x,y]) = x \, T([1,0]) + y \, T([0,1]) = x \, \mathrm{(column\ 1)} + y \, \mathrm{(column\ 2)}$$

is exactly the definition of matrix times vector. Try it.
"""

# ╔═╡ 23d8a45c-7a85-11eb-3a68-ef11e6f58cac
md"""
### Matrix multiply:  You know how to do it, but  why?
"""

# ╔═╡ 4a96d516-7a85-11eb-181c-63a6b461790b
md"""
Did you ever ask yourself why matrix multiply has that somewhat complicated multiplying and adding going on?
"""

# ╔═╡ 8206e1ee-7a8a-11eb-1f26-054f6b100076
let
	 A = randn(2,2)
	 B = randn(2,2)
	 v = rand(2)

	(lin(A) ∘ lin(B))(v) , lin(A * B)(v)
end

# ╔═╡ 7d803684-7a8a-11eb-33d2-89d5e2a05bcf
md"""
**Important:** The composition of the linear transformation is the linear transformation of the multiplied matrices!  There is only one definition of
matmul (matrix multiply) that realizes this fact.

To see what it is exactly, remember the first column of `lin(A) ∘ lin(B)` should
be the result of computing the two matrix times vectors  $y=A*[1,0]$ then $z=A*y$,
and the second column is the same for $[0,1]$.

This is worth writing out if you have never done this.
"""

# ╔═╡ 17281256-7aa5-11eb-3144-b72777334326
md"""
Let's try doing that with random matrices:
"""

# ╔═╡ 05049fa0-7a8e-11eb-283b-cb4753c4aaf0
begin

	P = randn(2, 2)
	Q = randn(2, 2)

	T₁ = lin(P) ∘ lin(Q)
	T₂ = lin(P*Q)

	lin(P*Q)((1, 0)), (lin(P)∘lin(Q))((1, 0))
end

# ╔═╡ 350f40f7-795f-4f33-89b8-ff9ba4819e1c
test_img = load(download(corgis_url));

# ╔═╡ 313cdcbd-5b11-41c8-9fcd-5aeaca3b8d24
test_pixels = 300;

# ╔═╡ 57848b42-7a8f-11eb-023a-cf247cb53819
md"""
`lin(P*Q)`
"""

# ╔═╡ da73d9f6-7a8d-11eb-2e6f-1b819bbb0185
begin
		[
			begin
			 x,y = transform_ij_to_xy(i,j, test_pixels)
			 X,Y =  T₁([x,y])
			 i,j = transform_xy_to_ij(test_img,X,Y)
			 getpixel(test_img,i,j)
			end

			for i = 1:test_pixels, j = 1:test_pixels
		]
end

# ╔═╡ 620ee7d8-7a8f-11eb-3888-356c27a2d591
md"""
`lin(P)∘lin(Q)`
"""

# ╔═╡ 30f522a0-7a8e-11eb-2181-8313760778ef
begin
		[
			begin
			 x,y = transform_ij_to_xy(i,j, test_pixels)
			 X,Y =  T₂([x,y])
			 i,j = transform_xy_to_ij(test_img,X,Y)
			 getpixel(test_img,i,j)
			end

			for i = 1:test_pixels, j = 1:test_pixels
		]
end

# ╔═╡ 04da7710-7a91-11eb-02a1-0b6e889150a2
md"""
# Coordinate transformations vs object transformations
"""

# ╔═╡ c2e0e032-7c4c-11eb-2b2a-27fe69c42a01
img;

# ╔═╡ c662e3d8-7c4c-11eb-0dcf-f9da2bd14baf
size(img)

# ╔═╡ d0e9a1e8-7c4c-11eb-056c-aff283c49c31
img[50,56]

# ╔═╡ 155cd218-7a91-11eb-0b4c-bd028507e925
md"""
If you want to move an object to the right, the first thing you might think of is adding 1 to the $x$ coordinate of every point.  The other thing you could do is to subtract one from the first coordinate of the coordinate system.  The latter is an example of a coordinate transform.
"""

# ╔═╡ fd25da12-7a92-11eb-20c0-995e7c46b3bc
md"""
### Coordinate transform of an array $(i, j)$ vs points $(x, y)$
"""

# ╔═╡ 1ab2265e-7c1d-11eb-26df-39c4c7289243
md"""
The original image has (1,1) in the upper left corner as an array but is thought
of as existing in the entire plane.
"""

# ╔═╡ 7c68c7b6-7a9e-11eb-3f7f-99bb10aedd95
Resource("https://raw.githubusercontent.com/mitmath/18S191/Spring21/notebooks/week3/coord_transform.png")

# ╔═╡ 7d0096ad-d89a-4ade-9679-6ee95f7d2044
begin
	function transform_xy_to_ij(img::AbstractMatrix, x::Float64, y::Float64)
	# convert coordinate system xy to ij
	# center image, and use "white" when out of the boundary

		rows, cols = size(img)
		m = max(cols, rows)

	    # function to take xy to ij
		xy_to_ij =  translate(rows/2, cols/2) ∘ swap ∘ flipy ∘ scale(m/2)

		# apply the function and "snap to grid"
		i,j = floor.(Int, xy_to_ij((x, y)))

	end

	function getpixel(img,i::Int,j::Int; circular::Bool=false, r::Real=200)
		#  grab image color or place default
		rows,cols = size(img)
		m = max(cols,rows)
		if circular
			c = (i-rows/2)^2 + (j-cols/2)^2 ≤ r*m^2/4
		else
			c = true
		end

		if 1 < i ≤ rows && 1 < j ≤ cols && c
			img[i, j]
		else
			#white(img[1, 1])
			black(img[1,1])
		end

	end

	# function getpixel(img,x::Float64,y::Float64)
	# 	i,j = transform_xy_to_ij(img,x,y)
	# 	getpixel(img,i,j)
	# end

	function transform_ij_to_xy(i::Int,j::Int,pixels)

	   ij_to_xy =  scale(2/pixels) ∘ flipy ∘ swap ∘ translate(-pixels/2,-pixels/2)
	   ij_to_xy([i,j])
	end

end

# ╔═╡ bf1954d6-7e9a-11eb-216d-010bd761e470
transform_ij_to_xy(1,1,400)

# ╔═╡ c1efc54a-7e9b-11eb-1e76-dbd0a66184a9
translate(-400,400)([1,1])

# ╔═╡ db4bc328-76bb-11eb-28dc-eb9df8892d01
md"""
# Inverses
"""

# ╔═╡ 0b8ed36c-7a1e-11eb-053c-63cf9ee0b16f
md"""
If $f$ is a function from 2-vectors to 2-vectors (say), we define the **inverse** of $f$, denoted
$f^{-1}$, to have the property that it "*undoes*" the effect of $f$, i.e.

$$f(f^{-1}(v))=v$$

and $f^{-1}(f(v))=v$.

This equation might be true for all $v$ or for some $v$ in a region.
"""

# ╔═╡ 7a4e785e-7a71-11eb-07fb-cfba453a117b
md"""
## Example: Scaling up and down
"""

# ╔═╡ 9264508a-7a71-11eb-1b7c-bf6e62788115
let
	v = rand(2)
	T = rotate(30)∘rotate(-30 )
	T(v) ,  v
end

# ╔═╡ e89339b2-7a71-11eb-0f97-971b2ed277d1
let
	  T = scale(0.5) ∘ scale(2)

	  v = rand(2)
	  T(v) .≈ v
end

# ╔═╡ 0957fd9a-7a72-11eb-0566-e93ef32fb626
md"""
We observe numerically that `scale(2)` and `scale(.5)` are mutually inverse transformations, i.e. each is the inverse of the other.
"""

# ╔═╡ c7cc412c-7aa5-11eb-2df1-d3d788047238
md"""
## Inverses: Solving equations
"""

# ╔═╡ ce620b8e-7aa5-11eb-370b-11e34b07d54d
md"""
What does an inverse really do?

Let's think about scaling again.
Suppose we scale an input vector $\mathbf{x}$ by 2 to get an output vector $\mathbf{x}$:

$$\mathbf{y} = 2 \mathbf{x}$$

Now suppose that you want to go backwards. If you are given $\mathbf{y}$, how do you find $\mathbf{x}$? In this particular case we see that $\mathbf{x} = \frac{1}{2} \mathbf{y}$.

If we have a *linear* transformation, we can write

$$\mathbf{y} = A \, \mathbf{x}$$

with a matrix $A$.

If we are given $\mathbf{y}$ and want to go backwards to find the $\mathbf{x}$ from that, we need to *solve a system of linear equations*.

*Usually*, but *not always*, we can solve these equations to find a new matrix $B$ such that

$$\mathbf{x} = B \, \mathbf{y},$$

i.e. $B$ *undoes* the effect of $A$. Then we have

$$\mathbf{x} = (B \, A) * \mathbf{x},$$

so that $B * A$ must be the identity matrix. We call $B$ the *matrix inverse* of $A$, and write

$$B = A^{-1}.$$

For $2 \times 2$ matrices we can write down an explicit formula for the matrix inverse, but in general we will need a computer to run an algorithm to find the inverse.

"""

# ╔═╡ 4f51931c-7aac-11eb-13ba-4b8768ac376f
md"""
### Inverting Linear Transformations
"""

# ╔═╡ 5ce799f4-7aac-11eb-0629-ebd8a404e9d3
let
	v = rand(2)
	A = randn(2,2)
    (lin(inv(A)) ∘ lin(A))(v) , v
end

# ╔═╡ 9b456686-7aac-11eb-3aa5-25e6c3c86aff
let
	 A = randn(2,2)
	 B = randn(2,2)
	 inv(A*B) ≈ inv(B) * inv(A)
end

# ╔═╡ c2b0a488-7aac-11eb-1d8b-edd6bd23d1fd
md"""
``A^{-1}
=
\begin{pmatrix} d & -b \\ -c & a  \end{pmatrix} / (ad-bc) \quad
``
if
``\ A \ =
\begin{pmatrix} a & b \\ c & d  \end{pmatrix} .
``
"""

# ╔═╡ 02d6b440-7aa7-11eb-1be0-b78dea91387f
md"""
### Inverting nonlinear transformations
"""

# ╔═╡ 0be9fb1e-7aa7-11eb-0116-c3e86ab82c77
md"""
What about if we have a *nonlinear* transformation $T$ -- can we invert it? In other words, if $\mathbf{y} = T(\mathbf{x})$, can we solve this to find $\mathbf{x}$ in terms of $\mathbf{y}$?

In general this is a difficult question! Sometimes we can do so analytically, but usually we cannot.

Nonetheless, there are *numerical* methods that can sometimes solve these equations, for example the [Newton method](https://en.wikipedia.org/wiki/Newton%27s_method).

There are several implementations of such methods in Julia, e.g. in the [NonlinearSolve.jl package](https://github.com/JuliaComputing/NonlinearSolve.jl). We have used that to write a function `inverse` that tries to invert nonlinear transformations of our images.
"""

# ╔═╡ 7609d686-7aa7-11eb-310a-3550509504a1
md"""
# The Big Diagram of Transforming Images
"""

# ╔═╡ 1b9faf64-7aab-11eb-1396-6fb89be7c445
Resource("https://raw.githubusercontent.com/mitmath/18S191/Spring21/notebooks/week3/comm2.png")

# ╔═╡ 5f0568dc-7aad-11eb-162f-0d6e26f17d59
md"""
Note that we are defining the map with the inverse of T so we can go pixel by pixel in the result.
"""

# ╔═╡ 8d32fff4-7c1b-11eb-1fa1-6ff2d87bfb73
md"""
## Collisions
"""

# ╔═╡ 80456168-7c1b-11eb-271c-83ef59a41102
Resource("https://raw.githubusercontent.com/mitmath/18S191/Spring21/notebooks/week3/collide.png")

# ╔═╡ 62a9201c-7938-11eb-144c-15690c06be94
begin
	function inverse(f, y, u0=@SVector[0.0, 0.0])
	    prob = NonlinearProblem{false}( (u, p) -> f(u, p) .- y, u0)
	    solver = solve(prob, NewtonRaphson(), tol = 1e-4)
	    return solver.u
	end

	inverse(f) = y -> inverse( (u, p) -> f(SVector(u...)), y )
end

# ╔═╡ 5227afd0-7641-11eb-0065-918cb8538d55
md"""

Check out
[Linear Map Wikipedia](https://en.wikipedia.org/wiki/Linear_map)

[Transformation Matrix Wikipedia](https://en.wikipedia.org/wiki/Transformation_matrix)
"""

# ╔═╡ 4c93d784-763d-11eb-1f48-81d4d45d5ce0
md"""
## Why are we doing this backwards?

If one moves the colors forward rather than backwards you have trouble dealing
with the discrete pixels.  You may have gaps.  You may have multiple colors going
to the same pixel.

An interpolation scheme or a newton scheme could work for going forwards, but very likely care would be neeeded for a satisfying general result.
"""

# ╔═╡ c536dafb-4206-4689-ad6d-6935385d8fdf
md"""
# Appendix
"""

# ╔═╡ fb509fb4-9608-421d-9c40-a4375f459b3f
det_A = det(A)

# ╔═╡ 40655bcc-6d1e-4d1e-9726-41eab98d8472
img_sources = [
	"https://user-images.githubusercontent.com/6933510/108605549-fb28e180-73b4-11eb-8520-7e29db0cc965.png" => "Corgis",
	"https://images.squarespace-cdn.com/content/v1/5cb62a904d546e33119fa495/1589302981165-HHQ2A4JI07C43294HVPD/ke17ZwdGBToddI8pDm48kA7bHnZXCqgRu4g0_U7hbNpZw-zPPgdn4jUwVcJE1ZvWQUxwkmyExglNqGp0IvTJZamWLI2zvYWH8K3-s_4yszcp2ryTI0HqTOaaUohrI8PISCdr-3EAHMyS8K84wLA7X0UZoBreocI4zSJRMe1GOxcKMshLAGzx4R3EDFOm1kBS/fluffy+corgi?format=2500w" => "Long Corgi",
"https://previews.123rf.com/images/camptoloma/camptoloma2002/camptoloma200200020/140962183-pembroke-welsh-corgi-portrait-sitting-gray-background.jpg"=>"Portrait Corgi",
	"https://www.eaieducation.com/images/products/506618_L.jpg"=>"Graph Paper"
]

# ╔═╡ 55898e88-36a0-4f49-897f-e0850bd2b0df
img = if show_grid
	with_gridlines(img_original;n=ngrid)
else
	img_original
end;

# ╔═╡ b754bae2-762f-11eb-1c6a-01251495a9bb
begin
	white(c::RGB) = RGB(1,1,1)
	white(c::RGBA) = RGBA(1,1,1,0.75)
	black(c::RGB) = RGB(0,0,0)
	black(c::RGBA) = RGBA(0,0,0,0.75)
end

# ╔═╡ 83d45d42-7406-11eb-2a9c-e75efe62b12c
function with_gridlines(img::Array{<:Any,2}; n = 10)
    n = 2n+1
	rows, cols = size(img)
	result = copy(img)
	# stroke = zero(eltype(img))#RGBA(RGB(1,1,1), 0.75)

	stroke = RGBA(1, 1, 1, 0.75)

	result[ floor.(Int,LinRange(1, rows, n) ), : ] .= stroke
	#result[ ceil.(Int,LinRange(1, rows, n) ), : ] .= stroke
	result[ : , floor.(Int,LinRange(1, cols, n) )] .= stroke
	#result[ : , ceil.(Int,LinRange(1, cols, n) )] .= stroke

    result[  rows ÷2    , :] .= RGBA(0,1,0,1)
	#result[  1+rows ÷2    , :] .= RGBA(0,1,0,1)
	result[ : ,  cols ÷2   ,] .= RGBA(1,0,0,1)
	#result[ : ,  1 + cols ÷2   ,] .= RGBA(1,0,0,1)
	return result
end

# ╔═╡ 2e8c4a48-d535-44ac-a1f1-4cb26c4aece6
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ Cell order:
# ╟─972b2230-7634-11eb-028d-df7fc722ec70
# ╟─b7895bd2-7634-11eb-211e-ef876d23bd88
# ╟─230b0118-30b7-4035-ad31-520165a76fcc
# ╠═6b473b2d-4326-46b4-af38-07b61de287fc
# ╟─230cba36-9d0a-4726-9e55-7df2c6743968
# ╠═96766502-7a06-11eb-00cc-29849773dbcf
# ╟─890d30b9-2cd0-4d3a-99f6-f7d3d7858fda
# ╟─85fba8fb-a9ea-444d-831b-ec6489b58b4f
# ╟─06beabc3-2aa7-4e78-9bae-dc4b37251aa2
# ╟─26dd0e98-7a75-11eb-2196-5d7bda201b19
# ╟─e0b657ce-7a03-11eb-1f9d-f32168cb5394
# ╟─005ca75a-7622-11eb-2ba4-9f450e71df1f
# ╟─23ade8ee-7a09-11eb-0e40-296c6b831d74
# ╠═58a30e54-7a08-11eb-1c57-dfef0000255f
# ╟─2efaa336-7630-11eb-0c17-a7d4a0141dac
# ╟─7f28ac40-7914-11eb-1403-b7bec34aeb94
# ╟─ce55beee-7643-11eb-04bc-b517703facff
# ╠═b76a5bd6-802f-11eb-0951-1f1092dee8de
# ╟─5d33f6ea-7e9c-11eb-2fb3-dbb7cb07c60c
# ╟─45dccdec-7912-11eb-01b4-a97e30344f39
# ╟─d2fb356e-7f32-11eb-177d-4f47d6c9e59b
# ╠═ca28189e-7e9a-11eb-21d6-bd819f3e0d3a
# ╠═ccea7244-7f2f-11eb-1b7b-b9b8473a8c74
# ╟─55b5fc92-7a76-11eb-3fba-854c65eb87f9
# ╟─85686412-7a75-11eb-3d83-9f2f8a3c5509
# ╟─a7df7346-79f8-11eb-1de6-71f027c46643
# ╟─044e6128-79fe-11eb-18c1-395ae857dc73
# ╟─78d61e28-79f9-11eb-0605-e77d206cda84
# ╟─aad4d6e4-79f9-11eb-0342-b900a41cfbaf
# ╟─d42aec08-76ad-11eb-361a-a1f2c90fd4ec
# ╟─d9115c1a-7aa0-11eb-38e4-d977c5a6b75b
# ╟─e965cf5e-79fd-11eb-201d-695b54d08e54
# ╟─1e11c1ec-79fe-11eb-1867-9da72b3f3bc4
# ╟─28ef451c-7aa1-11eb-340c-ab3a1193a3c4
# ╟─a0afe3ae-76b9-11eb-2301-cde7260ddd7f
# ╟─fc2deb7c-7aa1-11eb-019f-d3e3c80b9ff1
# ╠═d364f91a-76b9-11eb-1807-75e733940d53
# ╟─080d87e0-7aa2-11eb-18f5-2fb6a7a5bcb4
# ╠═15283aba-7aa2-11eb-389c-e9f215bd03e2
# ╟─2612d2c2-7aa2-11eb-085a-1f27b6174995
# ╟─a290d5e2-7a02-11eb-37db-41bf86b1f3b3
# ╠═b4cdd412-7a02-11eb-149a-df1888a0f465
# ╟─704a87ec-7a1e-11eb-3964-e102357a4d1f
# ╠═4b0e8742-7a70-11eb-1e78-813f6ad005f4
# ╟─44792484-7a20-11eb-1c09-95b27b08bd34
# ╟─f650b788-7a70-11eb-0b20-779d2f18f111
# ╟─c852d398-7aa2-11eb-2ded-ab2e5236e9b2
# ╟─061076c2-7aa3-11eb-0d04-b7cbc60e6cb2
# ╟─014c14a6-7a72-11eb-119b-f5cfc82085ca
# ╟─89f0bc54-76bb-11eb-271b-3190b4d8cbc0
# ╟─f70f7ea8-76b9-11eb-3bd7-87d40a2861b1
# ╟─bf28c388-76bd-11eb-08a7-af2671218017
# ╠═5655d2a6-76bd-11eb-3042-5b2dd3f6f44e
# ╠═56f1e4cc-7a03-11eb-187b-c5a917978eb9
# ╠═70dc4346-7a03-11eb-055e-111d2519a44c
# ╠═852592d6-76bd-11eb-1265-5f200e39113d
# ╠═8e36f4a2-76bd-11eb-2fda-9d1424752812
# ╟─a8bf7128-7aa5-11eb-3ee9-953b0b5ccd01
# ╟─ad700740-7a74-11eb-3369-15e5fd89194d
# ╟─e051259a-7a74-11eb-12fc-99c5dc867fbd
# ╟─1856ddae-7a78-11eb-3422-298e1103275b
# ╟─4b4fe818-7a78-11eb-2986-59e60063d346
# ╟─5d656494-7a78-11eb-12e8-d17856bd8c4d
# ╟─b0e6d1ac-7a7d-11eb-0a9e-1310dcb5957f
# ╟─7e4ad37c-7a84-11eb-1490-25090e133a7c
# ╟─96f47252-7a84-11eb-3d18-e3ba79dd20c2
# ╟─ae5b3a32-7a84-11eb-04c0-337a74105a58
# ╟─c9f2b61e-7a84-11eb-3841-33739a226ff9
# ╟─23d8a45c-7a85-11eb-3a68-ef11e6f58cac
# ╟─4a96d516-7a85-11eb-181c-63a6b461790b
# ╠═8206e1ee-7a8a-11eb-1f26-054f6b100076
# ╟─7d803684-7a8a-11eb-33d2-89d5e2a05bcf
# ╟─17281256-7aa5-11eb-3144-b72777334326
# ╠═05049fa0-7a8e-11eb-283b-cb4753c4aaf0
# ╠═350f40f7-795f-4f33-89b8-ff9ba4819e1c
# ╠═313cdcbd-5b11-41c8-9fcd-5aeaca3b8d24
# ╟─57848b42-7a8f-11eb-023a-cf247cb53819
# ╟─da73d9f6-7a8d-11eb-2e6f-1b819bbb0185
# ╟─620ee7d8-7a8f-11eb-3888-356c27a2d591
# ╟─30f522a0-7a8e-11eb-2181-8313760778ef
# ╟─04da7710-7a91-11eb-02a1-0b6e889150a2
# ╠═c2e0e032-7c4c-11eb-2b2a-27fe69c42a01
# ╠═c662e3d8-7c4c-11eb-0dcf-f9da2bd14baf
# ╠═d0e9a1e8-7c4c-11eb-056c-aff283c49c31
# ╟─155cd218-7a91-11eb-0b4c-bd028507e925
# ╟─fd25da12-7a92-11eb-20c0-995e7c46b3bc
# ╟─1ab2265e-7c1d-11eb-26df-39c4c7289243
# ╟─7c68c7b6-7a9e-11eb-3f7f-99bb10aedd95
# ╠═7d0096ad-d89a-4ade-9679-6ee95f7d2044
# ╠═bf1954d6-7e9a-11eb-216d-010bd761e470
# ╠═c1efc54a-7e9b-11eb-1e76-dbd0a66184a9
# ╟─db4bc328-76bb-11eb-28dc-eb9df8892d01
# ╟─0b8ed36c-7a1e-11eb-053c-63cf9ee0b16f
# ╟─7a4e785e-7a71-11eb-07fb-cfba453a117b
# ╠═9264508a-7a71-11eb-1b7c-bf6e62788115
# ╠═e89339b2-7a71-11eb-0f97-971b2ed277d1
# ╟─0957fd9a-7a72-11eb-0566-e93ef32fb626
# ╟─c7cc412c-7aa5-11eb-2df1-d3d788047238
# ╟─ce620b8e-7aa5-11eb-370b-11e34b07d54d
# ╟─4f51931c-7aac-11eb-13ba-4b8768ac376f
# ╠═5ce799f4-7aac-11eb-0629-ebd8a404e9d3
# ╠═9b456686-7aac-11eb-3aa5-25e6c3c86aff
# ╟─c2b0a488-7aac-11eb-1d8b-edd6bd23d1fd
# ╟─02d6b440-7aa7-11eb-1be0-b78dea91387f
# ╟─0be9fb1e-7aa7-11eb-0116-c3e86ab82c77
# ╟─7609d686-7aa7-11eb-310a-3550509504a1
# ╟─1b9faf64-7aab-11eb-1396-6fb89be7c445
# ╟─5f0568dc-7aad-11eb-162f-0d6e26f17d59
# ╟─8d32fff4-7c1b-11eb-1fa1-6ff2d87bfb73
# ╟─80456168-7c1b-11eb-271c-83ef59a41102
# ╠═62a9201c-7938-11eb-144c-15690c06be94
# ╟─5227afd0-7641-11eb-0065-918cb8538d55
# ╟─4c93d784-763d-11eb-1f48-81d4d45d5ce0
# ╟─c536dafb-4206-4689-ad6d-6935385d8fdf
# ╟─fb509fb4-9608-421d-9c40-a4375f459b3f
# ╠═40655bcc-6d1e-4d1e-9726-41eab98d8472
# ╠═55898e88-36a0-4f49-897f-e0850bd2b0df
# ╠═b754bae2-762f-11eb-1c6a-01251495a9bb
# ╠═83d45d42-7406-11eb-2a9c-e75efe62b12c
# ╟─2e8c4a48-d535-44ac-a1f1-4cb26c4aece6