text:
# Advent of Code 2020 - Day 20

using Combinatorics
using Logging

test_str = "Tile 2311:
..##.#..#.
##..#.....
#...##..#.
####.#...#
##.##.###.
##...#.###
.#.#.#..##
..#....#..
###...#.#.
..###..###

Tile 1951:
#.##...##.
#.####...#
.....#..##
#...######
.##.#....#
.###.#####
###.##.##.
.###....#.
..#.#..#.#
#...##.#..

Tile 1171:
####...##.
#..##.#..#
##.#..#.#.
.###.####.
..###.####
.##....##.
.#...####.
#.##.####.
####..#...
.....##...

Tile 1427:
###.##.#..
.#..#.##..
.#.##.#..#
#.#.#.##.#
....#...##
...##..##.
...#.#####
.#.####.#.
..#..###.#
..##.#..#.

Tile 1489:
##.#.#....
..##...#..
.##..##...
..#...#...
#####...#.
#..#.#.#.#
...#.#.#..
##.#...##.
..##.##.##
###.##.#..

Tile 2473:
#....####.
#..#.##...
#.##..#...
######.#.#
.#...#.#.#
.#########
.###.#..#.
########.#
##...##.#.
..###.#.#.

Tile 2971:
..#.#....#
#...###...
#.#.###...
##.##..#..
.#####..##
.#..####.#
#..#.#..#.
..####.###
..#.#.###.
...#.#.#.#

Tile 2729:
...#.#.#.#
####.#....
..#.#.....
....#..#.#
.##..##.#.
.#.####...
####.#.#..
##.####...
##..#.##..
#.##...##.

Tile 3079:
#.#.#####.
.#..######
..#.......
######....
####.#..#.
.#...#.##.
#.#####.##
..#.###...
..#.......
..#.###..."

test_list = split(test_str, '\n')

list = readlines("day20_input.txt")

"""
The operators are not always commutative, but when looking for all combinations,
it suffices to fix the order. Thus, in all what follows, it is assumed that rotations
are applied first, and then, flips.

In fact:
    * Commutative:
        * R0∘FH = FH∘R0 (R0 = identity)
        * R0∘FV = FV∘R0 (R0 = identity)
        * R180∘FH = FH∘R180
        * R180∘FV = FV∘R180
    * Non-commutative:
        * R90*FH = FV∘R90 = FH*R270
        * R90∘FV = FH∘R90 = FV∘R270
        * R270∘FH = FV∘R270 = FH∘R90
        * R270∘FV = FH∘R270 = FV∘R90

Notice from the above that there are some redundancies, even fixing the order:
    * FH∘R0 = FV∘R180
    * FH∘R90 = FV∘R270
    * FH∘R180 = FV∘R0
    * FH∘R270 = FV∘R90

The last three also follow from the first one by repeatedly applying R90 from the right.

From these, we see that we may consider only FH or only FV.

It is easy to deduce the relations above imagining a 2D square, but it is even
easier in the complex plane, with z=x+iy:
    * R0 = 1
    * R90 = multiplication by i
    * R180 = i^2 = -1
    * R270 = i^3 = -i
    * FN = 1
    * FV = conjugate (z -> z̄ )
    * FH = z -> -z̄ (in fact z = x+iy -> -x+iy = -(x-iy) = -z̄ ) (i.e. FH∘R0 = FV∘R180)
"""
const R0, R90, R180, R270 = Int8.((1, 2, 3, 4))
const FN, FH, FV = Int8.((1, 2, 3))
const RIGHT, LEFT, TOP, BOTTOM = Int8.((1, 2, 3, 4))

function get_side(tile, rotation, flip, side)
    """
        get_side(tile, orientation, side)

    Return the desired `side` for a given `tile` with the given `orientation`
    and `flip`

    The sides are, of course, `RIGHT`, `LEFT`, `TOP`, and `BOTTOM`

    The rotations are `R0`, `R90`, `R180`, and `R270` for,
    respectively 0ᵒ, 90ᵒ, 180ᵒ, and 270ᵒ anticlockwise rotations

    The possible flips are `FN`, `FH`, and `FV` for no-flip and horizontal
    and vertical flips, respectively.

    The returned side is from left to right, for `top` and `bottom` sides,
    and from top to bottom, for `LEFT` and `RIGHT` sides.
    """
    if (rotation, flip, side) ∈ ((R0, FN, RIGHT), (R0, FH, LEFT), (R90, FN, TOP), (R90, FV, BOTTOM), (R180, FV, LEFT), (R270, FH, BOTTOM))
        return tile[:,end]
    elseif (rotation, flip, side) ∈ ((R0, FV, RIGHT), (R90, FH, TOP), (R180, FN, LEFT), (R180, FH, RIGHT), (R270, FN, BOTTOM), (R270, FV, TOP))
        return tile[end:-1:1,end]
    elseif (rotation, flip, side) ∈ ((R0, FN, TOP), (R0, FV, BOTTOM), (R90, FV, LEFT), (R180, FH, BOTTOM), (R270, FN, RIGHT), (R270, FH, LEFT))
        return tile[1,:]
    elseif (rotation, flip, side) ∈ ((R0, FH, TOP), (R90, FN, LEFT), (R90, FH, RIGHT), (R180, FN, BOTTOM), (R180, FV, TOP), (R270, FV, RIGHT))
        return tile[1,end:-1:1]
    elseif (rotation, flip, side) ∈ ((R0, FN, LEFT), (R0, FH, RIGHT), (R90, FN, BOTTOM), (R90, FV, TOP), (R180, FV, RIGHT), (R270, FH, TOP))
        return tile[:,1]
    elseif (rotation, flip, side) ∈ ((R0, FV, LEFT), (R90, FH, BOTTOM), (R180, FN, RIGHT), (R180, FH, LEFT), (R270, FN, TOP), (R270, FV, BOTTOM))
        return tile[end:-1:1,1]
    elseif (rotation, flip, side) ∈ ((R0, FN, BOTTOM), (R0, FV, TOP), (R90, FV, RIGHT), (R180, FH, TOP), (R270, FN, LEFT), (R270, FH, RIGHT))
        return tile[end,:]
    elseif (rotation, flip, side) ∈ ((R0, FH, BOTTOM), (R90, FN, RIGHT), (R90, FH, LEFT), (R180, FN, TOP), (R180, FV, BOTTOM), (R270, FV, LEFT))
        return tile[end, end:-1:1]
    end
end

"""
    get_encoded_sides(tile)

Return an array with all the integer-encoded sides after all the possible transformations.

The sides of the tile are given as an array with dimensions 12 x 4,
representing all 12 transformations (4 rotations x 3 flips) and all 4 sides,
despite the fact that some transformations are redundant.

Each side is encoded into an integer, i.e. the decimal representation of
the binary number obtained when translating `"#" => true, "."` => false`.
"""
function get_encoded_sides(tile)
    sides = fill(Int16(0),4,3,4) # 4 rotations, 3 flips, 4 sides
    for s=1:4, f=1:3, r=1:4 # s=side, f=flip, r=rotation
        side = get_side(tile, r, f, s)
        sides[r,f,s] = [2^(k-1) for k=1:length(side)]' * side
    end
    return sides
end

"""
    get_tiles(list)

Return array with a binary representation of all the tiles in the `list`.

The translation of each tile is done according to `"#" => true, "."` => false`.
"""
function get_tiles(list)
    tile_side_length = findfirst(==(""), list) - 2 # discount header and empty line
    num_tiles = div(length(list)+1, tile_side_length+2) # add one since last tile is not followed by empty line
    tiles = Dict()
    tile = fill(false, tile_side_length, tile_side_length)
    for i in 1:num_tiles
        for j=1:tile_side_length
            tile[j,:] = map(==('#'), collect(list[(tile_side_length+2)*(i-1)+1+j]))
        end
        push!(tiles, parse(Int, match(r"Tile (\d*):", list[(tile_side_length+2)*(i-1) + 1])[1]) => copy(tile))
    end
    return tiles
end

"""
    get_sides_tiles(list)

Return a Dict with the array of sides of all the tiles.

The sides of each tile are assembled into an array with dimensions 12 x 4,
representing all 12 transformations (4 rotations x 3 flips) and all 4 sides,
despite reduntant transformations.

Each side is encoded into an integer, i.e. the decimal representation of
the binary number obtained when translating `"#" => true, "."` => false`.
"""
function get_sides_tiles(list)
    tile_side_length = findfirst(==(""), list) - 2 # discount header and empty line
    num_tiles = div(length(list)+1, tile_side_length+2) # add one since last tile is not followed by empty line
    sides_tiles = Dict()
    tile = fill(false, tile_side_length, tile_side_length)

    for n in 1:num_tiles
        for i=1:tile_side_length
            tile[i,:] = map(==('#'), collect(list[(tile_side_length+2)*(n-1)+1+i]))
        end
        push!(sides_tiles, parse(Int, match(r"Tile (\d*):",
            list[(tile_side_length+2)*(n-1) + 1])[1]) => get_encoded_sides(tile))
    end
    return sides_tiles
end

"""
solve_jigsaw_combinatorics(list, flipstep=1)

Solve jigsaw by finding all possible horizontal matches and matching them vertically.

It is very expensive computationally. It is feasible for the 3x3 test puzzle
but not for the 12x12 one.

With `flipstep=2`, it skips the redundant horizontal flip `FH`, while
if `flipstep=1`, it considers `FH` and `FV`, besides no flip `FN`.
"""
function solve_jigsaw_combinatorics(list, flipstep=1)
    sides_tiles = get_sides_tiles(list)
    num_tiles = length(sides_tiles)
    side_length = isqrt(num_tiles)

    horizontal_matches = Set()
    num_matches = 0
    for c in combinations(collect(keys(sides_tiles)), side_length)
        for p in permutations(c)
            for j in 0:4^side_length-1 # vary orientations for each tile
                r = digits(j, base=4, pad=side_length) .+ fill(1,side_length)
                for l in 0:flipstep:3^side_length-1 # vary flips for each tile, if flipstep=2, skip FH=R180∘FV=FV∘R180, it should find some suitable solutions
                    f = digits(l, base=3, pad=side_length) .+ fill(1,side_length)
                    matched = true
                    for k in 1:side_length-1
                        if sides_tiles[p[k]][r[k],f[k], 1] !=
                                sides_tiles[p[k+1]][r[k+1],f[k+1], 2]
                            matched = false
                            break
                        end
                    end
                    if matched
                        push!(horizontal_matches, collect(zip(p, r, f)))
                    end
                end
            end
        end
    end

    solutions = Set()
    for c in combinations(unique(horizontal_matches), side_length)
        if length(Set([c[i][j][1] for i=1:3,j=1:3])) == num_tiles
            for p in permutations(c)
                matched = true
                for j in 1:side_length-1, k in 1:side_length
                    if sides_tiles[p[j][k][1]][p[j][k][2],p[j][k][3],4] !=
                        sides_tiles[p[j+1][k][1]][p[j+1][k][2],p[j+1][k][3],3]
                        matched = false
                        break
                    end
                end
                if matched
                    push!(solutions, p)
                end
            end
        end
    end

    if length(solutions) == 0
        return nothing, nothing
    else
        corners_product =
            prod([first(solutions)[i][j][1] for i=[1;side_length], j=[1;side_length]])
        return solutions, corners_product
    end
end

