text:
type AminoAcid
    name
    length
    content
end

type Alignment
    seqX # string
    seqY # string
    startPos   # (x, y)
    endPos
    score
end

type Group
    startCoord
    endCoord
    members
    memberNum
end

aminoAcids = Set()
push!(aminoAcids,
"A",
"R",
"N",
"D",
"B",
"C",
"E",
"Q",
"Z",
"G",
"H",
"I",
"L",
"K",
"M",
"F",
"P",
"S",
"T",
"W",
"Y",
"V")

function loadData(dataPath)
    data = open(dataPath)
    lines = readlines(data)

    aminoAcidlist = []

    for i = 1:2:length(lines)
        lineEles = split(lines[i], ",")
        contend = split(lines[1+i], "")

        for j = contend
            if !checkElement(j, aminoAcids)
                print("Wrong element: ")
                print(j)
                print("\n")
                return
            end
        end

        aminoAcid = AminoAcid(
            split(lineEles[1])[2],
            split(lineEles[2])[2],
            split(lines[1+i], "")
        )
        push!(aminoAcidlist, aminoAcid)
    end

    return aminoAcidlist
end

function checkElement(target, ele)
    for i = ele
        if i == target
            return true
        end
    end
    return false
end

#=
A   ALA
R   ARG
N   ASN
D   ASP
C   CYS
Q   GLN
E   GLU
G   GLY
H   HIS
I   ILE
L   LEU
K   LYS
M   MET
F   PHE
P   PRO
S   SER
T   THR
W   TRP
Y   TYR
V   VAL

A	1
R	2
N	3
D	4
C	5
Q	6
E	7
G	8
H	9
I	10
L	11
K	12
M	13
F	14
P	15
S	16
T	17
W	18
Y	19
V	20

4	-1	-2	-2	0	-1	-1	0	-2	-1	-1	-1	-1	-2	-1	1	0	-3	-2	0
-1	5	0	-2	-3	1	0	-2	0	-3	-2	2	-1	-3	-2	-1	-1	-3	-2	-3
-2	0	6	1	-3	0	0	0	1	-3	-3	0	-2	-3	-2	1	0	-4	-2	-3
-2	-2	1	6	-3	0	2	-1	-1	-3	-4	-1	-3	-3	-1	0	-1	-4	-3	-3
0	-3	-3	-3	9	-3	-4	-3	-3	-1	-1	-3	-1	-2	-3	-1	-1	-2	-2	-1
-1	1	0	0	-3	5	2	-2	0	-3	-2	1	0	-3	-1	0	-1	-2	-1	-2
-1	0	0	2	-4	2	5	-2	0	-3	-3	1	-2	-3	-1	0	-1	-3	-2	-2
0	-2	0	-1	-3	-2	-2	6	-2	-4	-4	-2	-3	-3	-2	0	-2	-2	-3	-3
-2	0	1	-1	-3	0	0	-2	8	-3	-3	-1	-2	-1	-2	-1	-2	-2	2	-3
-1	-3	-3	-3	-1	-3	-3	-4	-3	4	2	-3	1	0	-3	-2	-1	-3	-1	3
-1	-2	-3	-4	-1	-2	-3	-4	-3	2	4	-2	2	0	-3	-2	-1	-2	-1	1
-1	2	0	-1	-3	1	1	-2	-1	-3	-2	5	-1	-3	-1	0	-1	-3	-2	-2
-1	-1	-2	-3	-1	0	-2	-3	-2	1	2	-1	5	0	-2	-1	-1	-1	-1	-1
-2	-3	-3	-3	-2	-3	-3	-3	-1	0	0	-3	0	6	-4	-2	-2	1	3	1
-1	-2	-2	-1	-3	-1	-1	-2	-2	-3	-3	-1	-2	-4	7	-1	-1	-4	-3	-2
1	-1	1	0	-1	0	0	0	-1	-2	-2	0	-1	-2	-1	4	1	-3	-2	-2
0	-1	0	-1	-1	-1	-1	-2	-2	-1	-1	-1	-1	-2	-1	1	5	-2	-2	0
-3	-3	-4	-4	-2	-2	-3	-2	-2	-3	-2	-3	-1	1	-4	-3	-2	11	2	-3
-2	-2	-2	-3	-2	-1	-2	-3	2	-1	-1	-2	-1	3	-3	-2	-2	2	7	-1
0	-3	-3	-3	-1	-2	-2	-3	-3	3	1	-2	-1	1	-2	-2	0	-3	-1	4

=#

using DataFrames

function getBlosum()

    datapath = joinpath(@__DIR__, "..", "materials")
    blosumpath = joinpath(datapath, "ezBlosum.csv")
    #blosum = readtable("../materials/ezBlosum.csv", separator = ',')
    blosum = readtable(blosumpath, separator = ',')
    return blosum
end

blosum = getBlosum()

blosumIndex = Dict(
    "A" => 1,
    "R"	=> 2,
    "N"	=> 3,
    "D"	=> 4,
    "C"	=> 5,
    "Q"	=> 6,
    "E"	=> 7,
    "G"	=> 8,
    "H"	=> 9,
    "I"	=> 10,
    "L"	=> 11,
    "K"	=> 12,
    "M"	=> 13,
    "F"	=> 14,
    "P"	=> 15,
    "S"	=> 16,
    "T"	=> 17,
    "W"	=> 18,
    "Y"	=> 19,
    "V"	=> 20
)

function getMappingScore(col, row, blosum)
    col = uppercase(col)
    row = uppercase(row)

    col = blosumIndex[col]
    row = blosumIndex[row]

    return blosum[col][row]
end

#=
This function searches for the maximum score left in the score matrix. If
it is >= the score threshold, then it will report the alignment corresponding
to it and return true, else false is returned.
@param matrix Score matrix
@param path Path matrix
@param xLen Length of seqX + 1
@param yLen Length of seqY + 1
@param alignmentList List of alignments to which we will append new alignment
@param threshold Lowest allowable score
@return True if alignment has sufficient score and is successfully added to
        list, else false
=#
function getAlignment(matrix, path, xLen, yLen, alignmentList, threshold, width, xAminoAcid, yAminoAcid)
    maxScore = 0
    col = 0
    row = 0
    ## First get the maxScore and its coordinates
    for r = 2:yLen
        for c = 2:xLen
            if matrix[r, c] > maxScore
                maxScore = matrix[r, c]
                col = c
                row = r
            end
        end
    end
    ## Check if maxScore >= threshold
    if maxScore < threshold
        return false ## score too low
    end
    ## Now retrieve the path
    seqX = ""
    seqY = ""
    x1 = 0
    x2 = col
    lastX = col
    y1 = 0
    y2 = row
    lastY = row

    #zero out the region around the maxScore
    for w in 0:width
        if row - w >= 2 && col - w >= 2
            matrix[row - w, col - w] = 0
        end
        if row - w >= 2 && col + w <= xLen
            matrix[row - w, col + w] = 0
        end
        if row + w <= yLen && col - w >= 2
            matrix[row + w, col - w] = 0
        end
        if row + w <= yLen && col + w <= xLen
            matrix[row + w, col + w] = 0
        end
    end
    pathVal = path[row, col] ## init path location
    while pathVal != 0
        ## First zero out the score element and nearby left/right region
        matrix[row, col] = 0
        for w in 1:width
            if col - w >= 2
                matrix[row, col - w] = 0
            end
            if col + w <= xLen
                matrix[row, col + w] = 0
            end
        end
        lastX = row
        lastY = col
        if pathVal == 1  # match/mismatch
            row -= 1
            col -= 1
            seqX = string(xAminoAcid.content[col], seqX)
            seqY = string(yAminoAcid.content[row], seqY)
        elseif pathVal == 2
            row -= 1
            seqX = string("-", seqX)
            seqY = string(yAminoAcid.content[row], seqY)
        else ## pathVal == 3
            col -= 1
            seqX = string(xAminoAcid.content[col], seqX)
            seqY = string("-", seqY)
        end
        pathVal = path[row, col]
    end
    x1 = lastX
    y1 = lastY
    ## Add this alignment to the list of alignments
    alignment = Alignment(
        seqX,
        seqY,
        (x1, y1),
        (y2, x2),
        maxScore
    )
    push!(alignmentList, alignment)
    return true
end

#=
This function returns a filtered list of alignments. If alignments are
overlapping (i.e. have the same start position), then only the one with
the highest score will be added.
@param alignmentList List of alignments
@return filtered list of non-overlapping alignments
=#
function filterAlignments(alignmentList)
    # sort the alignments first
    sort!(alignmentList, by = x -> x.startPos)
    filtList = []
    numFilt = 1
    i = length(alignmentList) - 1
    # push / pop the very first alignment
    push!(filtList, pop!(alignmentList))
    while i > 0
        if alignmentList[i].startPos == filtList[numFilt].startPos
            if alignmentList[i].score > filtList[numFilt].score
                pop!(filtList)
                push!(filtList,pop!(alignmentList))
            else
                pop!(alignmentList)
            end
        else # new start position
            push!(filtList,pop!(alignmentList))
            numFilt += 1
        end
        i -= 1
    end
    return filtList
end

#=
 = alignmentList is a list of Alignment object
 =#

function mergeAlignment(alignmentList)
    seqXGroups = [] # list of group x object
    seqYGroups = [] # list of group y object
    for alignment = alignmentList
        pushX = false
        for xGroup in seqXGroups
            if isCoordSimilar((xGroup.startCoord, xGroup.endCoord), (alignment.startPos[2], alignment.endPos[2]))
                push!(xGroup.members, alignment)
                xGroup.startCoord = (xGroup.startCoord * xGroup.memberNum + alignment.startPos[2]) / (xGroup.memberNum + 1)
                xGroup.endCoord = (xGroup.endCoord * xGroup.memberNum + alignment.endPos[2]) / (xGroup.memberNum + 1)
                xGroup.memberNum += 1
                pushX = true
            end
        end

        if pushX == false
            newSet = Set()
            push!(newSet, alignment)
            push!(
                seqXGroups,
                Group(
                    alignment.startPos[2],
                    alignment.endPos[2],
                    newSet,
                    1
                )
            )
        end

        pushY = false
        for yGroup in seqYGroups
            if isCoordSimilar((yGroup.startCoord, yGroup.endCoord), (alignment.startPos[1], alignment.endPos[1]))
                push!(yGroup.members, alignment)
                yGroup.startCoord = (yGroup.startCoord * yGroup.memberNum + alignment.startPos[1]) / (yGroup.memberNum + 1)
                yGroup.endCoord = (yGroup.endCoord * yGroup.memberNum + alignment.endPos[1]) / (yGroup.memberNum + 1)
                yGroup.memberNum += 1
                pushY = true
            end
        end

        if pushY == false
            newSet = Set()
            push!(newSet, alignment)
            push!(
                seqYGroups,
                Group(
                    alignment.startPos[1],
                    alignment.endPos[1],
                    newSet,
                    1
                )
            )
        end
    end
    return (seqXGroups, seqYGroups)
end

#=
    xCoord1 = (x1, x2) / (y1, y2)
    xCoord2 = (x1, x2) / (y1, y2)

    e.g.
    |xCoord1(x1) - xCoord2(x1)| <= distance (usually be 1)
    |xCoord1(x2) - xCoord2(x2)| <= distance (usually be 1)
=#
function isCoordSimilar(xCoord1, xCoord2)
    distance = 2
    if (abs(xCoord1[1] - xCoord2[1]) <= distance) && (abs(xCoord1[2] - xCoord2[2]) <= distance)
        return true
    else
        return false
    end
end

max_stars_repo_name: yguoogo/Suboptimal-Local-Alignment

max_issues_repo_name: yguoogo/Suboptimal-Local-Alignment

max_forks_repo_name: yguoogo/Suboptimal-Local-Alignment


text:
mutable struct image_metadata_t <: LCMType
    key::String
    n::Int32
    value::Vector{UInt8}
end

@lcmtypesetup(image_metadata_t,
    value => (n,)
)

ext: jl


text:
# Functions for the rotation-correlation method

"""
    rotation_correlation(t1, t2, window_start=starttime(t1), window_end=endtime(t1); nphi=$SPLIT_NPHI, dt_max=$SPLIT_DT_MAX) -> dt, phi, xcorr_map, phi_best, dt_best

Compute the rotation-correlation map between two traces `t1` and `t2`

`phi` is a set of orientations to test, whilst `dt_max` is the largest delay
time tested.

`xcorr_map[idt,iphi]` holds the normalised cross correlation at `dt[idt]`
 and `phi[iphi]`.  `phi_best` and `dt_best` give the value of the fast
 orientation and delay time, respectively, at the maximum correlation point.
"""
function rotation_correlation(t1::Seis.Trace{T,V}, t2::Seis.Trace{T,V},
                              window_start=Seis.starttime(t1), window_end=Seis.endtime(t1);
                              nphi=SPLIT_NPHI, dt_max=SPLIT_DT_MAX) where {T,V}
    t1, t2 = check_trace_order(t1, t2) # t2 is now clockwise of t1
    n, e = Seis.rotate_through(t1, t2, -t1.sta.azi)
    Seis.cut!.((n, e), window_start, window_end)
    phi = range(-90., stop=90., length=nphi+1)[1:nphi]
    _rotation_correlation!(n, e, phi, dt_max)
end

"""
    _rotation_correlation!(n, e, phi, dt_max) -> dt, phi, xcorr_map, phi_best, dt_best

Compute the rotation-correlation map between the two traces `n` and `e`.
`n` must point northwards and `e` eastwards.

`phi` is a set of orientations to test, whilst `dt_max` is the largest delay
time tested.

`xcorr_map[idt,iphi]` holds the normalised cross correlation at `dt[idt]`
 and `phi[iphi]`.  `phi_best` and `dt_best` give the value of the fast
 orientation and delay time, respectively, at the maximum correlation point.
"""
function _rotation_correlation!(n::Seis.Trace{T,V}, e::Seis.Trace{T,V}, phi, dt_max) where {T,V}
    npts = Seis.nsamples(n)
    Seis.nsamples(e) == npts || throw(ArgumentError("traces must be the same length"))
    nphi = length(phi)
    nxcorr = 2npts - 1
    # Correlation map in dt and phi
    xcorr_map = Array{T}(undef, nxcorr, nphi)
    xcorr = Vector{T}(undef, nxcorr)
    # Times for each cross-correlation point
    time = ((1-npts):(npts-1)).*n.delta
    # Create work arrays and FFT plan for (I)FFTs.  This is just a tuple of stuff to
    # splat into the call to _xcorr_fd!()
    work_arrays = _correlation_arrays(T, npts)
    # Calculate correlation map
    ϕ_prev = 0.0
    @inbounds for (iϕ, ϕ) in enumerate(phi)
        δϕ = ϕ - ϕ_prev
        rotate_traces!(n, e, -δϕ)
        _xcorr_fd!(xcorr, Seis.trace(n), Seis.trace(e), work_arrays...; norm=true)
        xcorr_map[:,iϕ] .= xcorr
        ϕ_prev = ϕ
    end
    # Orient map to correct location
    idt, iphi = maxabsloc2(xcorr_map)
    dt_best = abs(time[idt])
    # Flip map if fast and slow anticorrelated
    xcorr_map[idt,iphi] < 0 && (xcorr_map .*= -1)
    xcorr_map = if time[idt] >= 0
        # Fast wave arrives later, so take positive time half
        # of map and then rotate map 90°
        iphi = mod1(iphi + nphi÷2, nphi)
        xcorr_map = xcorr_map[0 .<= time .<= dt_max, :]
        circshift(xcorr_map, (0, nphi÷2))
    else
        # Fast wave arrives first, so take negative time half
        # of map
        xcorr_map = reverse(xcorr_map[-dt_max .<= time .<= 0, :], dims=1)
    end
    phi_best = phi[iphi]
    time[0 .<= time .<= dt_max], phi, xcorr_map, phi_best, dt_best
end

"""
    _correlation_arrays(T, n) -> (apad, bpad, Apad, Bpad, conv, fplan, bplan)

Return a tuple of work arrays to be used be `_xcorr_fd!` to perform
in-place cross-correlation of two traces of length `n` and element
type `T`.
"""
function _correlation_arrays(T, n)
    nxc = 2n - 1
    np2 = nxc > 1024 ? nextprod([2,3,5], nxc) : nextpow(2, nxc)
    np = np2÷2 + 1
    # Input vectors padded to power-of-2 length
    apad = Vector{T}(undef, np2)
    bpad = similar(apad)
    # Fourier transformed vectors
    Apad = Vector{Complex{T}}(undef, np)
    Bpad = similar(Apad)
    # Convolution of two vectors
    conv = similar(Apad)
    # FFTW transform plans
    fplan = FFTW.plan_rfft(apad)
    bplan = FFTW.plan_irfft(conv, np2)
    (apad, bpad, Apad, Bpad, conv, fplan, bplan)
end

"""
    _xcorr_fd!(xc, a, b, apad, bpad, Apad, Bpad, conv, fplan, bplan; norm=false) -> xc

Compute the cross correlation between arrays `a` and `b`, which must
have the same length.  `xc` is filled in-place and returned.

Remaining arguments are work arrays and FFT forward and backward
real-to-complex plans computed with a call to `_correlation_arrays`.

If `norm` is `true`, then the cross correlation is normalised by
`sqrt(ACF(a) * ACF(b))`, where `ACF` is the autocorrelation.

For speed, no bounds checks are performed in this function and
therefore behaviour is undefined if any arrays are unexpected sizes.
"""
function _xcorr_fd!(xc, a, b, apad, bpad, Apad, Bpad, conv, fplan, bplan; norm=false)
    n = length(a)
    nxc = 2n - 1
    @inbounds apad[1:n] .= a
    @inbounds for i in 1:n
        bpad[i] = b[end-i+1]
    end
    @inbounds apad[(n+1):end] .= 0
    @inbounds bpad[(n+1):end] .= 0
    LinearAlgebra.mul!(Apad, fplan, apad)
    LinearAlgebra.mul!(Bpad, fplan, bpad)
    @inbounds conv .= Apad .* Bpad
    LinearAlgebra.mul!(apad, bplan, conv)
    @inbounds xc .= apad[1:nxc]
    if norm
        norm_val = sum(x->x^2, a) * sum(x->x^2, b)
        norm_val = √norm_val
        xc ./= norm_val
    end
    xc
end

"""
    cross_correlation_td!(xc, a, b) -> xc

Fill the contents of `xc` with the cross-correlation of the two vectors `a` and `b`,
which must be the same length `n`.  `xc` must have length `2n - 1`.  The zero-lag
point is therefore at `xc[n]`.

This is done in the time domain; rough benchmarking suggests that this can be faster
than a frequency-domain approach when `a` and `b` are less than ~350 in length.
"""
function cross_correlation_td!(xc, a, b)
    n = length(a)
    @boundscheck length(b) == n ||
        throw(ArgumentError("two vectors must be the same length"))
    @boundscheck length(xc) == 2n - 1 ||
        throw(ArgumentError("output vector must be `2n - 1` long"))
    xc .= 0
    i = 0
    @inbounds for ilag in (1-n):(n-1)
        i += 1
        iover = n - abs(ilag)
        if ilag <= 0
            @simd for j in 1:iover
                xc[i] += a[j]*b[j-ilag]
            end
        else
            @simd for j in 1:iover
                xc[i] += a[j+ilag]*b[j]
            end
        end
    end
    xc
end

"""
    cross_correlation_td(a, b) -> xc

Return the cross-correlation of the two vectors `a` and `b`, which must be the same
length.

This is done in the time domain; rough benchmarking suggests that this can be faster
than a frequency-domain approach when `a` and `b` are less than ~350 in length.
"""
cross_correlation_td(a, b) = cross_correlation!(similar(a, 2length(a)-1), a, b)

"""
    cross_correlation_fd!(xc, a, b) -> xc

Fill `xc` with the cross-correlation of the two vectors `a` and `b`, which must be
the same length `n`.  `xc` must have length `2n - 1`.  This is done in the frequency domain.
"""
cross_correlation_fd!(xc, a, b) = (xc .= DSP.xcorr(a, b); xc)

"""
    cross_correlation_fd(a, b) -> xc

Return the cross-correlation of the two vectors `a` and `b`, which must be the same
length.  This is done in the frequency domain.
"""
cross_correlation_fd(a, b) = DSP.xcorr(a, b)

"""
    cross_correlation!(xc, a, b, cutoff=350) -> xc

Fill `xc` with the cross-correlation between the vectors `a` and `b`, which must be
the same length `n`.  `xc` must have length `2n - 1`.

When `n` is less than `cutoff`, a time-domain implementation is used, whilst
when `n` is `cutoff` or larger, a frequency-domain implementation is used.
"""
cross_correlation!(xc, a, b, cutoff=350) = length(a) < cutoff ?
                                               cross_correlation_td!(xc, a, b) :
                                               cross_correlation_fd!(xc, a, b)

"""
    cross_correlation(a, b, cutoff=350) -> xc

Return the cross-correlation  `xc` between the vectors `a` and `b`, which must be the same
length `n` and have the same element type.  `xc` has length `2n - 1`.

When `n` is less than `cutoff`, a time-domain implementation is used, whilst
when `n` is `cutoff` or larger, a frequency-domain implementation is used.
"""
cross_correlation(a::AbstractVector{T}, b::AbstractVector{T}, cutoff=350) where T =
    length(a) < cutoff ?
        cross_correlation_td!(Vector{T}(undef, 2*length(a) - 1), a, b) :
        cross_correlation_fd(a, b)

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
####################
## TURNOVER TESTS ##
####################
# projectdir = "C:/Users/Arman Angaji/OneDrive - Universität zu Köln/Dokumente/Uni-Köln/Masterarbeit/Workspace/Julia_Master/MasterProject_Julia"
# using Pkg; Pkg.activate(projectdir)
# include("Turnover.jl")

using Test

using .Turnover
@testset "Applying turnover methods" begin

    using DataFrames, Statistics

    b, d, μ = params = (b=1., d=0.8, μ=0.1)
    @testset "simulation" begin

        out = neutral_growth(100 ; b = b, d = d, μ = μ, return_obs=true, showprogress=false)
        neutral_growth!(out[:tumor], out[:obs], 120 ; params..., showprogress=false)
        neutral_growth!(deepcopy(out[:tumor]), 130 ; params..., showprogress=false)

        @test !(:obs in keys( neutral_growth(10 ; params..., return_obs=false, showprogress=false) ) )

        @test begin
            out_clones = neutral_growth_clones(100; params..., showprogress=false )
            t = DataFrame(out_clones.tumor)
            all( round.(mfreqs(t).*sum(t.n), digits=3) .== out_clones.clonesizes )
        end
    end

    out = neutral_growth(1000 ; params..., return_obs=true, showprogress=false)
    htumor = DataFrame(out[:tumor])

    @testset "estranged" begin

        estr_res = estranged(htumor) |> df -> (count(df.isestranged), nrow(df))
        estr_treeless_res = filter(h-> !iszero(h.n), htumor) |> estranged_treeless |> df -> (sum(df.isestranged), sum(df.isgreen))
        @test estr_res == round.(Int, estr_treeless_res)

        q = (d + b*(μ/2)^2) / (b*(1-(μ/2))^2)
        estranged_expected(htumor, q, out[:obs][:,1])

    end

    @testset "orphaned" begin
        orph_res = orphaned_green(htumor) |> df -> ( count( df.isorphaned ) , nrow(df) )
        orph_red_res = orphaned_red(htumor) |> df -> ( count( unique(df).isorphaned ) , length(unique(df.mutation)) )
        @test orph_res == orph_red_res

        orph_res = orphaned_green(htumor) |> df -> ( count( df.isorphaned ) , nrow(df) )
        orph_treeless_res = filter(h-> !iszero(h.n), htumor) |> orphaned_green_treeless |> df -> (sum(df.isorphaned), sum(df.isgreen))
        @test orph_res == round.(Int, orph_treeless_res)

        orph_res = orphaned_red(htumor) |> df -> ( count( df.isorphaned ) , nrow(df) )
        orph_treeless_res = filter(h-> !iszero(h.n), htumor) |> orphaned_red_treeless |> df -> (sum(df.isorphaned), sum(df.isgreen))
        @test orph_res == round.(Int, orph_treeless_res)

        orphaned_green_expected(htumor, d/b, out[:obs][:,2])

        htumor, clonesizes, nbirth, N, t = neutral_growth_clones(1000; params..., showprogress=false)
        @test isapprox( mean( (d/b) .^ vcat(nbirth...) ),
                        mean( orphaned_red_expected(d/b, nbirth).isorphaned ) )

    end
end

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
struct NotImplementedError <: Exception
    msg::String
end

function Base.showerror(io::IO, err::NotImplementedError)
    print(io, "NotImplementedError: ")
    print(io, err.msg)
end


text:
using DataFrames
using Gadfly
using BoteSalvatICX
using FFAST
using NeXLUncertainties
using NeXLCore
using NeXLMatrixCorrection
using NeXLSpectrum

include(joinpath(pkgdir(NeXLSpectrum), "test", "runtests.jl"))

specs = [ loadspectrum(joinpath(pkgdir(NeXLSpectrum), "test", "K412 spectra", "III-E K412[$i][4].msa")) for i in 0:4 ]
plot(specs..., autoklms=true, klms = [ n"C", ], edges = [ n"Fe K", n"Ca K" ], coincidences = [ n"O K-L3", ]) |> SVG(tempname(), 10inch, 3inch)

#=
using NeXLSpectrum, Gadfly, DataFrames # Load the libraries
using PackageCompiler
PackageCompiler.create_sysimage(
    [ "Gadfly", "DataFrames", "BoteSalvatICX", "FFAST", "NeXLUncertainties", "NeXLCore", "NeXLMatrixCorrection", "NeXLSpectrum" ];
    sysimage_path=joinpath(homedir(), ".julia", "NeXLSysimage.dll"),
    precompile_execution_file=joinpath(pkgdir(NeXLSpectrum),"scripts", "precompile.jl"))
=#

## To use the sysimage, start Julia with the `-J` option
# julia -JC:\Users\username\.julia\NeXLSysimage.dll"

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
function pmatmul(A::SparseMatrixCSC, v::Vector{Float64}, slc::Vector{UnitRange{Int64}})

    x = @distributed (+) for s in slc
        A[:,s] * v[s]
    end

    return x
end

function pdot(u::Vector{Float64}, v::Vector{Float64}, slc::Vector{UnitRange{Int64}})

    x = @distributed (+) for s in slc
        dot(u[s], v[s])
    end

    return x
end

function pnorm(u::Vector{Float64}, slc::Vector{UnitRange{Int64}})

    return sqrt(pdot(u, u, slc))
end

function pgmres(A::SparseMatrixCSC, x::Vector{Float64}, r::Vector{Float64}, m::Int64, slc::Vector{UnitRange{Int64}})

    n = length(x0)
    Q = zeros(Float64, (n,m+1))
    H = zeros(Float64, (m+2,m+1))
    sn = zeros(Float64, m)
    cn = zeros(Float64, m)
    e1 = zeros(Float64, m+1)
    e1[1] = 1.0

    rnorm = pnorm(r, slc)
    beta = rnorm* e1
    Q[:,1] = r / rnorm

    # begin iteration
    nk = 0
    for k=1:m
        nk += 1

        # do arnoldi iteration
        q = pmatmul(A, Q[:,k], slc)
        for i=1:k
            H[i,k] = pdot(q, Q[:,i], slc)
            q = q - H[i,k] * Q[:,i]
        end
        qnorm = pnorm(q, slc)
        H[k+1,k] = qnorm
        Q[:,k+1] = q / qnorm

        # apply givens rotation
        for i=1:k-1
            temp = cn[i] * H[i,k] + sn[i] * H[i+1,k]
            H[i+1,k] = -sn[i] * H[i,k] + cn[i] * H[i+1,k]
            H[i,k] = temp
        end

        # update sn and cn
        t = sqrt(H[k,k]^2 + H[k+1,k]^2)
        cn[k] = H[k,k] / t
        sn[k] = H[k+1,k] / t

        # eliminate H[k+1, k]
        H[k,k] = cn[k] * H[k,k] + sn[k] * H[k+1,k]
        H[k+1,k] = 0.0

        # update residual
        beta[k+1] = -sn[k]*beta[k]
        beta[k] = cn[k]*beta[k]
    end

    # compute correction to x
    y = H[1:nk,1:nk] \ beta[1:nk]
    return x + Q[:,1:nk] * y
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test,Random
using LinearAlgebra

@testset "Jacobi–Davidson" begin

@test_throws ErrorException jd()

@bench @testset "Betcke-Voss" begin
@info "Test: Betcke-Voss"

Random.seed!(0);
@info "Testing a PEP"
nep = nep_gallery("pep0",60)
TOL = 1e-11;
λ,u = jd_betcke(nep, tol=TOL, maxit=55, neigs = 2, logger=displaylevel, v=ones(size(nep,1)), λ=0 ,errmeasure=ResidualErrmeasure(nep))
@info " Smallest eigenvalue found: $λ"
Dc,Vc = polyeig(nep,DefaultEigSolver)
c = sortperm(abs.(Dc))
verify_lambdas(2, nep, λ, u, TOL)

Random.seed!(0);
@info "Testing SG as inner solver"
nep = nep_gallery("real_quadratic")
nep = SPMF_NEP(get_Av(nep), get_fv(nep))
TOL = 1e-10;
λ,u=jd_betcke(Float64, nep, tol=TOL, maxit=4, logger = displaylevel, projtype = :Galerkin, inner_solver_method = SGIterInnerSolver(), v=ones(size(nep,1)), λ=0)
verify_lambdas(1, nep, λ, u, TOL)

Random.seed!(0);
@info "Testing IAR Cheb as projected solver"
nep = nep_gallery("dep0_sparse",40)
TOL = 1e-10;
λ,u = jd_betcke(ComplexF64, nep, tol=TOL, maxit=30, logger = displaylevel, inner_solver_method = IARChebInnerSolver(), v=ones(size(nep,1)), λ=0)
verify_lambdas(1, nep, λ, u, TOL)
λ0 = λ[1] # store these for the next test
v0 = vec(u)

Random.seed!(0);
@info "Testing convergence before starting"
λ,u=jd_betcke(nep, tol=TOL, maxit=25, neigs=1, logger=displaylevel, λ=λ0, v=v0)
verify_lambdas(1, nep, λ, u, TOL)

@info "Testing errors thrown"
nep = nep_gallery("pep0",10)
# Throw error if iterating more than the size of the NEP
@test_throws ErrorException λ,u=jd_betcke(nep, tol=TOL, maxit=60, logger = displaylevel, v=ones(size(nep,1)))
# SG requires Galerkin projection type to keep Hermitian
@test_throws ErrorException λ,u=jd_betcke(Float64, nep, tol=TOL, maxit=4, projtype = :PetrovGalerkin, inner_solver_method = SGIterInnerSolver(), v=ones(size(nep,1)))
# An undefined projection type
@test_throws ErrorException λ,u=jd_betcke(nep, tol=TOL, maxit=4, projtype = :MYNOTDEFINED, v=ones(size(nep,1)))
# Too many required eigenvalues, will not converge and hence throw an exception
@test_throws NEPCore.NoConvergenceException λ,u=jd_betcke(nep, tol=TOL, maxit=4, λ=10.0, neigs=1000, v=ones(size(nep,1)))

end

@bench @testset "Effenberger" begin
@info "Test: Effenberger"

Random.seed!(0);
TOL = 1e-10
nep = nep_gallery("pep0",250)
λ, u = jd_effenberger(nep, neigs=5, logger=displaylevel, tol=TOL, maxit=80, λ=0.82+0.9im, v=ones(ComplexF64,size(nep,1)))
verify_lambdas(5, nep, λ, u, TOL)

Random.seed!(0);
TOL = 1e-10
nep = nep_gallery("dep0",60)
λ, u = jd_effenberger(nep, neigs=3, logger=displaylevel, tol=TOL, maxit=55, λ=0.6+0im, v=ones(ComplexF64,size(nep,1)))#, inner_solver_method = IARChebInnerSolver())
verify_lambdas(3, nep, λ, u, TOL)

@info "Testing convergence before starting"
λ,u=jd_effenberger(nep, neigs=1, logger=displaylevel, tol=TOL, maxit=55, λ=λ[1], v=vec(u[:,1]))
verify_lambdas(1, nep, λ, u, TOL)

@info "Testing errors thrown"
nep = nep_gallery("pep0",50)
# Throw error if iterating more than the size of the NEP
@test_throws ErrorException λ, u = jd_effenberger(nep, tol=TOL, maxit=(size(nep,1)+1), v=ones(size(nep,1)))
# SG not possible with Effenberger
@test_throws ErrorException λ, u = jd_effenberger(nep, tol=TOL, maxit=40, inner_solver_method = SGIterInnerSolver(), v=ones(size(nep,1)))
# Too many required eigenvalues, will not converge and hence throw an exception
@test_throws NEPCore.NoConvergenceException λ, u = jd_effenberger(nep, neigs=1000, tol=TOL, maxit=4, v=ones(size(nep,1)), λ=10.0)

end

end


text:
using IterTools: product
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(k::Int)::Int
    digits = collect(join(1:k))
    possible_numbers = Dict(4 => [[1, 1]], 5 => [[1, 2]], 6 => [[1, 2]], 7 => [[1, 3], [2, 3]],
                            8 => [[1, 3], [2, 3]], 9 => [[1, 4], [2, 3]])
    products = []
    for (i, j) ∈ possible_numbers[k]
        for n ∈ product(fill(digits, i)...)
            if length(∪(n)) == i
                for m ∈ product(fill(setdiff(digits, n), j)...)
                    if length(∪(m)) == j
                        candidate = parse(Int, join(n)) * parse(Int, join(m))
                        if sort(vcat(collect(n), collect(m), collect(string(candidate)))) == digits
                            append!(products, candidate)
                        end
                    end
                end
            end
        end
    end
    return sum(∪(products))
end

compute(9)

@benchmark compute(9)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
### A Pluto.jl notebook ###
# v0.16.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

# ╔═╡ e9576706-600e-11eb-1e10-e3bac02a254e
# edit the code below to set your name and UGent username

student = (name = "Hanne Janssen", email = "Hanne.Janssen@UGent.be");

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

# you might need to wait until all other cells in this notebook have completed running.
# scroll down the page to see what's up

# ╔═╡ fa42d7da-600e-11eb-13a4-7dfe5ebbafd0
begin
	using Pkg; Pkg.activate("../..")
	using DSJulia;
	using PlutoUI;
	tracker = ProgressTracker(student.name, student.email);
	md"""
	Submission by: **_$(student.name)_**
	"""
end

# ╔═╡ 4ec271b0-4e73-11eb-2660-6b8bd637d7ee
md"""
# Abstract and primitive types

All Julia objects, both those already defined as well as those you might invent yourself, have a type. The type system is the secret *sauce*, allowing Julia to be fast because code can be specialised for a particular combination of types. It is also supremely useful in conjunction with *multiple dispatch*, in which functions work differently depending on which types you feed into them. This notebook will show the basics of the type system. We will only consider abstract and primitive types here, leaving composite types for the next notebook. Don't worry if you do not yet understand what this means, you will discover it soon enough.
"""

# ╔═╡ a1f2d06e-4e73-11eb-3afd-1353def71700
md"""
## Checking the type

The type of objects can be assessed using the function `typeof`. For collections, `eltype` gives the types of individual elements. Try the following examples. Note that types are always capitalised!
"""

# ╔═╡ c0bfdf9e-4e73-11eb-3962-0b3c5d5424d7
a = 42; s = "mice"; n = 0.9; A = [1 2; 3 4];

# ╔═╡ bd994d64-600e-11eb-1ab3-ed6317b7c211
begin
	QuestionBlock(
		title=md"**Task:**",
		description = md"""
		Check the type of the `a`, `s`, `n` and `A`

		```julia
		typeof(a)

		typeof(s)

		typeof(n)

		typeof(A)

		```
		"""
	)
end

# ╔═╡ b844d568-4e73-11eb-3de9-4158b0bdca12

# ╔═╡ c662744a-4e73-11eb-1bfc-6daaf7282285

# ╔═╡ cae803e2-4e73-11eb-13e0-23abccf86bac

# ╔═╡ cc606026-4e73-11eb-3576-5d301a771a5a

# ╔═╡ d3803112-4e73-11eb-2018-f72ffb7f6ec6
md"These are all *concrete types*. Julia types are part of a hierarchical type system, forming a single, fully connected type graph. The concrete types are the leaves of this tree, whereas the inner nodes are *abstract types*. As hinted by the name, these are abstract and cannot be instantiated. They, however, help with conceptually ordering the type system."

# ╔═╡ e8fe2832-4e78-11eb-3fd6-a3f1f0c0892a
md"What is the type of `pi`?"

# ╔═╡ f533012c-4e78-11eb-0f45-3b47f088c9c6
typeofpi = missing

# ╔═╡ f69d89ba-4e73-11eb-3ab9-9179ea7e3217
md"Concrete types (should) have a well-defined memory layout, for example `Float64` is  encoded using 64 bits while `Float32` is encoded using 32 bits and hence some computations can be executed quicker but less precise by the former. Abstract types on the other hand mainly encode a semantic meaning, any `Real` should behave as a real number (e.g., addition, division are defined)."

# ╔═╡ de3de7fc-4e73-11eb-2ff6-1560481f7ee5
md"We can find the supertype (ancestor) of a concrete or abstract type using the function `supertype`."

# ╔═╡ 910c8e0e-600f-11eb-2558-57040714937b
begin
	QuestionBlock(
		title=md"**Task:**",
		description = md"""
		Check the supertype of the following concrete and abstract types,

		```julia
		supertype(Int8)

		supertype(Float64)

		supertype(AbstractFloat)

		supertype(Real)

		supertype(Number)

		supertype(Any)
		```
		"""
	)
end

# ╔═╡ e1c8cf4a-4e73-11eb-27be-d702064a0182

# ╔═╡ e56a46c6-4e73-11eb-1748-1b6fe5ab0376

# ╔═╡ e89adcaa-4e73-11eb-1ed8-e9c89ca633f6

# ╔═╡ ec2ab2be-4e73-11eb-1a22-010439761432

# ╔═╡ efa205b4-4e73-11eb-1647-e9dcab5f7b7a

# ╔═╡ f3b5a778-4e73-11eb-1d3c-11ae19713eca

# ╔═╡ a0cecb24-4e74-11eb-3634-cd8dd628e9ec
md"See how all the numbers are hierarchically represented? Note that any type is always a subtype of `Any`. We can check if an object is (sub)type using the function `isa` or use the `<:` operator."

# ╔═╡ 200d8d36-6012-11eb-0171-8bc2511e18f1
begin
	QuestionBlock(
		title=md"**Task:**",
		description = md"""
		Predict the outcome of the following statements,

		```julia
		Float64 <: AbstractFloat

		Float16 <: AbstractFloat

		AbstractFloat <: Number

		Int <: Number

		Int <: AbstractFloat

		Integer isa Int

		```
		""",
		hints = [
			hint(md"If you are confused by the last statement, read the next section.")
		]
	)
end

# ╔═╡ b31fe65a-4e74-11eb-0414-35f2be687c7f

# ╔═╡ c2ac0c48-4e74-11eb-10b0-91ad620fefcd

# ╔═╡ c5208cce-4e74-11eb-0615-135b510a9e8d

# ╔═╡ c817296a-4e74-11eb-0994-972871114f02

# ╔═╡ cb066442-4e74-11eb-35e7-ed38d4bd8bbf

# ╔═╡ ce3d5380-4e74-11eb-3d9d-5f34cbbae118

# ╔═╡ 66343826-6012-11eb-109c-17c7a582cbc8

# ╔═╡ 34ca6158-4e75-11eb-3d51-330952c9b3dd
md"We can check the entire subtree of a type using the function `subtypetree`"

# ╔═╡ 972154f6-6536-11eb-08a8-e3fc98623177
function subtypetree(roottype, level=1, last=false)
	level == 1 && println(split(string(roottype),".")[end])
	for (index, s) in enumerate(subtypes(roottype))

		if index < length(subtypes(roottype))
			println(join(fill("│  ", level-1-1*last))*
				join(fill("   ", 1*last)) * "├──" * split(string(s),".")[end])

			subtypetree(s, level + 1, false)
		else
			println(join(fill("│  ", level-1-1*last))*
				join(fill("   ", 1*last)) * "└──" * split(string(s),".")[end])
			subtypetree(s, level + 1, true)
		end

	end
end

# ╔═╡ e50c3f36-653c-11eb-0eb4-cf44af78e44c
@terminal subtypetree(Real)

# ╔═╡ 3d1db2de-6549-11eb-2649-2d96659813f6
md"It does not need to be complicated, though."

# ╔═╡ 46eaafee-6549-11eb-2e36-c9e566d3f3ba
function subtypetreesimple(roottype, level=1)
	level == 1 && println(roottype)
	for s in subtypes(roottype)
			println(" "^((level-1)*4) * string(s))
			subtypetreesimple(s, level + 1)
	end
end

# ╔═╡ 6cbd4874-6549-11eb-296e-af4a43d53223
@terminal subtypetreesimple(Real)

# ╔═╡ 4a01487c-4e78-11eb-1302-d9c6ec4ed6ab
md"""

## Converting types

Changing a variable `x` into a different type `T` can be done using the function `convert(T, x)`.
"""

# ╔═╡ 66292626-4e78-11eb-331b-0563b2110605
convert(Float64, 42)

# ╔═╡ 704b2ebe-4e78-11eb-1583-d10e0aeb2b8d
Float64(42)  # for most types this works as well

# ╔═╡ c3c40df2-4e78-11eb-3d4a-5fdfdf173da3
Int8(42)

# ╔═╡ cd32e96c-4e78-11eb-0b48-5767421c7875
Float32(π)

# ╔═╡ 40a761c2-5b24-11eb-09a8-a5cd0bc4ab95
md"We have seen that you can add any type of float with any type of integer (ditto for vectors and arrays with different types but how does this work? Julia uses *promotion* to cast two inputs in the more general type. For example, compare adding a `Float64` with an `Int`."

# ╔═╡ 99228eee-5b24-11eb-385e-7507ca20ae0e
promote(7.9, 79)

# ╔═╡ 03133d24-6621-11eb-03d3-bb8970ab19f7
promote_type(Float64, Int)  # find common type

# ╔═╡ ba39991a-5b24-11eb-260b-439bcde4c153
md"You see that `Float64` is the more general type, so both inputs are cast as floats and further processed by the function. That is why the their sum is a float: `7.9 + 79 = 86.9`.

Of course, this also works with more complex composite types, such as matrices:
"

# ╔═╡ 35e53434-5b25-11eb-10b7-e993e9477c8c
[1 2; 3 4] + [0.0 1.0; 2.0 3.0]

# ╔═╡ 92dc7a4e-5b25-11eb-1518-8182216f24ec
[0, 1.0, 2, 3]  # cast into a vector of floats

# ╔═╡ aa26b46c-4e78-11eb-24d8-7fdce7c94fff
md"""
When designing new types, one also often has to implement custom `convert` methods since it is not always clear how to convert different types.

When reading files, it is often useful to interpret parts of the strings as different datatypes such as numbers. This can be done using the function `parse` that works similarly to `convert`.
"""

# ╔═╡ 815b0436-4e78-11eb-13d4-0dc6531e34f2
parse(Int, "42")

# ╔═╡ 570c85bc-4e79-11eb-0249-891cf205d623
parse(Float64, "0.999")

# ╔═╡ 6756d6ac-4e79-11eb-21ab-4776195c9d3b
bunchofnumbers = "1.728002758512114, 0.45540258865644284, 1.4067738604851092, 1.6549474922755167, -0.5281073122489854, 2.219250973007533, 0.8195027302254512, 1.8833469318073521, 0.7429034224663096, -0.8127686064960085, -0.14337850083375886, -1.477193046160141, 0.024525761924498457, 0.16097115910472956, -0.39278880092280993, 1.3988081686729814, -1.3316370350161346, 0.2791510437718087, 1.9834455917052212, -0.8616791621501649
"

# ╔═╡ e6f31ad8-4e79-11eb-11f4-2936cb039f8d
sumofbunchofnumbers = missing

# ╔═╡ 4495802e-6013-11eb-22fa-03f635d30a7b
begin
	q_sparse = Question(
		validators=@safe[
			sumofbunchofnumbers == Solutions.bunchofnumbers_parser(bunchofnumbers)
		]
	)

	qb_sparse = QuestionBlock(
		title=md"**Question: string parsing**",
		description = md"""
		Below are a bunch of numbers in a text string. Can you compute their sum `(Float64)`?

		""",
		questions = [q_sparse],
		hints = [
			hint(md" `rstrip` is a handy function."),
			hint(md" `split` is a handy function."),
			hint(md" Maybe we can use a pipe `|>`?"),
		]
	)
	validate(qb_sparse, tracker)
end

# ╔═╡ 03766a5c-4e75-11eb-12ad-cb2e9468e0d2
md"""
## Methods and dispatch

When a function is run for the first time with a particular combination of input types, it gets compiled by the LLVM compiler. Such a specific function is referred to as a `method`. Every time a function is run with a new combination of types of arguments, a suitable method is compiled. This is noticeable when measuring the running time.
"""

# ╔═╡ 63166056-6014-11eb-09c4-e5a44d37095f
begin
	QuestionBlock(
		title=md"**Task:**",
		description = md"""
		Run the following examples in the terminal using `@time`.

		```julia
		@time mynewfun(1)  # returns an integer

		@time mynewfun(1.0)  # returns a Float64

		@time mynewfun(A)

		```
		"""
	)
end

# ╔═╡ 2dff8c88-4e75-11eb-050b-7152e82ac10d
mynewfun(x) = x^2 .+ x

# ╔═╡ 7c2b6dc0-4e76-11eb-1d78-553df82d9100

# ╔═╡ d2a4a32c-5b02-11eb-3839-8108c4965931

# ╔═╡ 32d64b6e-4e75-11eb-0a2a-27214f217f70

# ╔═╡ 861ba4c6-4e76-11eb-3d2b-bfabbd143df2
md"The known methods can be found using the function `methods`. For example, look how many methods are defined for sum:"

# ╔═╡ 8d5f7d8e-4e76-11eb-28ba-bdec03a3e150
methods(sum)

# ╔═╡ b4fac50a-6015-11eb-16d4-594858d78b3e

# ╔═╡ 9ed7cb5a-6014-11eb-0ae8-eba8d77867a2
begin
	QuestionBlock(
		title=md"**Did I break something?**",
		description = md"""

		check how many methods there are associated with the humble multiplication operator `*`. Just for your info, you might want to print this one to the terminal.
		"""
	)
end

# ╔═╡ b18d0532-4e76-11eb-2e8a-2bee580533cc

# ╔═╡ b3d15950-6015-11eb-1909-c127822a4a83

# ╔═╡ cd2eaafc-4e76-11eb-245a-e9898d3d57a4
md"The arguments a function can take can be restricted using the `::`-operator. Here, if we limit a function as `f(x::T)`, this means that `x` can be any type `<: T`. "

# ╔═╡ 002fdec6-6015-11eb-0e89-c7d020826cf9
begin
	QuestionBlock(
		title=md"**Question: ::**",
		description = md"""

		Can you explain the reasoning behind the following code? How does it process numbers? What does it do with strings?"

		```julia
		methods(twice)
		twice(10) # Int
		twice(10.0)  # Float64, also a Number but not an Int
		twice("A griffin! ")  # strings mean something else

		```

		"""
	)
end

# ╔═╡ db1bb4c8-4e76-11eb-2756-3f6ce778acc0
begin
	twice(x::Number) = 2x;
	twice(x::AbstractString) = x * x;
end

# ╔═╡ ff755bf8-4e76-11eb-205f-d52529ae50ed

# ╔═╡ 03932e5e-4e77-11eb-3769-635cc33c3c4d

# ╔═╡ 0b4f99ea-4e77-11eb-29fc-632788d179a3

# ╔═╡ 2a0b220a-4e77-11eb-1da7-2978422c11f4

# ╔═╡ bf91e40a-4e77-11eb-14f1-754b1ce5130e
md"> Julia will always select the method with the most specific type signature.
So, if we would define a function `twice(x::Float64)`, it would be chosen to process `Float64` inputs, even though these are also `Number`s."

# ╔═╡ e1a88a70-4e76-11eb-2486-e1d2f4211792
begin
	f(x, y) = "No life forms present";
	f(x::T, y::T) where {T} = x * y;  # short for {T <: Any}
	f(x::Integer, y::Real) = 2x + y;
	f(x::Int, y::Int) = 2x + 2y;
	f(x::Integer, y::Float64) = x + 2y;
	f(x::Float64, y::Real) = x - y;
	f(x::Float64, y::Float64) = 2x - y;
end

# ╔═╡ 622b8382-6015-11eb-17fb-3352c73a0d10
begin
	QuestionBlock(
		title=md"**Question: ::²**",
		description = md"""

		Predict the outcome of the following statements.

		```julia
		f(1, 2.0)

		f(1.0, 2)

		f(Int8(1), Int8(2))

		f(1.0, 2.0)

		f("one", 2)

		f("one", "two")

		f(1, Float32(2.0))

		f(1, 2)

		f([1 1; 1 1], [2.0 2.0; 2.0 2.0])

		f([1 1; 1 1], [2 2; 2 2])

		```

		"""
	)
end

# ╔═╡ 76fe9fc4-4e77-11eb-3bc7-2dfbdff8dfc8

# ╔═╡ 7aa14c94-4e77-11eb-25c7-fb0103267b06

# ╔═╡ 7f3a5336-4e77-11eb-2ad6-3d889dc75ac0

# ╔═╡ 822c01d4-4e77-11eb-1409-fbaf83c950b6

# ╔═╡ 85f186a6-4e77-11eb-19ca-5db29615ba97

# ╔═╡ 891c820e-4e77-11eb-1ebf-b3065e0d4211

# ╔═╡ 8d0d39c4-4e77-11eb-034d-07dc33ab6e9a

# ╔═╡ 901efaee-4e77-11eb-02d9-b5fe1f0931d5

# ╔═╡ 938d8b1e-4e77-11eb-03d3-9b88c7cab3c1

# ╔═╡ 96f6fef2-4e77-11eb-2ec4-399472d86a60

# ╔═╡ 812cfe48-4e7a-11eb-32e6-c918bbe3e602
md"""

## Extending the type system

Being aware of the type system is a first step but the fun is in extending it and creating your own types.

### Abstract types

Abstract types are defined using the following simple syntax:

```
abstract type «name» end
abstract type «name» <: «supertype» end
```

### Primitive types

*Primitive types* exist of simple bits. Examples are `Float64` and `Int16`. You can declare your own types, though this is likely not something many often do in practice. For example, there is a specific binary encoding possible for nucleotides to make bioinformatics [computations more efficient](https://medium.com/analytics-vidhya/bioinformatics-2-bit-encoding-for-dna-sequences-9b93636e90e2).

### Composite types
*Composite types* (records, structs, or objects) are more exciting. They are often containers for several objects set to behave in a certain way. We will study them in depth in the next notebook.
"""

# ╔═╡ 15bb7eb6-6016-11eb-01d2-0be88279db7e

# ╔═╡ 178a5f96-6016-11eb-0314-010203c5fedf
md"It takes some time to fully grasp the potential of julia's type system, especially if you have experience in object-oriented programming. In the following examples we will try to show you the uniqueness and power of this paradigm."

# ╔═╡ 4df6c0c0-4e7c-11eb-1d43-0d9bbf4896a7
md"""
## Case study 1: Mohs scale

![](https://i.imgur.com/WtyJ2Uq.png)
[source ](http://www.911metallurgist.com/blog/wp-content/uploads/2015/08/Mohs-Hardness-Test-Kit-and-Scale.jpg)
"""

# ╔═╡ d668a692-6017-11eb-2f28-e9b00762b92d
md"The Mohs scale of mineral hardness is a qualitative ordinal scale characterizing scratch resistance of various minerals through the ability of harder material to scratch softer material. [wikipedia](https://en.wikipedia.org/wiki/Mohs_scale_of_mineral_hardness)"

# ╔═╡ fb1d0b9a-6017-11eb-0ba1-67e204f150fd
md"First we define the abstract type tree structure,"

# ╔═╡ b487a776-4e7e-11eb-291b-e900e6e1a2f6
begin
	abstract type Mohs end

	abstract type Diamond <: Mohs end
	abstract type Corundum <: Mohs end
	abstract type Topaz <: Mohs end
	abstract type Quartz <: Mohs end
	abstract type Orthoclase <: Mohs end
	abstract type Apatite <: Mohs end
	abstract type Fluorite <: Mohs end
	abstract type Calcite <: Mohs end
	abstract type Gypsum <: Mohs end
	abstract type Talc <: Mohs end
end

# ╔═╡ 0884f752-6018-11eb-2eb3-d1cd317dceb3
md"You can see, it is a pretty flat hierarchy"

# ╔═╡ 24588124-6018-11eb-24d1-f9c7759f4c8f
@terminal subtypetree(Mohs)  # check the terminal!

# ╔═╡ 443ecf72-6018-11eb-1a7a-e75e9596e4bd
md"Next, let us define a function `mohs_scale` that dispatches on the different abstract types (minerals) and returns a hardness value"

# ╔═╡ 1aea83a8-4e7f-11eb-2d06-c3e550c4e1b9
begin
	mohs_scale(::Type{Diamond}) = 10
	mohs_scale(::Type{Corundum}) = 9
	mohs_scale(::Type{Topaz}) = 8
	mohs_scale(::Type{Quartz}) = 7
	mohs_scale(::Type{Orthoclase}) = 6
	mohs_scale(::Type{Apatite}) = 5
	mohs_scale(::Type{Fluorite}) = 4
	mohs_scale(::Type{Calcite}) = 3
	mohs_scale(::Type{Gypsum}) = 2
	mohs_scale(::Type{Talc}) = 1
end

# ╔═╡ 70390e06-6018-11eb-3bd8-97ff84985261
begin
	💎 = Diamond
	🔶 = Topaz
end;

# ╔═╡ 9b4cad6e-6018-11eb-0e41-2fea9c9219ed
supertype(💎)

# ╔═╡ ab1c59c4-6018-11eb-1147-8d441c25aa0b
mohs_scale(💎)

# ╔═╡ e36b2aa8-6018-11eb-0582-997ba8b7eae7
mohs_scale(🔶)

# ╔═╡ b5339e22-6018-11eb-38b6-352602a85cdd
md"Will 💎 scratch 🔶?"

# ╔═╡ fc0f2e92-6018-11eb-0b45-49fdec301609
mohs_scale(💎) > mohs_scale(🔶)

# ╔═╡ 02e57e74-6019-11eb-3e35-1d2b4d61b283
md"To make this more user-friendly, one can add a method to the `<`-operator or the `isless`-function to work directly on `mohs_scale`."

# ╔═╡ 6e429164-4e7f-11eb-1829-0582f1417815
Base.isless(m1::Type{<:Mohs}, m2::Type{<:Mohs}) = mohs_scale(m1) < mohs_scale(m2)

# ╔═╡ ad6e9d8e-4e7f-11eb-1e33-efcee699f2a0
isless(Diamond, Corundum)

# ╔═╡ 9d1c94ba-4e7f-11eb-060e-d1bc9683af92
Calcite < Fluorite < Corundum

# ╔═╡ 8dc6a6ce-6019-11eb-1cd5-d1c1c118f4ee
md"`>` is defined using the `isless` function, so this just works,"

# ╔═╡ 73f66dc6-6019-11eb-06b8-7db54f428e12
💎 > 🔶

# ╔═╡ ee370b50-4e7f-11eb-1ce0-d1bdb3e41ae2
rocks = [Gypsum, Orthoclase, Quartz, Corundum, Fluorite, Gypsum, Talc]

# ╔═╡ 1117c600-4e80-11eb-3231-4383d700f760
md"Even cooler, `sort` and `sort!` are also internally defined using the `isless` function, so these just work out of the box:"

# ╔═╡ 07a19574-4e80-11eb-38fa-8d3463dfd700
sort(rocks)

# ╔═╡ af5bdeb4-6019-11eb-2265-a9b81e01d9c0

# ╔═╡ 84a0f49c-4e7c-11eb-14f2-452e57f2e414
md"""
## Case study 2: rock-paper-scissors

![Fire, water, grass!](https://upload.wikimedia.org/wikipedia/commons/thumb/6/67/Rock-paper-scissors.svg/1200px-Rock-paper-scissors.svg.png)

We can easily implement the rock-paper-scissors rules using types.
"""

# ╔═╡ 99c4f3c8-4e7c-11eb-3d4a-33ba8d495eb2
begin
	abstract type Hand end

	abstract type Rock <: Hand end
	abstract type Paper <: Hand end
	abstract type Scissors <: Hand end
end

# ╔═╡ cba6d4cc-5b03-11eb-265d-3f08117b0e8d
md"Now we implement a function to play one hand against an opponent's hand."

# ╔═╡ d913cd92-4e7c-11eb-11e1-3d7539af7fed
begin
	play(h1::Type{Paper}, h2::Type{Rock}) = 1
	play(h1::Type{Rock}, h2::Type{Scissors}) = 1
	play(h1::Type{Scissors}, h2::Type{Paper}) = 1

	# this captures both when same inputs are given and
	# when the first person looses
	play(h1::Type{<:Hand}, h2::Type{<:Hand}) = h1 == h2 ? 0 : -1
end

# ╔═╡ 4f107d88-4e7d-11eb-3e49-f54ecf5163da
play(Rock, Rock)

# ╔═╡ 8331c8b0-4e7d-11eb-0690-8bbae3ed086a
play(Rock, Scissors)

# ╔═╡ 88a95ec0-4e7d-11eb-0a33-77ef82874f45
play(Scissors, Rock)

# ╔═╡ 925e2f40-4e7d-11eb-0bd2-f91913c5a23e
play(Scissors, Paper)

# ╔═╡ 0b2b40ae-601a-11eb-1d83-53ff171bb799
begin
	QuestionBlock(
		title=md"**Optional question: rock, paper, scissors, lizard, Spock**",
		description = md"""

		Can you extend the previous code so that it works with lizard an Spock?

		![](https://i.imgur.com/jbmAYKI.png)

		adapted from: [source](https://external-content.duckduckgo.com/iu/?u=http%3A%2F%2Fwordpress.morningside.edu%2Fcdl001%2Ffiles%2F2010%2F09%2FRockPaperScissorsLizardSpock.jpg&f=1&nofb=1)
		"""
	)
end

# ╔═╡ 269b934c-601b-11eb-00ad-5fec0e2c37e1

# ╔═╡ ce5d564a-f2f3-4b1e-aa70-85253b2ccf38
md"""## Answers:
If you would like to take a look at the answers, you can do so by checking the boxo the question you would like to see. The function will be shown just below the question you want to look at.

| Question | Show solution |
|-----|:---------:|
| Question string parsing | $(@bind answ_q1 CheckBox()) |
| Question Rock,Paper,Spock | $(@bind answ_q2 CheckBox()) |

"""

# ╔═╡ 9a33c306-653c-11eb-3373-01e90d12b246
if answ_q1 == true
	md"""
	```Julia
	sumofbunchofnumbers(bunchofnumbers) = parse.(Float64, split(rstrip(bunchofnumbers), ", ")) |> sum
	```
	"""
end

# ╔═╡ 3cf6dffd-d91e-465a-94bb-0ffe6b5152aa
if answ_q2 == true
	md"""
	```Julia
	abstract type Lizard <: Hand end
	abstract type Spock <: Hand end

	play(h1::Type{Paper}, h2::Type{Spock}) = 1
	play(h1::Type{Rock}, h2::Type{Lizard}) = 1
	play(h1::Type{Lizard}, h2::Type{Paper}) = 1
	play(h1::Type{Lizard}, h2::Type{Spock}) = 1
	play(h1::Type{Spock}, h2::Type{Scissors}) = 1
	play(h1::Type{Spock}, h2::Type{Rock}) = 1
	play(h1::Type{Scissors}, h2::Type{Lizard}) = 1
	```
	"""
end

# ╔═╡ Cell order:
# ╠═e9576706-600e-11eb-1e10-e3bac02a254e
# ╟─fa42d7da-600e-11eb-13a4-7dfe5ebbafd0
# ╟─4ec271b0-4e73-11eb-2660-6b8bd637d7ee
# ╟─a1f2d06e-4e73-11eb-3afd-1353def71700
# ╠═c0bfdf9e-4e73-11eb-3962-0b3c5d5424d7
# ╟─bd994d64-600e-11eb-1ab3-ed6317b7c211
# ╠═b844d568-4e73-11eb-3de9-4158b0bdca12
# ╠═c662744a-4e73-11eb-1bfc-6daaf7282285
# ╠═cae803e2-4e73-11eb-13e0-23abccf86bac
# ╠═cc606026-4e73-11eb-3576-5d301a771a5a
# ╟─d3803112-4e73-11eb-2018-f72ffb7f6ec6
# ╟─e8fe2832-4e78-11eb-3fd6-a3f1f0c0892a
# ╠═f533012c-4e78-11eb-0f45-3b47f088c9c6
# ╟─f69d89ba-4e73-11eb-3ab9-9179ea7e3217
# ╟─de3de7fc-4e73-11eb-2ff6-1560481f7ee5
# ╟─910c8e0e-600f-11eb-2558-57040714937b
# ╠═e1c8cf4a-4e73-11eb-27be-d702064a0182
# ╠═e56a46c6-4e73-11eb-1748-1b6fe5ab0376
# ╠═e89adcaa-4e73-11eb-1ed8-e9c89ca633f6
# ╠═ec2ab2be-4e73-11eb-1a22-010439761432
# ╠═efa205b4-4e73-11eb-1647-e9dcab5f7b7a
# ╠═f3b5a778-4e73-11eb-1d3c-11ae19713eca
# ╟─a0cecb24-4e74-11eb-3634-cd8dd628e9ec
# ╟─200d8d36-6012-11eb-0171-8bc2511e18f1
# ╠═b31fe65a-4e74-11eb-0414-35f2be687c7f
# ╠═c2ac0c48-4e74-11eb-10b0-91ad620fefcd
# ╠═c5208cce-4e74-11eb-0615-135b510a9e8d
# ╠═c817296a-4e74-11eb-0994-972871114f02
# ╠═cb066442-4e74-11eb-35e7-ed38d4bd8bbf
# ╠═ce3d5380-4e74-11eb-3d9d-5f34cbbae118
# ╟─66343826-6012-11eb-109c-17c7a582cbc8
# ╟─34ca6158-4e75-11eb-3d51-330952c9b3dd
# ╠═972154f6-6536-11eb-08a8-e3fc98623177
# ╠═e50c3f36-653c-11eb-0eb4-cf44af78e44c
# ╟─3d1db2de-6549-11eb-2649-2d96659813f6
# ╠═46eaafee-6549-11eb-2e36-c9e566d3f3ba
# ╠═6cbd4874-6549-11eb-296e-af4a43d53223
# ╟─4a01487c-4e78-11eb-1302-d9c6ec4ed6ab
# ╠═66292626-4e78-11eb-331b-0563b2110605
# ╠═704b2ebe-4e78-11eb-1583-d10e0aeb2b8d
# ╠═c3c40df2-4e78-11eb-3d4a-5fdfdf173da3
# ╠═cd32e96c-4e78-11eb-0b48-5767421c7875
# ╟─40a761c2-5b24-11eb-09a8-a5cd0bc4ab95
# ╠═99228eee-5b24-11eb-385e-7507ca20ae0e
# ╠═03133d24-6621-11eb-03d3-bb8970ab19f7
# ╟─ba39991a-5b24-11eb-260b-439bcde4c153
# ╠═35e53434-5b25-11eb-10b7-e993e9477c8c
# ╠═92dc7a4e-5b25-11eb-1518-8182216f24ec
# ╟─aa26b46c-4e78-11eb-24d8-7fdce7c94fff
# ╠═815b0436-4e78-11eb-13d4-0dc6531e34f2
# ╠═570c85bc-4e79-11eb-0249-891cf205d623
# ╟─4495802e-6013-11eb-22fa-03f635d30a7b
# ╠═6756d6ac-4e79-11eb-21ab-4776195c9d3b
# ╠═e6f31ad8-4e79-11eb-11f4-2936cb039f8d
# ╟─9a33c306-653c-11eb-3373-01e90d12b246
# ╟─03766a5c-4e75-11eb-12ad-cb2e9468e0d2
# ╟─63166056-6014-11eb-09c4-e5a44d37095f
# ╠═2dff8c88-4e75-11eb-050b-7152e82ac10d
# ╠═7c2b6dc0-4e76-11eb-1d78-553df82d9100
# ╠═d2a4a32c-5b02-11eb-3839-8108c4965931
# ╠═32d64b6e-4e75-11eb-0a2a-27214f217f70
# ╟─861ba4c6-4e76-11eb-3d2b-bfabbd143df2
# ╠═8d5f7d8e-4e76-11eb-28ba-bdec03a3e150
# ╟─b4fac50a-6015-11eb-16d4-594858d78b3e
# ╟─9ed7cb5a-6014-11eb-0ae8-eba8d77867a2
# ╠═b18d0532-4e76-11eb-2e8a-2bee580533cc
# ╟─b3d15950-6015-11eb-1909-c127822a4a83
# ╟─cd2eaafc-4e76-11eb-245a-e9898d3d57a4
# ╟─002fdec6-6015-11eb-0e89-c7d020826cf9
# ╠═db1bb4c8-4e76-11eb-2756-3f6ce778acc0
# ╠═ff755bf8-4e76-11eb-205f-d52529ae50ed
# ╠═03932e5e-4e77-11eb-3769-635cc33c3c4d
# ╠═0b4f99ea-4e77-11eb-29fc-632788d179a3
# ╠═2a0b220a-4e77-11eb-1da7-2978422c11f4
# ╟─bf91e40a-4e77-11eb-14f1-754b1ce5130e
# ╠═e1a88a70-4e76-11eb-2486-e1d2f4211792
# ╟─622b8382-6015-11eb-17fb-3352c73a0d10
# ╠═76fe9fc4-4e77-11eb-3bc7-2dfbdff8dfc8
# ╠═7aa14c94-4e77-11eb-25c7-fb0103267b06
# ╠═7f3a5336-4e77-11eb-2ad6-3d889dc75ac0
# ╠═822c01d4-4e77-11eb-1409-fbaf83c950b6
# ╠═85f186a6-4e77-11eb-19ca-5db29615ba97
# ╠═891c820e-4e77-11eb-1ebf-b3065e0d4211
# ╠═8d0d39c4-4e77-11eb-034d-07dc33ab6e9a
# ╠═901efaee-4e77-11eb-02d9-b5fe1f0931d5
# ╠═938d8b1e-4e77-11eb-03d3-9b88c7cab3c1
# ╠═96f6fef2-4e77-11eb-2ec4-399472d86a60
# ╟─812cfe48-4e7a-11eb-32e6-c918bbe3e602
# ╟─15bb7eb6-6016-11eb-01d2-0be88279db7e
# ╟─178a5f96-6016-11eb-0314-010203c5fedf
# ╟─4df6c0c0-4e7c-11eb-1d43-0d9bbf4896a7
# ╟─d668a692-6017-11eb-2f28-e9b00762b92d
# ╟─fb1d0b9a-6017-11eb-0ba1-67e204f150fd
# ╠═b487a776-4e7e-11eb-291b-e900e6e1a2f6
# ╟─0884f752-6018-11eb-2eb3-d1cd317dceb3
# ╠═24588124-6018-11eb-24d1-f9c7759f4c8f
# ╟─443ecf72-6018-11eb-1a7a-e75e9596e4bd
# ╠═1aea83a8-4e7f-11eb-2d06-c3e550c4e1b9
# ╠═70390e06-6018-11eb-3bd8-97ff84985261
# ╠═9b4cad6e-6018-11eb-0e41-2fea9c9219ed
# ╠═ab1c59c4-6018-11eb-1147-8d441c25aa0b
# ╠═e36b2aa8-6018-11eb-0582-997ba8b7eae7
# ╠═b5339e22-6018-11eb-38b6-352602a85cdd
# ╠═fc0f2e92-6018-11eb-0b45-49fdec301609
# ╟─02e57e74-6019-11eb-3e35-1d2b4d61b283
# ╠═6e429164-4e7f-11eb-1829-0582f1417815
# ╠═ad6e9d8e-4e7f-11eb-1e33-efcee699f2a0
# ╠═9d1c94ba-4e7f-11eb-060e-d1bc9683af92
# ╟─8dc6a6ce-6019-11eb-1cd5-d1c1c118f4ee
# ╠═73f66dc6-6019-11eb-06b8-7db54f428e12
# ╠═ee370b50-4e7f-11eb-1ce0-d1bdb3e41ae2
# ╟─1117c600-4e80-11eb-3231-4383d700f760
# ╠═07a19574-4e80-11eb-38fa-8d3463dfd700
# ╟─af5bdeb4-6019-11eb-2265-a9b81e01d9c0
# ╟─84a0f49c-4e7c-11eb-14f2-452e57f2e414
# ╠═99c4f3c8-4e7c-11eb-3d4a-33ba8d495eb2
# ╟─cba6d4cc-5b03-11eb-265d-3f08117b0e8d
# ╠═d913cd92-4e7c-11eb-11e1-3d7539af7fed
# ╠═4f107d88-4e7d-11eb-3e49-f54ecf5163da
# ╠═8331c8b0-4e7d-11eb-0690-8bbae3ed086a
# ╠═88a95ec0-4e7d-11eb-0a33-77ef82874f45
# ╠═925e2f40-4e7d-11eb-0bd2-f91913c5a23e
# ╟─0b2b40ae-601a-11eb-1d83-53ff171bb799
# ╠═269b934c-601b-11eb-00ad-5fec0e2c37e1
# ╟─3cf6dffd-d91e-465a-94bb-0ffe6b5152aa
# ╟─ce5d564a-f2f3-4b1e-aa70-85253b2ccf38


text:
using Base.Iterators: IteratorSize, HasShape, HasLength, SizeUnknown

@testset "DiscreteClock" begin
    c1 = DiscreteClock(3)
    c2 = DiscreteClock(0:3)
    c3 = DiscreteClock(0, 1:3)

    @testset "IteratorSize" begin
        @test IteratorSize(c1) == HasLength()
        @test IteratorSize(c2) == HasLength()
        @test IteratorSize(c3) == HasLength()
    end

    @testset "length" begin
        @test length(c1) == 3
        @test length(c2) == 3
        @test length(c3) == 3
    end

    @testset "eltype" begin
        @test eltype(c1) == Int
        @test eltype(c2) == Int
        @test eltype(c3) == Int
    end

    @testset "collect" begin
        @test collect(c1) == [1, 2, 3]
        @test collect(c2) == [1, 2, 3]
        @test collect(c3) == [1, 2, 3]
    end
end

@testset "ContinuousClock" begin
    c1 = ContinuousClock(3.0)
    c2 = ContinuousClock(3.0; max_epoch=2)
    c3 = ContinuousClock(3.0, 1.0)

    @testset "IteratorSize" begin
        @test IteratorSize(c1) == SizeUnknown()
        @test IteratorSize(c2) == SizeUnknown()
        @test IteratorSize(c3) == SizeUnknown()
    end

    @testset "eltype" begin
        @test eltype(c1) == Float64
        @test eltype(c2) == Float64
        @test eltype(c3) == Float64
    end

    @testset "comprehension" begin
        @test [(increase!(c1, 1); (i, currenttime(c1))) for i in c1] == collect(zip(1:3, 1.0:3.0))
        @test [(increase!(c2, 1); (i, currenttime(c2))) for i in c2] == collect(zip(1:2, 1.0:2.0))
        @test [(increase!(c3, 1); (i, currenttime(c3))) for i in c3] == collect(zip(1:2, 2.0:3.0))
    end

    @testset "auto-initialize" begin
        @test currenttime(c1) == start(c1)
        @test currenttime(c2) != start(c2)
        @test currenttime(c3) == start(c3)
    end
end


text:
export count_cross

function _orientation(a, b, c)
    A = ones(3, 3)
    A[1, 2:3] = a
    A[2, 2:3] = b
    A[3, 2:3] = c
    return sign(det(A))
end

function _cross(a, b, c, d)::Bool
    if _orientation(a, b, c) * _orientation(a, b, d) >= 0
        return false
    end
    if _orientation(c, d, a) * _orientation(c, d, b) >= 0
        return false
    end
    return true
end

"""
`count_cross(G::SimpleGraph)` returns the number of edge crossings
in the `xy`-embedding of this graph.
"""
function count_cross(G::SimpleGraph)::Int
    EE = elist(G)
    m = length(EE)
    count = 0
    xy = getxy(G)
    for i = 1:m-1
        e1 = EE[i]
        a = xy[e1[1]]
        b = xy[e1[2]]
        for j = i+1:m
            e2 = EE[j]
            if length(intersect(e1, e2)) == 0 # common end point doesn't count
                c = xy[e2[1]]
                d = xy[e2[2]]
                if _cross(a,b,c,d)
                    count += 1
                end
            end
        end

    end
    return count
end


text:
export
    Additive,
    Allocation,
    Conflicts,
    Constraint,
    Counts,
    OrderedCategory,
    Profile,
    Reduction,
    Submodular,
    agents,
    alloc_bb18_3,
    alloc_bkv18_1,
    alloc_bkv18_2,
    alloc_ef1,
    alloc_efx,
    alloc_ggi,
    alloc_ghss18_4,
    alloc_ghss18_4b,
    alloc_half_mms,
    alloc_hpps20_1,
    alloc_mm,
    alloc_mms,
    alloc_mnw,
    alloc_mnw_ef1,
    alloc_rand,
    bundle,
    category!,
    category,
    ceil_n,
    check,
    check_complete,
    check_ef,
    check_ef1,
    check_efx,
    check_partition,
    constraint,
    deny!,
    fill_even!,
    floor_n,
    give!,
    graph,
    items,
    mms,
    mms_alpha,
    na,
    nash_welfare,
    ni,
    normalize,
    owned,
    owner,
    profile,
    prop_alpha,
    required,
    threshold,
    utility,
    value,
    value!


text:
module CaseModule

using DataFrames

export
      Case,
      defineCase,
      setWeights!,
      setMisc!,
      setupResults,
      case2dfs,  # currently not working
      dataSet,
      Obs,
      defineObs
################################################################################
# Basic Types
################################################################################
############################### track info ########################################
type Track
  name
  func
  dir
  a
  b
  c
  y0
  Y
  X
end

function Track()
Track([],
      [],
      [],
      [],
      [],
      [],
      [],
      [],
      []);
end

function defineTrack(name)
  t=Track();
  if name==:path # test case for testPathFollowing.jl
    t.name=name;
    t.func=:poly;
    t.dir=:posY
    t.a=[0.0,0.0,2.33e-3,-6.43e-6,5.07e-9];
    t.Y=0:0.1:250;
    t.X=0:0.1:65;
  elseif name==:caseStudy  # TODO consider adding an @NLexpression for the path
    t.name=name;
    t.func=:fourier;
    t.dir=:posX
    dfs=dataSet(name,Pkg.dir("MAVs/examples/Cases/Tracks/"));
    t.a=dfs[1][:a][1:8];
    t.b=dfs[1][:b][1:8];
    t.c=dfs[1][:c][1:8];
    t.y0=dfs[1][:y0][1];
    t.Y=dfs[1][:Y];
    t.X=dfs[1][:X];
  end
  return t
end

############################### weights ########################################
type Weights
  goal  # should be zero when vehicle is not within goal distance during prediction horizon
  psi
  time
  haf
  Fz
  ce
  sa
  sr
  jx
  path
  driver
end

function Weights()
  Weights(1.,      # w_goal
          0.01,    # w_psi
          0.05,    # w_time
          1.0e-5,  # w_haf
          0.1,     # w_Fz  0.5 in paper
          1.,      # w_ce
          0.1,     # w_sa
          1.,      # w_sr
          0.01,    # w_jx
          1.0,     # path following
          1.0,     # follow driver steering angle
          );
end

function defineWeights(name)
  if name==:auto
    w=Weights();
  elseif name==:path
    w=Weights();
    w.sr=0.1;
    w.path=1.0;
    w.driver=1.0;
  elseif name==:m3

  elseif name!==:NA
    error("\n Pick a name for weights! \n")
  end
  return w
end

"""
setWeights!(c;sr=0.08,path=10.0,driver=0.5)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/12/2017, Last Modified: 4/12/2017 \n
--------------------------------------------------------------------------------------\n
"""
function setWeights!(c;
                      sr=c.w.sr,
                      path=c.w.path,
                      driver=c.w.driver
                      )
    c.w.sr=sr;
    c.w.path=path;
    c.w.driver=driver;
    return nothing
end
############################### obstacle info ########################################
type Obs
  name
  A
  B
  s_x
  s_y
  X0
  Y0
  status
end

function Obs()
  Obs([],
            [],
            [],
            [],
            [],
            [],
            [],
            []);
end

"""
o=defineObs(name)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 7/5/2017 \n
--------------------------------------------------------------------------------------\n
"""
function defineObs(name)
  o=Obs();
  if name==:auto
    o.name=name;
    o.A=[5.];
    o.B=[5.];
    o.s_x=[0.0];
    o.s_y=[0.0];
    o.X0=[200.];
    o.Y0=[75.];
    o.status=falses(length(o.X0));
  elseif name==:autoBench
    o.name=name;
    o.A=[10.];
    o.B=[5.];
    o.s_x=[0.];
    o.s_y=[0.];
    o.X0=[200.];
    o.Y0=[43.];
    o.status=falses(length(o.X0));
  elseif name==:autoGazebo
    o.name=name;
    o.A=[1.];
    o.B=[1.];
    o.s_x=[0.];
    o.s_y=[0.];
    o.X0=[200.];
    o.Y0=[50.];
    o.status=falses(length(o.X0));
  elseif name==:autoARC
    o.name=name;
    A=[2.5 1.5];      # small car
    B=[3.3 1.9];      # HUMMVEEs
    C=[9.8/2 3.65/2]; # tanks
    random=[0 0.6 1 0.25 0.6 0.5495 0.4852 0.8905 0.7990 1];
    o.X0=[225,225,220,170,220,165,200];
    o.Y0=[25,30,45,60,75,95,50];
    o.A=[B[1],B[1],B[1],C[1],A[1],C[1],C[2]];
    o.B=[B[2],B[2],B[2],C[2],A[2],C[2],C[1]];
    o.s_x=[-5,-5.5,-6,6,-5,4,-2];
    o.s_y=[0,0,0,0,0,0,5];
  elseif name==:path  # test case for testPathFollowing.jl
    o.name=name;
    o.A=[5.,4.,10.]
    o.B=[5.,5.,6.]
    o.s_x=[0.,0.,0.]
    o.s_y=[0.,0.,0.]
    o.X0=[7.0648272,17.37,39.98];
    o.Y0=[60.,130,173]
    o.status=falses(length(o.X0));
  elseif name==:o3  # test case for testSharedControl.jl
    o.name=name;
    o.A=[1.5,2.0,1.5]
    o.B=[1.5,2.0,1.5]
    o.s_x=[0.,0.,0.]
    o.s_y=[0.,0.,0.]
    o.X0=[0.0,-2.5,10.]
    o.Y0=[40.,30.,75.]
    o.status=falses(length(o.X0));
  elseif name==:o4  # test case for testPathFollowing.jl
    o.name=name;
    o.A=[5.,10.]
    o.B=[7.5,3.]
    o.s_x=[0.,0.]
    o.s_y=[0.,0.]
    o.X0=[17.,49.];
    o.Y0=[100.,200.];
    o.status=falses(length(o.X0));
  elseif name==:caseStudyPath # test case for testPathFollowing.jl
    o.name=:caseStudy;
    dfs=dataSet(o.name,Pkg.dir("MAVs/examples/Cases/Obstacles/"));
    o.X0=dfs[1][:X];
    o.Y0=dfs[1][:Y];
    o.A=dfs[1][:A];
    o.B=dfs[1][:B];
    o.status=dfs[1][:status];
    o.s_x=zeros(length(o.X0));
    o.s_y=zeros(length(o.X0));
    o.status=falses(length(o.X0));
  elseif name==:caseStudy # test case for testPathFollowing.jl
    o.name=:caseStudy;
    dfs=dataSet(o.name,Pkg.dir("MAVs/examples/Cases/Obstacles/"));
    o.X0=dfs[1][:X][1];
    o.Y0=dfs[1][:Y][1];
    o.A=dfs[1][:A][1];
    o.B=dfs[1][:B][1];
    o.status=dfs[1][:status][1];
    o.s_x=zeros(length(o.X0));
    o.s_y=zeros(length(o.X0));
    o.status=falses(length(o.X0));
  elseif name!==:NA
    error("\n Pick a name for obstacle data! \n")
  end
  return o
end

############################### goal info ########################################
type Goal
    name
    x_ref
    y_ref
    psi_ref
end

function Goal()
       Goal([],
            [],
            [],
            []);
end

"""
defineGoal(name)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""
function defineGoal(name)
  if name==:auto
    g=Goal();
    g.name=name;
    g.x_ref=200.;
    g.y_ref=125.;
    g.psi_ref=pi/2;
  elseif name==:autoBench
    g=Goal();
    g.name=name;
    g.x_ref=200.;
    g.y_ref=100.;
    g.psi_ref=pi/2;
  elseif name==:autoGazebo
    g=Goal();
    g.name=name;
    g.x_ref=200.;
    g.y_ref=100.;
    g.psi_ref=pi/2;
  elseif name==:path # test case for testPathFollowing.jl
    g=Goal();
    g.name=name;
    g.x_ref=65.0;
    g.y_ref=250.0;
    g.psi_ref=pi/2; #NA
  elseif name==:caseStudy
    g=Goal();
    g.name=name;
    g.x_ref=702.372760886434;
    g.y_ref=120.986699170193;
    g.psi_ref=pi/2;
  elseif name==:NA # g is already defined
    g=Goal();
    g.name=name;
  else
    error("\n Pick a name for goal data! \n")
  end
  return g
end

############################### misc info ########################################
type Misc
  name
  model              # function name
  X0                 # intial state (at the start of the simulation)
  Xlims              # limit for x
  Ylims              # limit for Y
  UX                 # vehicle speed (m/s)
  tp                 # prediction horizon (s)
  tex                # execution horizon (s)
  max_cpu_time       # maximum time the algorithm has
  sm                 # (m) distance to make sure we don't hit obstacle
  Lr                 # LiDAR range (m)
  L_rd               # relaxation distance to LiDAR range
  sigma              # 0.05 (m)small margin, if the vehicle is within this margin, then the target is considered to be reached
  Nck                # number of points per interval
  solver             # either :IPOPT or :KNITRO
  max_iter           # max evaluations in optimization
  mpc_max_iter       # an MPC parameter
  PredictX0          # if the state that the vehicle will be at when the optimization is finished is to be predicted
  FixedTp            # fixed final time for predicting X0
  activeSafety       # a bool if the system is only used to avoid collisions with obstacles
  followDriver       # a bool to try and follow the driver
  followPath         # a bool to follow the path
  NF                 # number to subtract off the constraint check vector-> used for activeSafety mode so driver sees obstacle before autonomy
  integrationScheme
end

function Misc()
       Misc([],
            [],
            [],
            [],
            [],
            0.0,
            0.0,
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            [],
            []
            );
end

"""
defineMisc(name)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 4/3/2017 \n
--------------------------------------------------------------------------------------\n
"""
function defineMisc(name)
  m=Misc();
  if name==:auto
    m.name=name;
    m.model=:ThreeDOFv2;
    m.X0=[200.0, 0.0, 0.0, 0.0, pi/2,0.0,7.0,0.0];
    m.Xlims=[-1., 400.]
    m.Ylims=[-1., 400.]
    m.tex=0.5;
    m.max_cpu_time=m.tex;
    m.sm=2.0;
    m.sigma=1.0;
    m.Nck=[12,10,8,6];
    m.solver=:KNITRO;
    m.max_iter=500;
    m.mpc_max_iter=600;
    m.PredictX0=true;
    m.FixedTp=true;
    m.Lr=50.;
    m.L_rd=5.;
    m.integrationScheme=:lgrExplicit
  elseif name==:autoBench
    m.name=name;
    m.model=:ThreeDOFv2;
    m.X0=[200.0, 0.0, 0.0, 0.0, pi/2,0.0,15.0,0.0];
    m.Xlims=[111.,250.]
    m.Ylims=[-1., 140.]
    m.tex=0.5;
    m.max_cpu_time=0.47;#0.41;
    m.sm=5.0;
    m.sigma=1.0;
    m.Nck=[10,8,6];#[12,10,8,6];
    m.solver=:Ipopt;
    m.max_iter=500;
    m.mpc_max_iter=60;
    m.PredictX0=true;
    m.FixedTp=true;
    m.Lr=150. # not in play
    m.L_rd=5.;
    m.integrationScheme=:lgrExplicit
  elseif name==:autoGazebo
    m.name=name;
    m.model=:ThreeDOFv2;
    m.X0=[200.0, 0.0, 0.0, 0.0, pi/2, 0.0, 0.0, 0.0];
    m.Xlims=[111.,250.]
    m.Ylims=[-1., 140.]
    m.tex=0.5;
    m.max_cpu_time=0.47;#0.41;
    m.sm=5.0;
    m.sigma=1.0;
    m.Nck=[10,8,6];#[12,10,8,6];
    m.solver=:KNITRO;
    m.max_iter=500;
    m.mpc_max_iter=60;
    m.PredictX0=true;
    m.FixedTp=true;
    m.Lr= 50. # not in play
    m.L_rd=5.;
    m.integrationScheme=:lgrExplicit
  elseif name==:autoARC
    m.name=name;
    m.model=:ThreeDOFv2;
    m.X0=[200.0, 0.0, 0.0, 0.0, pi/2,0.0,17.0,0.0];
    m.Xlims=[111.,250.]
    m.Ylims=[-1., 140.]
    m.tex=0.5;
    m.max_cpu_time=0.47;#0.41;
    m.sm=2.6;
    m.sigma=1.0;
    m.Nck=[10,8,6];#[12,10,8,6];
    m.solver=:KNITRO;
    m.max_iter=500;
    m.mpc_max_iter=600;
    m.PredictX0=true;
    m.FixedTp=true;
    m.Lr=50.
    m.L_rd=5.;
    m.integrationScheme=:lgrExplicit
  elseif name==:path
    m.name=name;
    m.model=:ThreeDOFv2;
    m.UX=10.0;
    m.X0=[0.0, 0.0, 0.0, 0.0, pi/2,0.0,m.UX,0.0];
    m.Xlims=[-1., 100.]
    m.Ylims=[-1., 300.]
    m.tp=6.0;
    m.tex=0.5;
    m.max_cpu_time=m.tex;
    m.sm=2.6;
    m.sigma=1.0;
    m.Nck=[10,8,6];
    m.solver=:KNITRO;
    m.max_iter=500;
    m.mpc_max_iter=30;
    m.PredictX0=true;
    m.FixedTp=true;
    m.integrationScheme=:lgrExplicit
  elseif name==:caseStudy
    m.name=name;
    m.model=:ThreeDOFv2;
    m.UX=10.0;
    m.X0=[0.0, 0.0, 0.0, 0.0, 1.2037,0.0,m.UX,0.0];
    m.Xlims=[-10., 750.]
    m.Ylims=[-10., 200.]
    m.tp=6.0;
    m.tex=0.3;
    m.max_cpu_time=0.25;  #NOTE:this needs to match Matlab model!
    m.sm=2.0;
    m.sigma=1.0;
    m.Nck=[10,8,6];
    m.solver=:KNITRO;
    m.max_iter=450;
    m.mpc_max_iter=600;
    m.PredictX0=true;
    m.FixedTp=true;
    m.activeSafety=true;
    m.followPath=false;
    m.followDriver=false;
    m.NF=0;
    m.integrationScheme=:lgrExplicit
  elseif name==:caseStudyPath # test case for testPathFollowing.jl with other model
    m.name=:caseStudy;
    m.model=:ThreeDOFv2;
    m.UX=10.0;
    m.X0=[0.0,0.0, 0.0, 0.0,1.2037,0.0,m.UX,0.0];
    m.Xlims=[-10., 730.]
    m.Ylims=[-10., 200.]
    m.tp=7.0;
    m.tex=0.3;
    m.max_cpu_time=0.25;
    m.sm=2.0;
    m.sigma=1.0;
    m.Nck=[10,8,6];
    m.solver=:KNITRO;
    m.max_iter=300;
    m.mpc_max_iter=30;
    m.PredictX0=true;
    m.FixedTp=true;
    m.followPath=true;
    m.integrationScheme=:lgrExplicit
  elseif name!==:NA
    error("\n Pick a name for misc data! \n")
  end
  return m
end

"""
setMisc!(c;activeSafety=true,followPath=false,followDriver=false,predictX0=false,fixedTp=false,tp=5.0,tex=0.3,max_cpu_time=0.26,Ni=3,Nck=[10,8,6]);
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/28/2017, Last Modified: 4/7/2017 \n
--------------------------------------------------------------------------------------\n
"""
function setMisc!(c;
                Xlims=c.m.Xlims,
                Ylims=c.m.Ylims,
                tp=c.m.tp,
                tex=c.m.tex,
                max_cpu_time=c.m.max_cpu_time,
                sm=c.m.sm,
                Lr=c.m.Lr,
                L_rd=c.m.L_rd,
                sigma=c.m.sigma,
                Nck=c.m.Nck,
                solver=c.m.solver,
                max_iter=c.m.max_iter,
                mpc_max_iter=c.m.mpc_max_iter,
                PredictX0=c.m.PredictX0,
                FixedTp=c.m.FixedTp,
                activeSafety=c.m.activeSafety,
                followPath=c.m.followPath,
                followDriver=c.m.followDriver,
                NF=c.m.NF,
                integrationScheme=c.m.integrationScheme);
    c.m.Xlims=Xlims;
    c.m.Ylims=Ylims;
    c.m.tp=tp;
    c.m.tex=tex;
    c.m.max_cpu_time=max_cpu_time;
    c.m.sm=sm;
    c.m.Lr=Lr;
    c.m.L_rd=L_rd;
    c.m.sigma=sigma;
    c.m.Nck=Nck;
    c.m.solver=solver;
    c.m.max_iter=max_iter;
    c.m.mpc_max_iter=mpc_max_iter;
    c.m.PredictX0=PredictX0;
    c.m.FixedTp=FixedTp;
    c.m.activeSafety=activeSafety;
    c.m.followPath=followPath;
    c.m.followDriver=followDriver;
    c.m.NF=NF;
    c.m.integrationScheme=integrationScheme;
    return nothing
end

################################################################################
# Model Class
################################################################################
abstract type AbstractCase end
type Case <: AbstractCase
 name
 g::Goal        # goal data
 o::Obs   # obstacle data
 w::Weights     # weight data
 t::Track       # track data
 m::Misc        # miscelaneous data
end

"""
c = Case();
# default constructor
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/11/2017, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""
function Case()
 Case(Any,
      Goal(),
      Obs(),
      Weights(),
      Track(),
      Misc()
    );
end

"""
# a basic example of autonomousControl
c=defineCase(;(mode=>:auto));

# a basic example of PathFollowing
c=defineCase(;(mode=>:path));

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/11/2017, Last Modified: 7/4/2017 \n
--------------------------------------------------------------------------------------\n
"""
function defineCase(;name::Symbol=:auto,
                    goal::Symbol=:auto,
                    obstacles::Symbol=:auto,
                    weights::Symbol=:auto,
                    track::Symbol=:NA,
                    misc::Symbol=:auto, kwargs...)

  kw = Dict(kwargs);
  if !haskey(kw,:mode); kw_ = Dict(:mode => :auto); mode = get(kw_,:mode,0);
  else; mode = get(kw,:mode,0);
  end

   c=Case();
   if mode==:auto
     c.name=mode;
     c.g=defineGoal(c.name);
     c.o=defineObs(c.name);
     c.w=defineWeights(c.name);
     c.t=defineTrack(:NA);
     c.m=defineMisc(c.name);
   elseif mode==:autoBench
     c.name=mode;
     c.g=defineGoal(c.name);
     c.o=defineObs(c.name);
     c.w=defineWeights(:auto);
     c.t=defineTrack(:NA);
     c.m=defineMisc(c.name);
   elseif mode==:autoGazebo
     c.name=mode;
     c.g=defineGoal(c.name);
     c.o=defineObs(c.name);
     c.w=defineWeights(:auto);
     c.t=defineTrack(:NA);
     c.m=defineMisc(c.name);
   elseif mode==:autoARC
     c.name=:auto;
     c.g=defineGoal(c.name);
     c.o=defineObs(:autoARC);
     c.w=defineWeights(c.name);
     c.t=defineTrack(:NA);
     c.m=defineMisc(:autoARC);
   elseif mode==:path
     c.name=mode;
     c.g=defineGoal(c.name);
     c.o=defineObs(c.name);
     c.w=defineWeights(c.name);
     c.t=defineTrack(c.name);
     c.m=defineMisc(c.name);
   elseif mode==:caseStudy
     c.name=mode;
     c.g=defineGoal(c.name);
     c.o=defineObs(c.name);
     c.w=defineWeights(:path);
     c.t=defineTrack(c.name);
     c.m=defineMisc(c.name);
   elseif mode==:caseStudyPath
     c.name=:caseStudy;
     c.g=defineGoal(c.name);
     c.o=defineObs(:caseStudyPath);
     c.w=defineWeights(:path);
     c.t=defineTrack(c.name);
     c.m=defineMisc(:caseStudyPath);
   else
     name=:user;
     c.name=name;
     c.g=defineGoal(goal);
     c.o=defineObs(obstacles);
     c.w=defineWeights(weights);
     c.t=defineTrack(:NA);
     c.m=defineMisc(misc);
   end

   return c
end

################################################################################
# data stuff
################################################################################
"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/7/2016, Last Modified: 3/28/2017 \n
--------------------------------------------------------------------------------------\n
"""
function dataSet(set,path)
  main_dir=pwd();
  dfs=Vector{DataFrame}(1) # create am empty DataFrame
  cd(path)
  dfs[1]=readtable(string(set,".csv"))
  cd(main_dir)
  return dfs
end

"""
case2dfs(c)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 3/11/2017, Last Modified: 3/11/2017 \n
--------------------------------------------------------------------------------------\n
"""
# save senario data  # TODO add all options and settings here?
function case2dfs(c)
# TODO think about just saving n,r,s,params, and mdl
 ###########
 # goal data
 x_ref_data = DataFrame(ID=1, x_ref = c.x_ref)
 y_ref_data = DataFrame(ID=1, y_ref = c.y_ref)
 psi_ref_data = DataFrame(ID=1, psi_ref = c.psi_ref)

 s_data = join(x_ref_data,y_ref_data,on=:ID)
 s_data = join(s_data,psi_ref_data,on=:ID)

 ###############
 # obstacle data
 A_data = DataFrame(ID = 1:Q, A = c.A);
 B_data = DataFrame(ID = 1:Q, B = c.B);
 s_x_data = DataFrame(ID = 1:Q, s_x = c.s_x);
 s_y_data = DataFrame(ID = 1:Q, s_y = c.s_y);
 X0_data = DataFrame(ID = 1:Q, X0 = c.X0_obs);
 Y0_data = DataFrame(ID = 1:Q, Y0 = c.Y0_obs);

 obs_data = join(A_data,B_data, on = :ID);
 obs_data = join(obs_data,s_x_data,on =:ID);
 obs_data = join(obs_data,s_y_data,on =:ID);
 obs_data = join(obs_data,X0_data,on =:ID);
 obs_data = join(obs_data,Y0_data,on=:ID);

 return case_data, obs_data
end

end # module


text:
macro workunit(title, code, savelist)
    expr = quote
        is_cached = all(i->isfile(i), [$((i.args[3] for i in savelist.args)...)])

        is_cached && println("SKIPPING: ", $title)

        is_cached && return

        println("STARTING: ", $title)

        $code

        println("SAVING: ", $title)

        map(i->save(i[2], i[1]), [$((:($(i.args[2]), $(i.args[3])) for i in savelist.args)...)])

        println("FINISHED: ", $title)
    end

    thunk = esc(:(()->($expr)))
    var = esc(Base.sync_varname)
    spawncall = :($(Distributed.spawn_somewhere)($thunk))
    quote
        local ref = $spawncall
        if $(Expr(:islocal, var))
            put!($var, ref)
        end
        ref
    end
end


text:
using BayesLaplace
using BayesLaplace.Curvature
using BayesLaplace.Data
using Flux

@testset "Construction" begin

    # One layer:
    nn = Chain(Dense(2,1))
    la = Laplace(nn)
    @test la.n_params == 3

    # Multi-layer:
    nn = Chain(Dense(2,2,σ),Dense(2,1))
    la = Laplace(nn)
    @test la.n_params == 3
    la = Laplace(nn; subset_of_weights=:all)
    @test la.n_params == 9

end

@testset "Parameters" begin

    nn = Chain(Dense(2,2,σ),Dense(2,1))
    la = Laplace(nn)
    @test BayesLaplace.get_params(la) == collect(Flux.params(nn))[(end-1):end]

    la = Laplace(nn; subset_of_weights=:all)
    @test BayesLaplace.get_params(la) == collect(Flux.params(nn))

end

# We know the analytical expression for the gradient of logit binary cross entropy loss for a single-layer neural net with sigmoid activation just corresponds to the gradient in logistic regression (see for example: https://www.paltmeyer.com/blog/posts/bayesian-logit/): ∇ℓ=(μ-y)x. We can use this analytical expression to see if we get the expected results.

@testset "Hessian" begin

    nn = Chain(Dense([0 0]))
    la = Laplace(nn)

    # (always ignoring constant)
    @testset "Empirical Fisher - full" begin

        target = 1
        x = [1,1]
        grad = [-0.5,-0.5] # analytical solution for gradient
        @test BayesLaplace.hessian_approximation(la, (x,target))[1:2,1:2] == grad * grad'

        x = [-1,-1]
        grad = [0.5,0.5] # analytical solution for gradient
        @test BayesLaplace.hessian_approximation(la, (x,target))[1:2,1:2] == grad * grad'

        x = [0,0]
        grad = [0,0] # analytical solution for gradient
        @test BayesLaplace.hessian_approximation(la, (x,target))[1:2,1:2] == grad * grad'

    end

end

# We know the analytical expression for the Hessian of logit binary cross entropy loss for a single-layer neural net with sigmoid activation just corresponds to the Hessian in logistic regression (see for example: https://www.paltmeyer.com/blog/posts/bayesian-logit/): ∇ℓ=(μ-y)(μ(1-μ)xx'). With a weight-penalty (Gaussian prior), the Hessian becomes: ∇ℓ=∑(μ-y)(μ(1-μ)xx')+H₀. We can use this analytical expression to see if we get the expected results.

@testset "Fitting" begin

    nn = Chain(Dense([0 0]))
    la = Laplace(nn)

    hessian_exact(x,target) = (nn(x).-target).*(nn(x).*(1 .- nn(x)).*x*x') + la.H₀

    @testset "Empirical Fisher - full" begin

        target = [1]
        x = [[0,0]]
        fit!(la,zip(x,target))
        @test la.H[1:2,1:2] == hessian_exact(x[1],target[1])

    end

end