"""
attach_tile(n, tableau, tile_key, sides_of_tiles, side_length, flipstep)

Return a new set of tableaux with all possible attachments of given tile to the given tableau.

Look for all possible attachments of tile `sides_of_tiles[tile_key]` to the
given `tableau`.

If `flipstep=2`, it skips the redundant horizontal flip `FH` of the given tile,
while if `flipstep=1`, all `FN`, `FH` and `FV` are considered.

`n` is the number of tiles already attached to the tableau. It could
"""
function attach_tile(n, tableau, tile_key, sides_of_tiles, side_length, flipstep)
    tile_sides = sides_of_tiles[tile_key]
    new_tableaux = Set()

    n = count(!=((0,0,0)), tableau) + 1
    side_length = size(tableau)[1]

    i, j = divrem(n-1,side_length) .+ (1,1)
    for r=1:4, f=1:flipstep:3 # rotations and FN and FV flips (if flipstep = 2, skip FH since FH = FV∘R180)
        side_matches = 0
        for (di, dj, s1, s2) in ((-1, 0, TOP, BOTTOM), (0, -1, LEFT, RIGHT))
            if 1 ≤ i+di ≤ side_length && 1 ≤ j+dj ≤ side_length && tableau[i+di,j+dj] != (0,0,0)
                if tile_sides[r,f,s1] == sides_of_tiles[tableau[i+di,j+dj][1]][tableau[i+di,j+dj][2],tableau[i+di,j+dj][3],s2]
                    side_matches += 1
                end
            else
                side_matches += 1
            end
        end
        if side_matches == 2
            new_tableau = copy(tableau)
            new_tableau[i,j] = (tile_key,r,f)
            push!(new_tableaux, new_tableau)
        end
    end

    return new_tableaux
end

function solve_jigsaw(list, flipstep=2, all=false)
    sides_of_tiles = get_sides_tiles(list)
    num_tiles = length(sides_of_tiles)
    side_length = isqrt(num_tiles)

    tableaux = Set()
    if all # This find all solutions accounting for symmetries (eight total)
        rrange = 1:4
        frange = 1:flipstep:3 # it suffices to consider FV, may skip FH, which is FV∘R180
    else # By symmetry, it suffices to put the 1st tile in the original position to find a solution
        rrange = 1:1
        frange = 1:1
    end
    for tk in keys(sides_of_tiles)
        for r ∈ rrange, f ∈ frange
            tableau = fill((0,0,0), side_length, side_length)
            tableau[1, 1] = (tk, r, f)
            push!(tableaux, tableau)
        end
    end
    for n in 2:num_tiles
        new_tableaux = Set()
        for tableau in tableaux
            for tk in keys(sides_of_tiles)
                if tk ∉ Set([k for (k,o) in tableau])
                    union!(new_tableaux, attach_tile(n, tableau, tk, sides_of_tiles, side_length, flipstep))
                end
            end
        end
        tableaux = copy(new_tableaux)
    end

    if length(tableaux) == 0
        return nothing, nothing
    else
        corners_product =
            prod([first(tableaux)[i,j][1] for i=[1;side_length], j=[1;side_length]])
        return tableaux, corners_product
    end
end

# solution given in AoC's page
aoc_test_sol = [
    (1951,R0,FV) (2311,R0,FV) (3079,R0,FN);
    (2729,R0,FV) (1427,R0,FV) (2473,R270,FV);
    (2971,R0,FV) (1489,R0,FV) (1171,R180,FV);
]

# AoC solution flipped vertically
aoc_test_sol_FV = [
    (2971,R0,FN) (1489,R0,FN) (1171,R180,FN);
    (2729,R0,FN) (1427,R0,FN) (2473,R270,FN);
    (1951,R0,FN) (2311,R0,FN) (3079,R0,FV);
]

"""
    all_prod_tree(list, flipstep=2)

Return a set with the product of the corners of all solutions.

Aim to check whether all solutions yield the same corner's product.
"""
function all_prod_tree(list, flipstep=2)
    tableaux, = solve_jigsaw(list, flipstep)
    prods = Set()
    side_length = isqrt(length(first(tableaux)))
    for tb in tableaux
        corners_product =
            prod([tb[i,j][1] for i=[1;side_length], j=[1;side_length]])
        push!(prods, corners_product)
    end
    return prods
end

@show solve_jigsaw_combinatorics(test_list, 1)[2] == 20899048083289
@show solve_jigsaw(test_list, 1)[2] == 20899048083289
@show aoc_test_sol_FV ∈ solve_jigsaw(test_list, 1)[1]
@show aoc_test_sol ∈ solve_jigsaw(test_list, 1, true)[1]
@show length(solve_jigsaw(test_list, 2, true)[1]) == 8
@show solve_jigsaw(list, 2)[2] == 22878471088273

#=
```julia
julia> @btime solve_jigsaw_combinatorics(test_list, 1); # no skipping FH
  983.581 ms (8762301 allocations: 498.44 MiB)

julia> @btime solve_jigsaw_combinatorics(test_list, 2); # skipping FH
  181.198 ms (2170275 allocations: 138.39 MiB)

julia> @btime solve_jigsaw(test_list, 1); # no skipping FH
  803.124 μs (8118 allocations: 505.52 KiB)

julia> @btime solve_jigsaw(test_list, 2); # skipping FH
  626.993 μs (6072 allocations: 443.05 KiB)

julia> @btime solve_jigsaw(list, 2);
  1.312 s (6099477 allocations: 1.05 GiB)
``` =#

# Part 2

# One does not need to find all solutions with `all=true` in `solve_jigsaw(list, flipstep=2, all=true)`,
# it suffices (it is faster) to find just one and then transform this one (i.e. transform
# the whole jigsaw.

"""
    transform(tile, rotation, flip)

First rotate and then flip a `tile`, with the given `rotation` and `flip` arguments.
"""
function transform(tile, rotation, flip)
    side_length = isqrt(length(tile))

    if rotation == R0
        rotated_tile = copy(tile)
    elseif rotation == R90
        rotated_tile = [tile[j, i] for i=side_length:-1:1, j=1:side_length]
    elseif rotation == R180
        rotated_tile = [tile[i, j] for i=side_length:-1:1, j=side_length:-1:1]
    elseif rotation == R270
        rotated_tile = [tile[j, i] for i=1:side_length, j=side_length:-1:1]
    end

    if flip == FN
        transformed_tile = copy(rotated_tile)
    elseif flip == FH
        transformed_tile = [rotated_tile[i, j] for i=1:side_length, j=side_length:-1:1]
    elseif flip == FV
        transformed_tile = [rotated_tile[i, j] for i=side_length:-1:1, j=1:side_length]
    end
    return transformed_tile
end

"""
    puzzle_roughness(list)

Return the water roughness of the puzzle assembled from the `list` of tiles.
"""
function puzzle_roughness(list)
    tiles = get_tiles(list)
    num_tiles = length(tiles)
    isqrt_num_tiles = isqrt(num_tiles)
    side_length = size(first(tiles)[2])[1]

    monster_str = split(
"                  #
#    ##    ##    ###
 #  #  #  #  #  #   ", '\n')
    monster_bitarray = map(x->x=='#', [monster_str[i][j] for i=1:3, j=1:20])

    solution = first(solve_jigsaw(list, 2, false)[1])
    jigsaw =
        vcat(
            [
                hcat(
                    [transform(tiles[solution[i,j][1]], solution[i,j][2], solution[i,j][3])[2:end-1,2:end-1]
                        for j=1:isqrt_num_tiles]...
                )
            for i=1:isqrt_num_tiles]...
        )

    isqrt_len_jigsaw = isqrt_num_tiles*(side_length-2)
    len_jigsaw = isqrt_len_jigsaw^2

    for r=1:4, f=1:2:3
        tj = transform(jigsaw, r, f)
        tj_stretched = reshape(tj, (len_jigsaw, 1))
        roughness = sum(tj)

        for j=1:isqrt_len_jigsaw-19, i=1:isqrt_len_jigsaw-2
            if tj[i:i+2, j:j+19].*monster_bitarray == monster_bitarray
                roughness -= sum((monster_bitarray) .* tj[i:i+2, j:j+19])
            end
        end
        if roughness < sum(tj)
            return roughness
        end
    end
    return nothing
end

@show puzzle_roughness(test_list) == 273
@show puzzle_roughness(list) == 1680

aoc_test_sol_mm = [ # = aoc_test_sol_FV_R90 position of the test puzzle found with the sea monster
    (1951,R90,FN) (2729,R270,FN) (2971,R270,FN);
    (2311,R270,FN) (1427,R270,FN) (1489,R270,FN);
    (3079,R90,FV) (2473,R0,FN) (1171,R270,FN);
]

function jigsaw_test(list, jigsaw_sol, inner = false)
    tiles = get_tiles(list)
    offset = inner ? 1 : 0
    isqrt_num_tiles = isqrt(length(tiles))
    return vcat(
            [
            hcat(
                [transform(tiles[jigsaw_sol[i,j][1]], jigsaw_sol[i,j][2], jigsaw_sol[i,j][3])[1+offset:end-offset,1+offset:end-offset]
                    for j=1:isqrt_num_tiles]...
            )
            for i=1:isqrt_num_tiles]...
        )
end

aoc_image_str =
".#.#..#.##...#.##..#####
###....#.#....#..#......
##.##.###.#.#..######...
###.#####...#.#####.#..#
##.#....#.##.####...#.##
...########.#....#####.#
....#..#...##..#.#.###..
.####...#..#.....#......
#..#.##..#..###.#.##....
#.####..#.####.#.#.###..
###.#.#...#.######.#..##
#.####....##..########.#
##..##.#...#...#.#.#.#..
...#..#..#.#.##..###.###
.#.#....#.##.#...###.##.
###.#...#..#.##.######..
.#.#.###.##.##.#..#.##..
.####.###.#...###.#..#.#
..#.#..#..#.#.#.####.###
#..####...#.#.#.###.###.
#####..#####...###....##
#.##..#..#...#..####...#
.#.###..##..##..####.##.
...###...##...#...#..###"

aoc_image_sol = split(aoc_image_str, '\n')

jigsaw_sol = jigsaw_test(test_list, aoc_test_sol, true)
image_sol = [join(map(x-> x ? '#' : '.', jigsaw_sol[i,:])) for i=1:isqrt(length(jigsaw_sol))]

@show image_sol == aoc_image_sol

transformed_with_sea_monster =
".####...#####..#...###..
#####..#..#.#.####..#.#.
.#.#...#.###...#.##.O#..
#.O.##.OO#.#.OO.##.OOO##
..#O.#O#.O##O..O.#O##.##
...#.#..##.##...#..#..##
#.##.#..#.#..#..##.#.#..
.###.##.....#...###.#...
#.####.#.#....##.#..#.#.
##...#..#....#..#...####
..#.##...###..#.#####..#
....#.##.#.#####....#...
..##.##.###.....#.##..#.
#...#...###..####....##.
.#.##...#.##.#.#.###...#
#.###.#..####...##..#...
#.###...#.##...#.##O###.
.O##.#OO.###OO##..OOO##.
..O#.O..O..O.#O##O##.###
#.#..##.########..#..##.
#.#####..#.#...##..#....
#....##..#.#########..##
#...#.....#..##...###.##
#..###....##.#...##.##.#"

# count(==('#'), transformed_with_sea_monster)

nothing


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GenericBattery <: Storage
        name::String
        available::Bool
        bus::Bus
        primemover::PrimeMovers.PrimeMover
        energy::Float64
        capacity::Min_Max
        rating::Float64
        activepower::Float64
        inputactivepowerlimits::Min_Max
        outputactivepowerlimits::Min_Max
        efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
        reactivepower::Float64
        reactivepowerlimits::Union{Nothing, Min_Max}
        basepower::Float64
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data structure for a generic battery

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `primemover::PrimeMovers.PrimeMover`: PrimeMover Technology according to EIA 923
- `energy::Float64`: State of Charge of the Battery p.u.-hr, validation range: (0, nothing), action if invalid: error
- `capacity::Min_Max`: Maximum and Minimum storage capacity in p.u.-hr, validation range: (0, nothing), action if invalid: error
- `rating::Float64`
- `activepower::Float64`
- `inputactivepowerlimits::Min_Max`, validation range: (0, nothing), action if invalid: error
- `outputactivepowerlimits::Min_Max`, validation range: (0, nothing), action if invalid: error
- `efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}`, validation range: (0, 1), action if invalid: warn
- `reactivepower::Float64`, validation range: reactivepowerlimits, action if invalid: warn
- `reactivepowerlimits::Union{Nothing, Min_Max}`
- `basepower::Float64`: Base power of the unit in system base per unit, validation range: (0, nothing), action if invalid: warn
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct GenericBattery <: Storage
    name::String
    available::Bool
    bus::Bus
    "PrimeMover Technology according to EIA 923"
    primemover::PrimeMovers.PrimeMover
    "State of Charge of the Battery p.u.-hr"
    energy::Float64
    "Maximum and Minimum storage capacity in p.u.-hr"
    capacity::Min_Max
    rating::Float64
    activepower::Float64
    inputactivepowerlimits::Min_Max
    outputactivepowerlimits::Min_Max
    efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
    reactivepower::Float64
    reactivepowerlimits::Union{Nothing, Min_Max}
    "Base power of the unit in system base per unit"
    basepower::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GenericBattery(name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, basepower, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    GenericBattery(name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, basepower, services, dynamic_injector, ext, forecasts, InfrastructureSystemsInternal(), )
end

function GenericBattery(; name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, basepower, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    GenericBattery(name, available, bus, primemover, energy, capacity, rating, activepower, inputactivepowerlimits, outputactivepowerlimits, efficiency, reactivepower, reactivepowerlimits, basepower, services, dynamic_injector, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function GenericBattery(::Nothing)
    GenericBattery(;
        name="init",
        available=false,
        bus=Bus(nothing),
        primemover=PrimeMovers.BA,
        energy=0.0,
        capacity=(min=0.0, max=0.0),
        rating=0.0,
        activepower=0.0,
        inputactivepowerlimits=(min=0.0, max=0.0),
        outputactivepowerlimits=(min=0.0, max=0.0),
        efficiency=(in=0.0, out=0.0),
        reactivepower=0.0,
        reactivepowerlimits=(min=0.0, max=0.0),
        basepower=0.0,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::GenericBattery) = value.name
"""Get GenericBattery available."""
get_available(value::GenericBattery) = value.available
"""Get GenericBattery bus."""
get_bus(value::GenericBattery) = value.bus
"""Get GenericBattery primemover."""
get_primemover(value::GenericBattery) = value.primemover
"""Get GenericBattery energy."""
get_energy(value::GenericBattery) = value.energy
"""Get GenericBattery capacity."""
get_capacity(value::GenericBattery) = value.capacity
"""Get GenericBattery rating."""
get_rating(value::GenericBattery) = value.rating
"""Get GenericBattery activepower."""
get_activepower(value::GenericBattery) = value.activepower
"""Get GenericBattery inputactivepowerlimits."""
get_inputactivepowerlimits(value::GenericBattery) = value.inputactivepowerlimits
"""Get GenericBattery outputactivepowerlimits."""
get_outputactivepowerlimits(value::GenericBattery) = value.outputactivepowerlimits
"""Get GenericBattery efficiency."""
get_efficiency(value::GenericBattery) = value.efficiency
"""Get GenericBattery reactivepower."""
get_reactivepower(value::GenericBattery) = value.reactivepower
"""Get GenericBattery reactivepowerlimits."""
get_reactivepowerlimits(value::GenericBattery) = value.reactivepowerlimits
"""Get GenericBattery basepower."""
get_basepower(value::GenericBattery) = value.basepower
"""Get GenericBattery services."""
get_services(value::GenericBattery) = value.services
"""Get GenericBattery dynamic_injector."""
get_dynamic_injector(value::GenericBattery) = value.dynamic_injector
"""Get GenericBattery ext."""
get_ext(value::GenericBattery) = value.ext

InfrastructureSystems.get_forecasts(value::GenericBattery) = value.forecasts
"""Get GenericBattery internal."""
get_internal(value::GenericBattery) = value.internal

InfrastructureSystems.set_name!(value::GenericBattery, val::String) = value.name = val
"""Set GenericBattery available."""
set_available!(value::GenericBattery, val::Bool) = value.available = val
"""Set GenericBattery bus."""
set_bus!(value::GenericBattery, val::Bus) = value.bus = val
"""Set GenericBattery primemover."""
set_primemover!(value::GenericBattery, val::PrimeMovers.PrimeMover) = value.primemover = val
"""Set GenericBattery energy."""
set_energy!(value::GenericBattery, val::Float64) = value.energy = val
"""Set GenericBattery capacity."""
set_capacity!(value::GenericBattery, val::Min_Max) = value.capacity = val
"""Set GenericBattery rating."""
set_rating!(value::GenericBattery, val::Float64) = value.rating = val
"""Set GenericBattery activepower."""
set_activepower!(value::GenericBattery, val::Float64) = value.activepower = val
"""Set GenericBattery inputactivepowerlimits."""
set_inputactivepowerlimits!(value::GenericBattery, val::Min_Max) = value.inputactivepowerlimits = val
"""Set GenericBattery outputactivepowerlimits."""
set_outputactivepowerlimits!(value::GenericBattery, val::Min_Max) = value.outputactivepowerlimits = val
"""Set GenericBattery efficiency."""
set_efficiency!(value::GenericBattery, val::NamedTuple{(:in, :out), Tuple{Float64, Float64}}) = value.efficiency = val
"""Set GenericBattery reactivepower."""
set_reactivepower!(value::GenericBattery, val::Float64) = value.reactivepower = val
"""Set GenericBattery reactivepowerlimits."""
set_reactivepowerlimits!(value::GenericBattery, val::Union{Nothing, Min_Max}) = value.reactivepowerlimits = val
"""Set GenericBattery basepower."""
set_basepower!(value::GenericBattery, val::Float64) = value.basepower = val
"""Set GenericBattery services."""
set_services!(value::GenericBattery, val::Vector{Service}) = value.services = val
"""Set GenericBattery ext."""
set_ext!(value::GenericBattery, val::Dict{String, Any}) = value.ext = val

InfrastructureSystems.set_forecasts!(value::GenericBattery, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set GenericBattery internal."""
set_internal!(value::GenericBattery, val::InfrastructureSystemsInternal) = value.internal = val

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light,md:markdown
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.0
#     language: julia
#     name: julia-1.0
# ---

using DataFrames
using GraphDataFrameBridge
using MetaGraphs
using CSV
using LightGraphs
using GraphPlot
#using SemanticModels

# ### Objective 1:
#
# Build a toy knowledge graph that can represent the component pieces of our flu example, and (when traversed) can output a sequence of function calls capable of reproducing our manually produced metamodel.

# In metagraphs you can make a DataFrame with columns you want like ```src,dst,eprop1,eprop2,eprop3```... with an arbitrary number of columns and then ```g = MetaGraph(dataframe)``` will populate all the edges with those properties. Then we could add the vertex properties with a loop to set those vertex properties from a second dataframe ```vertex, vprop1, vprop2, vprop3```...
#
# The properties are stored in a dictionary so we can change them at run time without making new types. The knowledge graph will have

# ##### Vertex Types:
#
# - Code symbols
# - Functions
# - Variables
# - Symbols
#
# ##### Math/Science Concepts:
# - Math Expressions
# - Concept
# - Unit

# ##### Edge Types: SourceType -> DestinationType
# - IsCalledBy: function -> function
# - Co-occurs with: Any -> Any
# - IsComponentOf: {Symbol,Function,Variable} -> {Function, Expression}
# - IsMeasuredIn: Any -> Unit
# - Implements: Function -> {Math Expression, Concept}
# - VERB: Any -> Any
# - IsSubClassOf: SpecificConcept -> More General Concept

'''
function create_sir(m, solns)
        sol = solns[1]
        initialS = u"10000person"
        initialI = u"1person"
        initialpop = [initialS, initialI, u"0.0person"]
        β = u"1.0/18"/u"d*C" * sol(sol.t[end-2])[1] #infectiousness
        @show β
        sirprob = SIRSimulation(initialpop, #initial_conditions S,I,R
                                (u"0.0d", u"20d"), #time domain
                                SIRParams(β, u"40.0person/d")) # parameters β, γ
        return sirprob
    end
'''

function create_bidirectional_edges(edges_df)
    for edge in eachrow(edges_df)
        if edge.CoOccursWith == 1
            push!(edges_df, [edge.dst_vid,
                    edge.dst_name,
                    edge.src_vid,
                    edge.src_name,
                    edge.IsCalledBy,
                    edge.CoOccursWith,
                    edge.IsComponentOf,
                    edge.IsMeasuredIn,
                    edge.Implements,
                    edge.IsSubClassOf,
                    edge.Verb,
                    edge.VerbToken])
        end
    end
    return edges_df
end

function generate_kg(vertices, edges)
    n_vertices = size(vertices,1)
    g = MetaDiGraph(n_vertices)

    for vertex in eachrow(vertices)
        set_props!(g, vertex.v_id,
            Dict(:v_name=>vertex.v_name, :v_type=>vertex.v_type, :v_math_sci_concept=> vertex.v_math_sci_concept))
    end

    for edge in eachrow(edges)
        add_edge!(g, edge.src_vid, edge.dst_vid)

        edge_attrs = Dict(
            :IsCalledBy=>edge.IsCalledBy,
            :CoOccursWith=>edge.CoOccursWith,
            :IsComponentOf=>edge.IsComponentOf,
            :IsMeasuredIn=>edge.IsMeasuredIn,
            :Implements=>edge.Implements,
            :IsSubClassOf=>edge.IsSubClassOf,
            :Verb=>edge.Verb,
            :VerbToken=>edge.VerbToken)

        set_props!(g, Edge(edge.src_vid, edge.dst_vid), edge_attrs)

    end

    nodelabel = [vertices.v_name]

    return g
end

schema = CSV.read("./data/kg_schema.csv")

vertices = CSV.read("./data/kg_vertices.csv")

edges = CSV.read("./data/kg_edges.csv")
edges = edges[1:5, [:src_vid, :src_name, :dst_vid, :dst_name, :IsCalledBy, :CoOccursWith, :IsComponentOf, :IsMeasuredIn, :Implements, :IsSubClassOf, :Verb, :VerbToken]]

edges_full = create_bidirectional_edges(edges)

g = generate_kg(vertices, edges_full)
display(g)

gplot(g, nodelabel=vertices.v_name, layout=circular_layout)


text:
"""
Subtypes are different approaches for solving an optimization problem. We use
"approach" to mean a configured algorithm: that is, the algorithm (e.g.
simplex), along with fixed values for all the algorithm hyperparemeters.
"""
abstract type AbstractApproach end

"""
Subtypes contain all the information for an individual instance of an
optimization problem. They either contain all of the requisite data themselves,
or a mechanism to somehow retreive it from, e.g., a website.
"""
abstract type AbstractInstance end

"""
Subtypes encapsulate an single approach applied to a single instance. Callers
can ask to optimize a model and then query the results, but the workings of the
approach are otherwise treated as a black box.
"""
abstract type AbstractModel end

"""
Subtypes contain all the useful information needed for understanding how a
particular approach performed on a particular instance. Typical fields stored
may include solve time, or node count for a MIP-based approach.
"""
abstract type AbstractResult end

"""
Subtypes store the results for experiments across multiple approaches and
multiple instances. Backends could be as simple as a CSV file, or something
more complex like HDF5.
"""
abstract type AbstractResultTable end

"""
    prep_instance!(instance::AbstractInstance)::Nothing

Prepare an instance for experimentation. This is useful if there is a
lightweight way to represent an instance (e.g. a string name in an instance
library), but collecting the data is relatively costly (e.g. download something
from the web and converting it to the right file format). Additionally, this
prep work can do more complex pre-processing: if you are evaluating methods for
dual bounds, this function could run a set of heuristics to produce a good
primal solution. This method is only called once per instance, meaning that the
resulting computations can be shared across all approaches run on the
particular instance.

!!! note
    Implementing this method is optional. If you don't need to do any prep work
    as all the requisite data is already available, simply defer to the
    available fallback method.
"""
function prep_instance! end
prep_instance!(::AbstractInstance) = nothing

"""
    build_model(
        approach::AbstractApproach,
        instance::AbstractInstance,
        config::Dict{String,Any}=Dict{String,Any}()
    )::AbstractModel

Given a particular `approach` and a particular `instance`, build a model that
encapsulates the two. The optional `config` argument contains parameters that
are shared across all approaches and instances in a experiment run: for
example, a time limit, a node limit, a constraint violation tolerance, etc.
"""
function build_model end

"""
    optimize!(model::AbstractModel)

Solve the instance using the approach specified by `model`.
"""
function optimize! end

"""
    tear_down(model::AbstractModel, instance::AbstractInstance, T::Type{<:AbstractResult}))

Record the results of in individual experiment of a particular approach on a
particular instance. The results will be stored in an object of type `T`.

!!! note
    This function should only be called after `optimize!(model)` has been
    called.
"""
function tear_down end

"""
    record_result!(table::AbstractResultTable, result::AbstractResult)

Given the results from an individual experiment, record it in `table`.
"""
function record_result! end


text:
using AgentBasedModels
using Base.Test

# Mock simulation model
type Foo
    bar::Int
end
init(; start=1, kwargs...) = Foo(start)
step!(s::Foo; kwargs...) = (s.bar += 2; nothing)

model = AgentBasedModel(init, step!)
AgentBasedModels.init!(model)
AgentBasedModels.step!(model)
@test isa(state(model), Foo)
@test state(model).bar == 3

init_kwargs = [:start=>2]
simulate(model, steps=200, init_kwargs = init_kwargs)
@test state(model).bar == 402

init_output(f::Foo) = [:start => f.bar, :bar => Int[f.bar]]
output!(r::Result{Foo}) = push!(r.output[:bar], r.abm.state.bar)
res = simulate(model; steps=10, init_kwargs=init_kwargs, init_output=init_output, output! = output!)
@test res.output[:bar] == [2:2:22;]


text:
mutable struct pointcloud_t <: LCMType
    utime::Int64
    seq::Int32
    frame_id::String
    n_points::Int32
    points::Matrix{Float32}

    n_channels::Int32
    channel_names::Vector{String}
    channels::Matrix{Float32}
end

@lcmtypesetup(pointcloud_t,
    points => (n_points, 3),
    channel_names => (n_channels,),
    channels => (n_channels, n_points),
)


text:
# Given a project on disk, loop through mutations and runs of unit tests.

using Pkg
using TestReports
using SourceWalk
using Vimes

"""Replace a mutated file with an original file"""
function reset(dir, tmp, f)
    relative_f = f[length(tmp) + 2:end]
    rm(f)
    cp(joinpath(dir, relative_f), f)
    return
end

function mutate_and_reset(callback::Function, dir, tmp, idx)
    f = nothing
    while isempty(SourceWalk.diff(dir, tmp))
      f = Vimes.mutate(tmp, idx)
    end
    callback()
    reset(dir, tmp, f)
end

function initialise_noclean(dir)
    (isfile(joinpath(dir, "Project.toml")) && isdir(joinpath(dir, "src"))) ||
      error("No Julia project found at $dir")
    tmp = joinpath(tempdir(), "vimes-$(rand(UInt64))")
    mkdir(tmp)
    for path in readdir(dir)
        if !startswith(path, ".")
            cp(joinpath(dir, path), joinpath(tmp, path))
        end
    end
    return tmp
end

function setup_working(project_dir)
    pkgname = splitext(basename(project_dir))[1]
    copy_dir = initialise_noclean(expanduser(project_dir))
    # TestReports expects the package to be installed in order to test it.
    if !isnothing(Base.find_package(pkgname))
        Pkg.rm(pkgname)
    end
    Pkg.develop(path = copy_dir)
    (copy_dir, pkgname)
end

function teardown_working(copy_dir, pkgname)
    Pkg.rm(pkgname)
    rm(copy_dir, recursive=true)
end

"""
For some reason, this creates a test log only once, the first time it is called.
It makes that log in the directory of the project itself. It makes another, empty
log in the current directory. The empty log always says there were no failures.
I can't figure out why it won't write another log.
"""
function run_testreports(pkgname)
    try
        run(`$(Base.julia_cmd()) --project=. -e "using TestReports; TestReports.test($(repr(pkgname)))"`)
    catch e
        if e isa ProcessFailedException
            false
        else
            throw(e)
        end
    end
    true
end

next_log(file_list) =
    maximum(map(x -> parse(Int, match(r"[0-9]+", x).match),
                filter(x->endswith(x, ".xml"), file_list))) + 1

function generate_mutation_reports(project_dir, tmp, cnt::UnitRange, pkgname, logdir = pwd())
    mutant_index = Vimes.indices(joinpath(tmp, "src"), Vimes.defaults)
    begin_idx = next_log(readdir(logdir))
    for mutation_idx in cnt .+ begin_idx
        `$(Base.julia_cmd()) --project=$tmp -e 'using Pkg; Pkg.test()'`
        Vimes.mutate_and_reset(project_dir, tmp, mutant_index) do
            testlog = joinpath(tmp, joinpath(logdir, "log$(lpad(mutation_idx, 5, "0")).xml"))
            runonce(tmp, testlog, pkgname)
        end
    end
end

function generate_reports(project_dir, cnt)
    project_dir = expanduser(project_dir)
    isdir(project_dir) || error("Cannot find directory $(project_dir)")
    project_copy, pkgname = setup_working(project_dir)
    println("cnt is $cnt")
    generate_mutation_reports(project_dir, project_copy, 1:cnt, pkgname)
    teardown_working(project_copy, pkgname)
end

"""
This is a small equivalent to TestReports.test(["BijectiveHilbert"]).
It's here because that function seems to write only one report
and won't write a second one.
"""
function runonce(tmp, logfile, pkgname)
    runtests = joinpath(tmp, "test", "runtests.jl")
    runner_code = """
using Test
using TestReports
using $(pkgname)

append!(empty!(ARGS), String[])

ts = @testset ReportingTestSet "" begin
    include($(repr(runtests)))
end

write($(repr(logfile)), report(ts))
"""
    println(runner_code)
    test_process = open(`$(Base.julia_cmd()) --project=. -e $(runner_code)`, read = false)
    # Let the existence of the XML report tell us success or failure.
    wait(test_process)
end


text:
using Oceananigans, Printf, PyPlot

include("utils.jl")

function informative_message(model, u, w, ℕ)
    return @sprintf("""
        This is an informative message.

         model vertical resolution : %d
                   model iteration : %d
        kinetic + potential energy : %.2e
                    kinetic energy : %.2e
               relative error in w : %.2e
               relative error in u : %.2e
        """, model.grid.Nz, model.clock.iteration, total_energy(model, ℕ),
        total_kinetic_energy(model), w_relative_error(model, w), u_relative_error(model, u))
end

function makeplot(axs, model, w)
    w_num = model.velocities.w

    w_ans = FaceFieldZ(w.(xnodes(w_num), ynodes(w_num), znodes(w_num),
                        model.clock.time), model.grid)

    wmax = maximum(abs.(w_num.data))

    sca(axs[1, 1])

    PyPlot.plot(data(w_ans)[1, 1, :], view(znodes(w_num), 1, 1, :),  "-", linewidth=2, alpha=0.4)
    PyPlot.plot(data(w_num)[1, 1, :], view(znodes(w_num), 1, 1, :), "--", linewidth=1)

    xlim(-wmax, wmax)

    axs[1, 1].spines["top"].set_visible(false)
    axs[1, 1].spines["right"].set_visible(false)

    xlabel(L"w")
    ylabel(L"z")

    sca(axs[2, 1])
    PyPlot.plot(data(w_ans)[1, 1, :] .- data(w_num)[1, 1, :], view(znodes(w_num), 1, 1, :), "-")

    xlim(-wmax, wmax)

    axs[2, 1].spines["top"].set_visible(false)
    axs[2, 1].spines["right"].set_visible(false)

    xlabel(L"\mathrm{mean} \left [ \left ( w_\mathrm{model} - w_\mathrm{analytical} \right )^2 \right ] / \mathrm{mean} \left ( w_\mathrm{analytical}^2 \right )")
    ylabel(L"z")

    sca(axs[1, 2])
    plotxzslice(w_num)
    axis("off")
    title("Model vertical velocity")
    xlabel(L"x")
    ylabel(L"z")

    sca(axs[2, 2])
    plotxzslice(w_ans)
    axis("off")
    title("Analytical vertical velocity")
    xlabel(L"x")
    ylabel(L"z")

    return nothing
end

# Internal wave parameters
 ν = κ = 1e-9
 L = 2π
z₀ = -L/3
 δ = L/20
a₀ = 1e-3
 m = 16
 k = 1
 f = 0.2
 ℕ = 1.0
 σ = sqrt( (ℕ^2*k^2 + f^2*m^2) / (k^2 + m^2) )

# Numerical parameters
 N = 128
Δt = 0.01 * 1/σ

@show σ/f
cᵍ = m * σ / (k^2 + m^2) * (f^2/σ^2 - 1)
 U = a₀ * k * σ   / (σ^2 - f^2)
 V = a₀ * k * f   / (σ^2 - f^2)
 W = a₀ * m * σ   / (σ^2 - ℕ^2)
 Θ = a₀ * m * ℕ^2 / (σ^2 - ℕ^2)

a(x, y, z, t) = exp( -(z - cᵍ*t - z₀)^2 / (2*δ)^2 )

u(x, y, z, t) =           a(x, y, z, t) * U * cos(k*x + m*z - σ*t)
v(x, y, z, t) =           a(x, y, z, t) * V * sin(k*x + m*z - σ*t)
w(x, y, z, t) =           a(x, y, z, t) * W * cos(k*x + m*z - σ*t)
T(x, y, z, t) = ℕ^2 * z + a(x, y, z, t) * Θ * sin(k*x + m*z - σ*t)

u₀(x, y, z) = u(x, y, z, 0)
v₀(x, y, z) = v(x, y, z, 0)
w₀(x, y, z) = w(x, y, z, 0)
T₀(x, y, z) = T(x, y, z, 0)

# Create a model where temperature = buoyancy.
model = BasicModel(N=(N, 1, N), L=(L, L, L), ν=ν, κ=κ,
                    eos=LinearEquationOfState(βT=1.0),
                    constants=PlanetaryConstants(f=f, g=1.0))

set_ic!(model, u=u₀, v=v₀, w=w₀, T=T₀)
println(informative_message(model, u, w, ℕ))

fig, axs = subplots(ncols=2, nrows=2, figsize=(8, 8))

for Nt in (1, 10, 100)
    time_step!(model, Nt, Δt)
    makeplot(axs, model, w)
    println(informative_message(model, u, w, ℕ))
end

gcf()


text:
##############################################################################################
module runtests
##############################################################################################
##############################################################################################
# imports and exports
using Test
using JLD2

##############################################################################################
##############################################################################################
# include test sub-modules
include("TestDeSolvers.jl")
include("TestL96.jl")
include("TestTimeSeriesGeneration.jl")
include("TestIEEE39bus.jl")
include("TestFilterExps.jl")
include("TestClassicSmootherExps.jl")
include("TestIterativeSmootherExps.jl")
include("TestSingleIterationSmootherExps.jl")

##############################################################################################
##############################################################################################
# Run tests

# test set 1: Calculate the order of convergence for standard integrators
@testset "Calculate Order Convergence" begin
    @test TestDeSolvers.testEMExponential()
    @test TestDeSolvers.testRKExponential()
end

# test set 2: test L96 model equations for known behavior
@testset "Lorenz-96" begin
    @test TestL96.EMZerosStep()
    @test TestL96.EMFStep()
end

# test set 3: Test time series generation, saving output to default directory and loading
@testset "Time Series Generation" begin
    @test TestTimeSeriesGeneration.testGenL96()
    @test TestTimeSeriesGeneration.testLoadL96()
    @test TestTimeSeriesGeneration.testGenIEEE39bus()
    @test TestTimeSeriesGeneration.testLoadIEEE39bus()
end

# test set 4: test the model equations for known behavior
@testset "IEEE 39 Bus" begin
    @test TestIEEE39bus.test_synchrony()
end

# test set 5: test filter state and parameter experiments
@testset "Filter Experiments" begin
    @test TestFilterExps.run_filter_state_L96()
    @test TestFilterExps.analyze_filter_state_L96()
    @test TestFilterExps.run_filter_param_L96()
    @test TestFilterExps.analyze_filter_param_L96()
    @test TestFilterExps.run_filter_state_IEEE39bus()
    @test TestFilterExps.analyze_filter_state_IEEE39bus()
end

# test set 6: test classic smoother state and parameter experiments
@testset "Classic Smoother Experiments" begin
    @test TestClassicSmootherExps.run_smoother_state_L96()
    @test TestClassicSmootherExps.analyze_smoother_state_L96()
    @test TestClassicSmootherExps.run_smoother_param_L96()
    @test TestClassicSmootherExps.analyze_smoother_param_L96()
end

# test set 7: test IEnKS smoother state and parameter experiments
@testset "Iterative Smoother Experiments" begin
    @test TestIterativeSmootherExps.run_sda_smoother_state_L96()
    @test TestIterativeSmootherExps.analyze_sda_smoother_state_L96()
    @test TestIterativeSmootherExps.run_sda_smoother_param_L96()
    @test TestIterativeSmootherExps.analyze_sda_smoother_param_L96()
    @test TestIterativeSmootherExps.run_sda_smoother_state_L96()
    @test TestIterativeSmootherExps.analyze_sda_smoother_state_L96()
    @test TestIterativeSmootherExps.run_sda_smoother_param_L96()
    @test TestIterativeSmootherExps.analyze_sda_smoother_param_L96()
end

# test set 8: test SIEnKS smoother state and parameter experiments
@testset "Single Iteration Smoother Experiments" begin
    @test TestSingleIterationSmootherExps.run_sda_smoother_state_L96()
    @test TestSingleIterationSmootherExps.analyze_sda_smoother_state_L96()
    @test TestSingleIterationSmootherExps.run_sda_smoother_param_L96()
    @test TestSingleIterationSmootherExps.analyze_sda_smoother_param_L96()
    @test TestSingleIterationSmootherExps.run_mda_smoother_state_L96()
    @test TestSingleIterationSmootherExps.analyze_mda_smoother_state_L96()
    @test TestSingleIterationSmootherExps.run_mda_smoother_param_L96()
    @test TestSingleIterationSmootherExps.analyze_mda_smoother_param_L96()
end

##############################################################################################
# end module

end

max_stars_repo_name: cgrudz/DataAssimilationBenchmarks

max_issues_repo_name: cgrudz/DataAssimilationBenchmarks


text:
using StaticArrays
using BenchmarkTools

const  ⊗ = kron

function cost(𝛉::AbstractArray, 𝒞::Tuple{AbstractArray, Vararg{AbstractArray}}, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ = collect(𝒟)
    Λ₁, Λ₂ = collect(𝒞)
    Jₐₘₗ = 0.0
    N = length(𝒟[1])
    𝚲ₙ = @MMatrix zeros(4,4)
    # 𝐞₁ = @SVector [1.0, 0.0, 0.0]
    # 𝐞₂ = @SVector [0.0, 1.0, 0.0]
    𝐞₁ = @SMatrix [1.0; 0.0; 0.0]
    𝐞₂ = @SMatrix [0.0; 1.0; 0.0]
    @inbounds for n = 1:N
        index = SVector(1, 2)
        𝚲ₙ[1:2,1:2] .= Λ₁[n][index, index]
        𝚲ₙ[3:4,3:4] .= Λ₂[n][index, index]
        𝐦 = ℳ[n]
        𝐦ʹ= ℳʹ[n]
        𝐔ₙ = (𝐦 ⊗ 𝐦ʹ)
        ∂ₓ𝐮ₙ =  [(𝐞₁ ⊗ 𝐦ʹ) (𝐞₂ ⊗ 𝐦ʹ) (𝐦 ⊗ 𝐞₁) (𝐦 ⊗ 𝐞₂)]
        𝐁ₙ =  ∂ₓ𝐮ₙ * 𝚲ₙ * ∂ₓ𝐮ₙ'
        𝚺ₙ = 𝛉' * 𝐁ₙ * 𝛉
        Jₐₘₗ += 𝛉' * 𝐔ₙ * (𝚺ₙ \ 𝐔ₙ') * 𝛉
    end
    Jₐₘₗ
end

# Some sample data
N = 3376821
ℳ = [@MMatrix(rand(3,1)) for i = 1:N]
ℳʹ = [@MMatrix(rand(3,1)) for i = 1:N]
Λ₁ =  [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(ℳ)]
Λ₂ =  [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(ℳ)]
t = @MVector rand(9)
𝒞 = (Λ₁,Λ₂)
𝒟 = (ℳ, ℳʹ)

cost(t,𝒞 , 𝒟)
@time cost(t,𝒞 , 𝒟)
@btime cost($t,$𝒞 , $𝒟)


text:
module GibbsSeaWater
import Base.beta

export add_barrier, add_mean, adiabatic_lapse_rate_from_ct, alpha,
       alpha_on_beta, alpha_wrt_t_exact, beta_const_t_exact, beta, c_from_sp,
       cabbeling, ct_freezing, ct_from_pt, ct_from_t, deltasa_atlas,
       deltasa_from_sp, dynamic_enthalpy, enthalpy, enthalpy_sso_0_p,
       enthalpy_t_exact, cp_t_exact, entropy_from_t, entropy_part,
       entropy_part_zerop, fdelta, gibbs, gibbs_pt0_pt0, grav,
       hill_ratio_at_sp2, indx, internal_energy, ipv_vs_fnsquared_ratio, kappa,
       kappa_t_exact, latentheat_evap_ct, latentheat_evap_t,
       latentheat_melting, nsquared, pot_rho_t_exact, pt0_from_t, pt_from_ct,
       pt_from_t, rho, rho_first_derivatives, rho_t_exact, saar, sa_from_rho,
       sa_from_sp_baltic, sa_from_sp, sa_from_sstar, sigma0, sigma1, sigma2,
       sigma3, sigma4, sound_speed, sound_speed_t_exact, specvol_anom, specvol,
       specvol_sso_0_p, specvol_t_exact, sp_from_c, sp_from_sa_baltic,
       sp_from_sa, sp_from_sk, sp_from_sr, sp_from_sstar, sr_from_sp,
       sstar_from_sa, sstar_from_sp, t_freezing, t_from_ct, thermobaric,
       turner_rsubrho, xinterp1, z_from_p

ext = @linux? (".so" : @osx? ( ".dylib" : @windows? ( ".dll" : "" )))
if ext == ""
    error("Platform not linux, OS X, or Windows")
end

const _MODULEPATH = Pkg.dir("GibbsSeaWater")
const _LIBGSWPATH = joinpath(_MODULEPATH, "deps/build/libgswteos-10"*ext)

function add_barrier(input_data, lon, lat, long_grid, lat_grid, dlong_grid, dlat_grid, output_data)
    return ccall((:gsw_add_barrier, _LIBGSWPATH), Void, (Ptr{Float64}, Float64, Float64, Float64, Float64, Float64, Float64, Ptr{Float64}), pointer(input_data), lon, lat, long_grid, lat_grid, dlong_grid, dlat_grid, pointer(output_data))
end

function add_mean(data_in, lon, lat, data_out)
    return ccall((:gsw_add_mean, _LIBGSWPATH), Void, (Ptr{Float64}, Float64, Float64, Ptr{Float64}), pointer(data_in), lon, lat, pointer(data_out))
end

function adiabatic_lapse_rate_from_ct(sa, ct, p)
    return ccall((:gsw_adiabatic_lapse_rate_from_ct, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function alpha(sa, ct, p)
    return ccall((:gsw_alpha, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function alpha_on_beta(sa, ct, p)
    return ccall((:gsw_alpha_on_beta, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function alpha_wrt_t_exact(sa, t, p)
    return ccall((:gsw_alpha_wrt_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function beta_const_t_exact(sa, t, p)
    return ccall((:gsw_beta_const_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function beta(sa, ct, p)
    return ccall((:gsw_beta, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function c_from_sp(sp, t, p)
    return ccall((:gsw_c_from_sp, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sp, t, p)
end

function cabbeling(sa, ct, p)
    return ccall((:gsw_cabbeling, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function ct_freezing(sa, p, saturation_fraction)
    return ccall((:gsw_ct_freezing, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, p, saturation_fraction)
end

function ct_from_pt(sa, pt)
    return ccall((:gsw_ct_from_pt, _LIBGSWPATH), Float64, (Float64, Float64), sa, pt)
end

function ct_from_t(sa, t, p)
    return ccall((:gsw_ct_from_t, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function deltasa_atlas(p, lon, lat)
    return ccall((:gsw_deltasa_atlas, _LIBGSWPATH), Float64, (Float64, Float64, Float64), p, lon, lat)
end

function deltasa_from_sp(sp, p, lon, lat)
    return ccall((:gsw_deltasa_from_sp, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sp, p, lon, lat)
end

function dynamic_enthalpy(sa, ct, p)
    return ccall((:gsw_dynamic_enthalpy, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function enthalpy(sa, ct, p)
    return ccall((:gsw_enthalpy, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function enthalpy_sso_0_p(p)
    return ccall((:gsw_enthalpy_sso_0_p, _LIBGSWPATH), Float64, (Float64,), p)
end

function enthalpy_t_exact(sa, t, p)
    return ccall((:gsw_enthalpy_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function cp_t_exact(sa, t, p)
    return ccall((:gsw_cp_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function entropy_from_t(sa, t, p)
    return ccall((:gsw_entropy_from_t, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function entropy_part(sa, t, p)
    return ccall((:gsw_entropy_part, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function entropy_part_zerop(sa, pt0)
    return ccall((:gsw_entropy_part_zerop, _LIBGSWPATH), Float64, (Float64, Float64), sa, pt0)
end

function fdelta(p, lon, lat)
    return ccall((:gsw_fdelta, _LIBGSWPATH), Float64, (Float64, Float64, Float64), p, lon, lat)
end

function gibbs(ns, nt, np, sa, t, p)
    return ccall((:gsw_gibbs, _LIBGSWPATH), Float64, (Int64, Int64, Int64, Float64, Float64, Float64), ns, nt, np, sa, t, p)
end

function gibbs_pt0_pt0(sa, pt0)
    return ccall((:gsw_gibbs_pt0_pt0, _LIBGSWPATH), Float64, (Float64, Float64), sa, pt0)
end

function grav(lat, p)
    return ccall((:gsw_grav, _LIBGSWPATH), Float64, (Float64, Float64), lat, p)
end

function hill_ratio_at_sp2(t)
    return ccall((:gsw_hill_ratio_at_sp2, _LIBGSWPATH), Float64, (Float64,), t)
end

function indx(x, n, z)
    return ccall((:gsw_indx, _LIBGSWPATH), Int64, (Ptr{Float64}, Int64, Float64), pointer(x), n, z)
end

function internal_energy(sa, ct, p)
    return ccall((:gsw_internal_energy, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function ipv_vs_fnsquared_ratio(sa, ct, p, nz, ipv_vs_fnsquared_ratio, p_mid)
    return ccall((:gsw_ipv_vs_fnsquared_ratio, _LIBGSWPATH), Void, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Float64}), pointer(sa), pointer(ct), pointer(p), nz, pointer(ipv_vs_fnsquared_ratio), pointer(p_mid))
end

function kappa(sa, ct, p)
    return ccall((:gsw_kappa, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function kappa_t_exact(sa, t, p)
    return ccall((:gsw_kappa_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function latentheat_evap_ct(sa, ct)
    return ccall((:gsw_latentheat_evap_ct, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function latentheat_evap_t(sa, t)
    return ccall((:gsw_latentheat_evap_t, _LIBGSWPATH), Float64, (Float64, Float64), sa, t)
end

function latentheat_melting(sa, p)
    return ccall((:gsw_latentheat_melting, _LIBGSWPATH), Float64, (Float64, Float64), sa, p)
end

function nsquared(sa, ct, p, lat, nz, n2, p_mid)
    return ccall((:gsw_nsquared, _LIBGSWPATH), Void, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Float64}), pointer(sa), pointer(ct), pointer(p), pointer(lat), nz, pointer(n2), pointer(p_mid))
end

function pot_rho_t_exact(sa, t, p, p_ref)
    return ccall((:gsw_pot_rho_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sa, t, p, p_ref)
end

function pt0_from_t(sa, t, p)
    return ccall((:gsw_pt0_from_t, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function pt_from_ct(sa, ct)
    return ccall((:gsw_pt_from_ct, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function pt_from_t(sa, t, p, p_ref)
    return ccall((:gsw_pt_from_t, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sa, t, p, p_ref)
end

function rho(sa, ct, p)
    return ccall((:gsw_rho, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function rho_first_derivatives(sa, ct, p, drho_dsa, drho_dct, drho_dp)
    return ccall((:gsw_rho_first_derivatives, _LIBGSWPATH), Void, (Float64, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), sa, ct, p, pointer(drho_dsa), pointer(drho_dct), pointer(drho_dp))
end

function rho_t_exact(sa, t, p)
    return ccall((:gsw_rho_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function saar(p, lon, lat)
    return ccall((:gsw_saar, _LIBGSWPATH), Float64, (Float64, Float64, Float64), p, lon, lat)
end

function sa_from_rho(rho, ct, p)
    return ccall((:gsw_sa_from_rho, _LIBGSWPATH), Float64, (Float64, Float64, Float64), rho, ct, p)
end

function sa_from_sp_baltic(sp, lon, lat)
    return ccall((:gsw_sa_from_sp_baltic, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sp, lon, lat)
end

function sa_from_sp(sp, p, lon, lat)
    return ccall((:gsw_sa_from_sp, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sp, p, lon, lat)
end

function sa_from_sstar(sstar, p, lon, lat)
    return ccall((:gsw_sa_from_sstar, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sstar, p, lon, lat)
end

function sigma0(sa, ct)
    return ccall((:gsw_sigma0, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function sigma1(sa, ct)
    return ccall((:gsw_sigma1, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function sigma2(sa, ct)
    return ccall((:gsw_sigma2, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function sigma3(sa, ct)
    return ccall((:gsw_sigma3, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function sigma4(sa, ct)
    return ccall((:gsw_sigma4, _LIBGSWPATH), Float64, (Float64, Float64), sa, ct)
end

function sound_speed(sa, ct, p)
    return ccall((:gsw_sound_speed, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function sound_speed_t_exact(sa, t, p)
    return ccall((:gsw_sound_speed_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function specvol_anom(sa, ct, p)
    return ccall((:gsw_specvol_anom, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function specvol(sa, ct, p)
    return ccall((:gsw_specvol, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function specvol_sso_0_p(p)
    return ccall((:gsw_specvol_sso_0_p, _LIBGSWPATH), Float64, (Float64,), p)
end

function specvol_t_exact(sa, t, p)
    return ccall((:gsw_specvol_t_exact, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, t, p)
end

function sp_from_c(c, t, p)
    return ccall((:gsw_sp_from_c, _LIBGSWPATH), Float64, (Float64, Float64, Float64), c, t, p)
end

function sp_from_sa_baltic(sa, lon, lat)
    return ccall((:gsw_sp_from_sa_baltic, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, lon, lat)
end

function sp_from_sa(sa, p, lon, lat)
    return ccall((:gsw_sp_from_sa, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sa, p, lon, lat)
end

function sp_from_sk(sk)
    return ccall((:gsw_sp_from_sk, _LIBGSWPATH), Float64, (Float64,), sk)
end

function sp_from_sr(sr)
    return ccall((:gsw_sp_from_sr, _LIBGSWPATH), Float64, (Float64,), sr)
end

function sp_from_sstar(sstar, p, lon, lat)
    return ccall((:gsw_sp_from_sstar, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sstar, p, lon, lat)
end

function sr_from_sp(sp)
    return ccall((:gsw_sr_from_sp, _LIBGSWPATH), Float64, (Float64,), sp)
end

function sstar_from_sa(sa, p, lon, lat)
    return ccall((:gsw_sstar_from_sa, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sa, p, lon, lat)
end

function sstar_from_sp(sp, p, lon, lat)
    return ccall((:gsw_sstar_from_sp, _LIBGSWPATH), Float64, (Float64, Float64, Float64, Float64), sp, p, lon, lat)
end

function t_freezing(sa, p, saturation_fraction)
    return ccall((:gsw_t_freezing, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, p, saturation_fraction)
end

function t_from_ct(sa, ct, p)
    return ccall((:gsw_t_from_ct, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function thermobaric(sa, ct, p)
    return ccall((:gsw_thermobaric, _LIBGSWPATH), Float64, (Float64, Float64, Float64), sa, ct, p)
end

function turner_rsubrho(sa, ct, p, nz, tu, rsubrho, p_mid)
    return ccall((:gsw_turner_rsubrho, _LIBGSWPATH), Void, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), pointer(sa), pointer(ct), pointer(p), nz, pointer(tu), pointer(rsubrho), pointer(p_mid))
end

function xinterp1(x, y, n, x0)
    return ccall((:gsw_xinterp1, _LIBGSWPATH), Float64, (Ptr{Float64}, Ptr{Float64}, Int64, Float64), pointer(x), pointer(y), n, x0)
end

function z_from_p(p, lat)
    return ccall((:gsw_z_from_p, _LIBGSWPATH), Float64, (Float64, Float64), p, lat)
end

end # module


text:
using JLD2

@testset "Test issue #83" begin
    nxrec = 120
    nyrec = 100
    xrec = range(50f0, stop=1150f0, length=nxrec)
    yrec = range(100f0, stop=900f0, length=nyrec)
    zrec = 50f0

    # Construct 3D grid from basis vectors
    (x3d, y3d, z3d) = setup_3D_grid(xrec, yrec, zrec)

    for k=1:nyrec
        s = (k - 1) * nxrec + 1
        e = s + nxrec - 1
        @test all(x3d[k:nxrec:end] .== xrec[k])
        @test all(y3d[s:e] .== yrec[k])
    end
    @test all(z3d .== zrec)
end

@testset "Test backward comp issue" begin
    datapath = joinpath(dirname(pathof(JUDI)))*"/../data/"

    # Load file with old judiVector type and julia <1.7 StepRangeLen
    @load "$(datapath)backward_comp.jld" dat model_true

    @test typeof(dat) == judiVector{Float32, Matrix{Float32}}
    @test typeof(dat.geometry) == GeometryIC{Float32}
    @test typeof(dat.geometry.xloc) == Vector{Vector{Float32}}
    @test typeof(model_true) == Model
    @test typeof(model_true.m) == PhysicalParameter{Float32}
end


text:
__precompile__()
module st16SseqJuliaTools
    using Statistics
    using StatsBase
    using FastaIO
    using AutoHashEquals
    using DataStructures
    using DataFrames
    using Conda
    using BioSequences
    using FASTX
    using JLD2
    using Distributions
    using BioTools.BLAST
    using IntervalSets
    using ProgressMeter
    using Artifacts
    using BioAlignments

    #Alignment analysis stuff
    export Alignment
    export read_fasta
    export write_to_fasta
    export sub_alignment
    export vsearh_cluster!
    export get_max_cl_length
    export char_array_entropy
    export aln_conservatio
    export eval_recomb

    #Export BLAST parsing
    export split_by_query
    export add_new_interval
    export covered_length
    export unique_length_fraction
    export get_coverage_maximising_hits

    #Export TaxonKitWrap
    export taxid_to_species
    export taxid_to_genus

    #HMM parsing
    export extract_rRNA_names

    #export workflow specific staff
    export parse_mapping
    export parse_dada_results
    export check_pair
    export fasta_to_dict
    export dict_to_fasta
    export read_salmon

    #export clustering output parsing functinos
    export read_jc_clusters
    export uc2jc
    export mergejc
    export mark_proper_size_clusters

    include("utils.jl"	)
    include("BLASTparsing.jl")
    include("TaxonKitWrap.jl")
    include("HMMutils.jl")
    include("UCparsing.jl")
    include("AlignmentMatrix.jl")
end


text:
getexperiment() = loadtable(experiment_file_name, indexcols = :experiment)

function register_experiment()
    experiment = getexperiment()
    @label startname
    println("Give a descriptive and unique name for the experiment")
    name = strip(readline(stdin))
    if isempty(name)
        @warn "The experiment must have a name"
        @goto startname
    end
    if name ∈ select(experiment, :experiment)
        @warn "You have already registered an experiment with that name:" select(experiment, :experiment)
        @goto startname
    end
    @label startdescription
    println("Describe the experiment as best you can")
    description = strip(readline(stdin))
    if isempty(description)
        @warn "The experiment must have a description"
        @goto startdescription
    end
    [(experiment = name, experiment_description = description)] |> CSV.write(experiment_file_name, append = true)
end


text:
using MPI
using ClimateMachine
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays

MPI.Initialized() || MPI.Init()
mpicomm = MPI.COMM_WORLD

# set up domain
topl = StackedBrickTopology(
    mpicomm,
    (0:10, 0:10, 0:3);
    periodicity = (false, false, false),
    boundary = ((1, 1), (1, 2), (1, 2)),
)

@show MPI.Comm_rank(mpicomm) length(topl.realelems)

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


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

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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 736856ce-f490-11eb-3349-057c86edfe7e
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using MLJ, MLJLinearModels, Plots, DataFrames, Random
end

# ╔═╡ c693088f-7f80-4cdd-b9b5-65a50da732ac
begin
    using MLCourse
    import MLCourse: poly, Polynomial
    MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ e1077092-f72a-42af-b6e0-a616f938cba8
using PlutoUI; PlutoUI.TableOfContents()

# ╔═╡ 7eb6a060-f948-4d85-881a-4909e74c15bd
md"# Validation Set Approach

In the following cell we define a `data_generator` that creates data sets with a polynomial relationship of degree 3 between the input x and the average output f(x).
"

# ╔═╡ 7dd7e9a7-9245-4c64-af0c-8f7d2f62b2bf
begin
    f(x) = 0.3 + 2x - 0.8x^2 - 0.4x^3
    function data_generator(; n = 500, seed = 12)
        rng = MersenneTwister(seed)
        x = randn(rng, n)
        DataFrame(x = x, y = f.(x) .+ randn(rng, n))
    end
end;

# ╔═╡ 5049407c-bf93-41a5-97bd-863a69d7016a
md"## Test Error Estimation

We will estimate here the test error with the validation set approach for fixed hyper-parameters."

# ╔═╡ 46827e21-0eea-4ec9-83ed-05e41dda1502
data1 = data_generator(seed = 1);

# ╔═╡ 6cfc32b3-40c0-4149-ba67-2ec67b3938f3
md"split seed = $(@bind split_seed Slider(1:20, show_value = true))"

# ╔═╡ 1320b9b0-d6fc-4e61-8d62-de3988b8b21d
idxs = shuffle(MersenneTwister(split_seed), 1:500) # shuffle the indices

# ╔═╡ 5dc21667-d94f-4c2a-b217-80b9a2262d8c
begin
	data1_train = data1[idxs[1:250], :]
	data1_test = data1[idxs[251:end], :]
	mach = machine(Polynomial(degree = 8) |> LinearRegressor(),
	               select(data1_train, :x), data1_train.y)
	fit!(mach, verbosity = 0)
	(training_error = rmse(predict(mach), data1_train.y),
	 test_error = rmse(predict(mach, select(data1_test, :x)), data1_test.y))
end

# ╔═╡ 1fe3f6d1-6336-417e-885e-de3fc5821bdf
let
	scatter(data1_train.x, data1_train.y, label = "training set")
	scatter!(data1_test.x, data1_test.y, label = "test set")
	xgrid = -4:.1:4
	pred = predict(mach, DataFrame(x = xgrid))
	plot!(xgrid, pred, label = "fit", ylim = (-6, 5), xlim = (-4, 4))
end

# ╔═╡ b79fae90-87f3-4f89-933b-d82e76c94d81
md"We see that the resulting curve and the test error depends strongly on the training set, i.e. the indices that are selected for the training set. The training error is always close to the irreducible error in this case."

# ╔═╡ 7d30f043-2915-4bba-aad5-cb7dbe76a3e5
md"## Hyper-parameter tuning

Now we will use a nested validation set approach to find the optimal degree in polynomial regression and estimate the test error. For this we will split the data into training, validation and test set. For each hyper-parameter value we will fit the parameters on the training set and estimate the test error on the validation set. When we have found the best hyper-parameters based on the test error computed on the validation set, we will estimate the test error of the best model using the test set."

# ╔═╡ ffffa3d0-bffb-4ff8-9966-c3312f952ac5
begin
	function data_split(data;
                        shuffle = false,
                        idx_train = 1:50,
                        idx_valid = 51:100,
                        idx_test = 101:500)
        idxs = if shuffle
                randperm(size(data, 1))
            else
                1:size(data, 1)
            end
        (train = data[idxs[idx_train], :],
         valid = data[idxs[idx_valid], :],
         test = data[idxs[idx_test], :])
    end
    function fit_and_evaluate(model, data)
        mach = fit!(machine(model, select(data.train, :x), data.train.y),
                    verbosity = 0)
        (train = rmse(predict(mach, select(data.train, :x)), data.train.y),
         valid = rmse(predict(mach, select(data.valid, :x)), data.valid.y),
         test = rmse(predict(mach, select(data.test, :x)), data.test.y))
    end
end

# ╔═╡ 2e02ac0d-c4d0-47ba-be57-445adeb6ab8b
losses1 = [fit_and_evaluate(Polynomial(; degree) |> LinearRegressor(),
	                        data_split(data1))
           for degree in 1:10]

# ╔═╡ 91eecd2b-af18-4a63-9684-28950e604d1a
let validlosses = getproperty.(losses1, :valid), i = argmin(validlosses)
    plot(1:10, validlosses, label = nothing)
    scatter!([i], [validlosses[i]], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ b45b32b1-8a65-4823-b3bb-f0b7cc57604b
md"For this seed of the random number generator the optimal degree (x coordinate of the red point in the figure above) found by the validation set approach is close to the actual degree of the data generating process. The estimated test loss of the model with degree $(argmin((x -> x.valid).(losses1))) is $(losses1[argmin((x -> x.valid).(losses1))].test)."

# ╔═╡ 23e91f72-0124-48d4-8b53-bd5da67b7ac7
losses1[4]

# ╔═╡ 4461d5b0-bc9b-4c89-aa76-70524a5caa7c
md"Let us repeat the process above for other splits of the data."

# ╔═╡ ac3c7e84-6c47-4dc1-b862-de4cfb05dad9
losses = [fit_and_evaluate(Polynomial(; degree) |> LinearRegressor(),
                           data_split(data1, shuffle = true))
          for degree in 1:10, _ in 1:20]

# ╔═╡ 26362233-b006-423d-8fb5-7cd9150405b4
let validlosses = getproperty.(losses, :valid), i = argmin(validlosses, dims = 1)
    plot(validlosses, label = nothing, ylims = (0, 10))
    scatter!((x -> x[1]).(i), validlosses[i], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 5cff19c7-3426-4d91-b555-f059f5f41886
md"Here we see a large variability of the optimal degree (x coordinates of all the points in the figure above).

In the following cell we compare the validation set estimate of the test error to the test set estimation of the test error. We can confirm that the validation set estimate is much lower than the test set estimate."

# ╔═╡ df16313f-76be-4d8b-88e7-7cc618ff49d4
let validlosses = getproperty.(losses, :valid),
	testlosses = getproperty.(losses, :test),
	winners = argmin(validlosses, dims = 1)
	(mean_valid_winner = mean(validlosses[winners]),
	 mean_test_winner = mean(testlosses[winners]))
end

# ╔═╡ 46bffe01-1a7c-405b-b3ec-bfe7570b8a3c
md"# Cross-Validation

Here we would like to find the optimal hyper-parameters with cross-validation.
We will not estimate the final test error of the model with the best hyper-parameters.

In the following cell we implement cross-validation."

# ╔═╡ b3aed705-5072-4d9b-bb8f-865ac1561bf6
begin
    function cross_validation_sets(idx, K)
        n = length(idx)
        r = n ÷ K
        [let idx_valid = idx[(i-1)*r+1:(i == K ? n : i*r)]
             (idx_valid = idx_valid, idx_train = setdiff(idx, idx_valid))
         end
         for i in 1:K]
    end
    function cross_validation(model, data; K = 5, shuffle = false)
		idxs = 1:size(data, 1)
		if shuffle
            idxs = Random.shuffle(idxs)
		end
        losses = [fit_and_evaluate(model,
                                   data_split(data;
									          idx_test = [], # no test set
									          idxs...)) # training and validation
                  for idxs in cross_validation_sets(idxs, K)]
        (train = mean(getproperty.(losses, :train)),
         valid = mean(getproperty.(losses, :valid)))
    end
end

# ╔═╡ 6916273a-d116-4173-acaa-2bcac1d1753b
cross_validation_sets(1:100, 4)

# ╔═╡ c14c9e49-a993-456e-9115-97da86f8e498
losses1_cv10 = [cross_validation(Polynomial(; degree) |> LinearRegressor(),
		                         data1,
		                         K = 10) for degree in 1:10]

# ╔═╡ 4b435518-2f12-4921-bb1f-fdd049ddfaed
let validlosses = getproperty.(losses1_cv10, :valid), i = argmin(validlosses)
    plot(1:10, validlosses, label = nothing)
    scatter!([i], [validlosses[i]], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 1e584a38-2fef-4877-87f6-92237d71c4b3
losses_cv10 = [cross_validation(Polynomial(; degree) |> LinearRegressor(),
                                data1,
                                K = 10, shuffle = true)
               for degree in 1:10, _ in 1:20]

# ╔═╡ a7c88b3f-92cb-4253-a889-c78683722c1d
let validlosses = getproperty.(losses_cv10, :valid), i = argmin(validlosses, dims = 1)
    plot(validlosses, label = nothing, ylims = (0.5, 2))
    scatter!((x -> x[1]).(i), validlosses[i], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ 0b81c3d5-277a-4fe6-889b-550e2f83c39d
md"With cross-validation we find often degree = 3 as the optimal one."

# ╔═╡ 5747d116-1e61-45f0-a87b-89372c6f270f
md"### Tuning Hyper-parameters with cross-validation and estimating the test error with the validation set approach

Instead of using all data for hyper-parameter search using cross-validation, we could set some data aside as a test set and use this to estimate the test error of the model with the optimal hyper-parameters."

# ╔═╡ 8204d9da-855d-4c66-b4dd-2c18a17b539b
losses2_cv10 = [cross_validation(Polynomial(; degree) |> LinearRegressor(),
		                         data1[1:100, :],
		                         K = 10) for degree in 1:10]

# ╔═╡ 96b3c070-0c72-44c7-ac98-ff6cb7e82380
md"The optimal degree found with cross-validation on the first 100 data points is:"

# ╔═╡ bae635b1-6eba-4ba4-80e9-547770612843
argmin((x -> x.valid).(losses2_cv10))

# ╔═╡ d9b04582-3a7a-4e19-a0da-d363810883f7
md"Let us refit the model with the optimal degree on all data that was used for cross-validation and estimate it's test error on data points 101 to 500."

# ╔═╡ 9042e811-9a34-42b8-8069-877f3b3f1e75
let
	model = Polynomial(degree = 4) |> LinearRegressor()
	mach = machine(model, select(data1[1:100, :], :x), data1[1:100, :y])
	fit!(mach, verbosity = 0)
	testset = data1[101:end, :]
	(test_error = rmse(predict(mach, select(testset, :x)), testset.y),)
end

# ╔═╡ e2658128-4053-4484-8e1e-229eceb755ab
md"# Resampling Strategies with MLJ"

# ╔═╡ eaacf529-1727-4743-941b-360c53088b1d
md"Instead of our own cross-validation function we can also use the builtin functions of MLJ. The `train_test_pairs` function does basically the same as our `cross_validation_sets` function."

# ╔═╡ 967a6e08-f0e6-45b7-988b-d0df237f3ddf
MLJ.MLJBase.train_test_pairs(CV(nfolds = 10), 1:100)

# ╔═╡ 535a11d9-833e-4237-949e-9a3e289c600b
md"`MLJ` has the very useful function `evaluate!`. Have a look at the Live docs to learn more about this function. In the following you see an example where we perform 10-fold cross-validation on the `rmse` measure to estimate the expected test error for a polynomial regressor of degree 4."

# ╔═╡ 8ae790b3-3987-4f41-8e21-adbb71081eb9
evaluate!(machine(Polynomial(degree = 4) |> LinearRegressor(),
	              select(data1, :x), data1.y),
          resampling = CV(nfolds = 10), measure = rmse)

# ╔═╡ 134df37f-b737-4c80-a5f9-1149aeec970c
md"We can use this now to find the best degree with 10-fold cross-validation for 100 different seeds. You will see that - even though not perfect - the hyper-parameters found with 10-fold cross-validation have lower variance and are usually closer to the true value than with the validation set approach."

# ╔═╡ abb71af3-8aae-4806-9d5a-d144c15d22ef
losses_mlj_cv10 = [evaluate!(machine(Polynomial(; degree) |> LinearRegressor(),
                                     select(data1, :x), data1.y),
                             resampling = CV(nfolds = 10, shuffle = true),
                             measure = rmse,
                             verbosity = 0).measurement[]
                   for degree in 1:10, seed in 1:20]

# ╔═╡ 6599d2b4-68f4-4c22-8e40-bf3722597692
let validlosses  = losses_mlj_cv10, i = argmin(validlosses, dims = 1)
    plot(validlosses, label = nothing, ylims = (0.5, 2))
    scatter!((x -> x[1]).(i), validlosses[i], label = nothing,
             xlabel = "degree", ylabel = "validation loss")
end

# ╔═╡ fc07009a-eeac-4c6b-9559-1e99ae68a6c3
md"MLJ implements the following resampling strategies. With `Holdout` we use the validation set approach."

# ╔═╡ 42511b31-6d70-495e-8362-01a29143b96e
subtypes(ResamplingStrategy)

# ╔═╡ 29683c99-6a6a-4f65-bea2-d592895d887e
md"# Model Tuning

Finding good hyper-parameters (tuning) is such an important step in the process of finding good machine learning models that there exist some nice utility functions to tune the hyper-parameters ([tuning section in the MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/)). In the cell below you see an example where the degree of polynomial regression is automatically tuned on a given dataset by performing 10-fold cross-validation on all degrees on a \"grid\", i.e. of all degrees from 1 to 17 are tested."

# ╔═╡ f93f20db-4fed-481f-b085-ca744b68fa8f
begin
    model = Polynomial() |> LinearRegressor()
    data2 = data_generator(seed = 2, n = 100)
    self_tuning_model = TunedModel(model = model,
                                   resampling = CV(nfolds = 10),
                                   tuning = Grid(),
                                   range = range(model, :(polynomial.degree),
                                                 values = 1:17),
                                   measure = rmse)
    self_tuning_mach = machine(self_tuning_model, select(data2, :x), data2.y)
	fit!(self_tuning_mach, verbosity = 0)
end

# ╔═╡ d72609a1-f93f-4ca6-8759-727662233e97
md"The self-tuned model actually found the true degree 3 for this data set, as we can see from the report of the tuned machine. The estimated test error can also be found in the report under `best_history_entry.measurement`."

# ╔═╡ 59245b3e-ddfc-46c4-ba44-86ce191672ae
report(self_tuning_mach)

# ╔═╡ f8b48f50-09cb-498a-89d1-9ac9b5722d0c
md"The result of the tuned machine looks quite good here: the reducible error (average distance between red and green curve) looks small."

# ╔═╡ 47a2c14b-a02d-43d7-ac8a-4d95a0d91fa8
let x = -3:.1:3
    scatter(data2.x, data2.y, label = "data", legend = :bottomleft)
    plot!(f, label = "generator", w = 2, xlims = (-3, 3))
    plot!(x, predict(self_tuning_mach, DataFrame(x = x)), label = "self tuning fit", w = 2)
end

# ╔═╡ cf3841ba-3963-4716-8fc9-8cce1dc4a5fa
md"# Nested Cross-Validation

If we want to tune the hyper-parameters with cross-validation and care about the estimated test error of the model found with cross-validation, we can use nested cross-validation, where hyper-parameters are optimized with cross-validation for each fold of the outer cross-validation. This happens when we `evaluate!` with cross-validation a self-tuning machine that itself uses cross-validation for hyper-parameter tuning. If we used `resampling = Holdout(fraction_train = 0.5)` we would use the validation set approach to estimate the test error with 50% of the data in the training set. If we used `resampling = Holdout(fraction_train = 0.5)` in the self-tuning machine, we would use the validation set approach for tuning the hyper-parameters.
"

# ╔═╡ 3e9c6f1c-4cb6-48d8-8119-07a4b03c2e4b
nested_cv = evaluate!(machine(self_tuning_model, select(data1, :x), data1.y),
                      resampling = CV(nfolds = 5), measure = rmse)

# ╔═╡ eba90af0-bf67-48f8-8d3c-f831b4fc289a
md"Have a look at the best models found for each fold:"

# ╔═╡ c8d26ab6-86ec-4b37-9540-0f785fd8cdc2
nested_cv.report_per_fold

# ╔═╡ ee89c448-1e69-4fd1-a4b8-7297a09f2685
md"# Exercises

## Conceptual
1. We review k-fold cross-validation.
    - Explain how k-fold cross-validation is implemented.
    - What are the advantages and disadvantages of k-fold cross-validation relative to:
        - The validation set approach?
        - LOOCV?
3. You are given a model with unknown hyper-parameters and the goal is to find the optimal hyper-parameters and estimate the test error of the optimal model. Suppose of colleague of yours argues that this should be done in the following way: \"Take all the data and run cross-validation to find the best hyper-parameters. Taking all data is better, than taking only a subset, because one has more data to estimate the parameters and estimate the test error on the validation sets. Once the best hyper-parameters are found: estimate the test error by running again cross-validation on all the data with the hyper-parameter fixed to the best value.\" Do you agree with your colleague? If not, explain where you disagree and where you agree.
2. Suppose you receive the following email of a colleague of yours. Write an answer to this email.
```
    Hi

    In my internship I am doing this research project with a company that
    manufactures constituents for medical devices. We need to assure high
    quality and detect the pieces that are not fully functional. To do so we
    perform 23 measurements on each piece. My goal is to find a machine
    learning method that detects defective pieces automatically. Both, my
    training set and my test set consist of 4000 pieces that were fine and
    30 defective ones.  I ran logistic regression on that data and found a
    training error of 0.1% and a test error of 0.3%, which is already pretty
    good, I think, no? But with kNN classification it is even better. For
    k = 7 I found a test error of 0.05% which is by far lower than all test
    errors I obtained with other values of k. I was really impressed. Now
    we have a method that predicts with 99.95% accuracy whether a piece is
    defective or not!

    Because you are taking this machine learning class now, I wanted to ask you
    for advice. Does it all make sense to you what I described? If not, do you
    have any suggestion to get even better results?

    Ã bientôt
    Jamie
```

## Applied
1. Take the `classification_data` in our notebook on
   \"flexibility and bias-variance-decomposition notebook\" and find with 10-fold
   cross-validation the optimal number ``k`` of neighbors of kNN
   classification, using the AUC measure. Hint: `MLJ` has the builtin function `auc`.
   Plot the validation AUC for ``k = 1, \ldots, 50``.
3. With the same data as in the previous exercise, use the `MLJ` function `estimate!` and a self tuning machine to estimate with the validation set approach the test error of kNN classifier whose hyper-parameter is tuned with 5 fold cross-validation. Use one quarter of the data for the test set.
2. In this exercise you apply our \"recipe for supervised learning\" (see slides). The goal is to predict the miles a car can drive per gallon fuel (mpg) as a function of its horsepower. You can download a dataset with `using OpenML; cars = DataFrame(OpenML.load(455))`. In the cleaning step we will remove all rows that contain missing values (you can use the function `dropmissing`). We select the machine learning methods polynomial regression and k nearest neighbors regression and we take as measure the `rmse`. Make sure to go trough the steps 2, 5, 9 of the recipe. Plot the predictions of the best method you found.
"

# 1. Perform k-nearest neighbors regression on data generated with our `data_generator`
#    defined in the first cell of this notebook and find the optimal number k of neighbors
#    with k-fold cross validation.

# ╔═╡ 0651292e-3f4e-4263-8235-4caa563403ec
MLCourse.footer()

# ╔═╡ Cell order:
# ╠═736856ce-f490-11eb-3349-057c86edfe7e
# ╟─7eb6a060-f948-4d85-881a-4909e74c15bd
# ╠═7dd7e9a7-9245-4c64-af0c-8f7d2f62b2bf
# ╟─5049407c-bf93-41a5-97bd-863a69d7016a
# ╠═46827e21-0eea-4ec9-83ed-05e41dda1502
# ╟─6cfc32b3-40c0-4149-ba67-2ec67b3938f3
# ╠═1320b9b0-d6fc-4e61-8d62-de3988b8b21d
# ╠═5dc21667-d94f-4c2a-b217-80b9a2262d8c
# ╟─1fe3f6d1-6336-417e-885e-de3fc5821bdf
# ╟─b79fae90-87f3-4f89-933b-d82e76c94d81
# ╟─7d30f043-2915-4bba-aad5-cb7dbe76a3e5
# ╠═ffffa3d0-bffb-4ff8-9966-c3312f952ac5
# ╠═2e02ac0d-c4d0-47ba-be57-445adeb6ab8b
# ╟─91eecd2b-af18-4a63-9684-28950e604d1a
# ╟─b45b32b1-8a65-4823-b3bb-f0b7cc57604b
# ╠═23e91f72-0124-48d4-8b53-bd5da67b7ac7
# ╟─4461d5b0-bc9b-4c89-aa76-70524a5caa7c
# ╠═ac3c7e84-6c47-4dc1-b862-de4cfb05dad9
# ╟─26362233-b006-423d-8fb5-7cd9150405b4
# ╟─5cff19c7-3426-4d91-b555-f059f5f41886
# ╠═df16313f-76be-4d8b-88e7-7cc618ff49d4
# ╟─46bffe01-1a7c-405b-b3ec-bfe7570b8a3c
# ╠═b3aed705-5072-4d9b-bb8f-865ac1561bf6
# ╠═6916273a-d116-4173-acaa-2bcac1d1753b
# ╠═c14c9e49-a993-456e-9115-97da86f8e498
# ╟─4b435518-2f12-4921-bb1f-fdd049ddfaed
# ╠═1e584a38-2fef-4877-87f6-92237d71c4b3
# ╟─a7c88b3f-92cb-4253-a889-c78683722c1d
# ╟─0b81c3d5-277a-4fe6-889b-550e2f83c39d
# ╟─5747d116-1e61-45f0-a87b-89372c6f270f
# ╠═8204d9da-855d-4c66-b4dd-2c18a17b539b
# ╟─96b3c070-0c72-44c7-ac98-ff6cb7e82380
# ╠═bae635b1-6eba-4ba4-80e9-547770612843
# ╟─d9b04582-3a7a-4e19-a0da-d363810883f7
# ╠═9042e811-9a34-42b8-8069-877f3b3f1e75
# ╟─e2658128-4053-4484-8e1e-229eceb755ab
# ╟─eaacf529-1727-4743-941b-360c53088b1d
# ╠═967a6e08-f0e6-45b7-988b-d0df237f3ddf
# ╟─535a11d9-833e-4237-949e-9a3e289c600b
# ╠═8ae790b3-3987-4f41-8e21-adbb71081eb9
# ╟─134df37f-b737-4c80-a5f9-1149aeec970c
# ╠═abb71af3-8aae-4806-9d5a-d144c15d22ef
# ╟─6599d2b4-68f4-4c22-8e40-bf3722597692
# ╟─fc07009a-eeac-4c6b-9559-1e99ae68a6c3
# ╠═42511b31-6d70-495e-8362-01a29143b96e
# ╟─29683c99-6a6a-4f65-bea2-d592895d887e
# ╠═f93f20db-4fed-481f-b085-ca744b68fa8f
# ╟─d72609a1-f93f-4ca6-8759-727662233e97
# ╠═59245b3e-ddfc-46c4-ba44-86ce191672ae
# ╟─f8b48f50-09cb-498a-89d1-9ac9b5722d0c
# ╟─47a2c14b-a02d-43d7-ac8a-4d95a0d91fa8
# ╟─cf3841ba-3963-4716-8fc9-8cce1dc4a5fa
# ╠═3e9c6f1c-4cb6-48d8-8119-07a4b03c2e4b
# ╟─eba90af0-bf67-48f8-8d3c-f831b4fc289a
# ╠═c8d26ab6-86ec-4b37-9540-0f785fd8cdc2
# ╟─ee89c448-1e69-4fd1-a4b8-7297a09f2685
# ╟─c693088f-7f80-4cdd-b9b5-65a50da732ac
# ╟─e1077092-f72a-42af-b6e0-a616f938cba8
# ╟─0651292e-3f4e-4263-8235-4caa563403ec

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse