text:
module HBM

include("Computation.jl")
using .Computation

include("Format.jl")
using .Format

include("FC.jl")
using .FC

export hbm, vhbm

# Determines an upper bound using Hard-Buddy Match method
# Author: Antara Hebbar
function hbm(m::Int64, s::Int64; output::Int64 = 2)
    output > 0 && printHeader(center("HARD BUDDY MATCH"))

    if !(m > s > 0)
        if output > 0
            printf("Hard Buddy Match does not apply", line = true)
            printEnd()
        end
        return 1 // 1
    elseif m % s == 0
        if output > 0
            printf("Since $m % $s = 0, f($m, $s) = 1.", line = true)
            printEnd()
        end
        return 1 // 1
    else
        V = Int64(ceil(2m / s))
        if V == 3
            #important variables for algorithm hbm
            d = m - s
            i = 0
            while 3 * d * i < s
                i += 1
            end
            k = i - 1
            a = s - (3 * d * k)
            Xsolutions = Xsol(a, d, k)

            if length(Xsolutions) == 0
                if output > 0
                    printf("No X solutions", line = true)
                    printEnd()
                end
                return 1 // 1
            elseif a == (2 * d)
                alpha = fc(m, s, output = 0)
                if output > 0
                    printf(
                        "Since a=2d, hbm($m, $s) = fc($m, $s). Therefore, upperbound of hbm($m, $s) is $alpha.",
                        line = true,
                    )
                    printEnd()
                end
                return alpha
            else
                X = minimum(Xsolutions)

                #special cases
                if a == 10 && d == 7 && m == (21 * k + 17) && s == (21 * k + 10)
                    alpha = (7 * k + 4) // (21 * k + 10)
                elseif a == 4 && d == 3 && m == (9 * k + 7) && s == (9 * k + 4)
                    alpha = (9 * k + 5) // (27 * k + 12)
                else
                    alpha = (d * k + X) // (3 * d * k + a)
                end

                if vhbm(a, d, k, X, output = output)
                    if output > 0
                        printfT(
                            "Deriving alpha",
                            "The derived value for X is $X. f($m, $s) ≤ (dk+X)/(3dk+a).",
                            "",
                            "Once plugging in d, k, a, and X to HBM algorithm, upper bound $alpha is derived.",
                        )
                        printEnd()
                    end
                    return alpha
                else
                    if output > 0
                        printf(
                            "VHBM does not verify the alpha derived by Hard-Buddy Match.",
                            line = true,
                        )
                        printEnd()
                    end
                    return 1 // 1
                end
            end
        else
            if output > 0
                printf(
                    "Hard-buddy match method requires that students are given 3-shares and 2-shares.
                    Since in f($m, $s), students are given $V and $(V-1) shares, HBM does not apply here.",
                    line = true,
                )
                printEnd()
            end
            return 1 // 1
        end
    end
end

# Verifies that with an input of a,d,k,X, f(3dk+a+d, 3dk+a) ≤ (dk+X)/(3dk+a), ouputs proof if wanted
function vhbm(a, d, k, X; output::Int64 = 2)
    #pre-processing stage: checks for a bad input
    pre = true
    alpha = (d * k + X) // (3 * d * k + a)
    muffins = 3 * d * k + a + d
    students = 3 * d * k + a
    V = Int64(ceil(2 * muffins / students))
    comden = 3 * d * k + a

    #formatting
    alphaS = formatFrac(alpha)

    output > 1 && printHeader("PREPROCESSING STAGE: ")

    if a == (2 * d) #use FC
        alpha = fc(muffins, students, output = output)
        if output > 1
            printf("Inputs are verified.")
        end
    end

    if X >= a // 3 && a > (2 * d + 1) && a < (3 * d - 1) #theorem 10.11.1 in muffins book, we assume a is in that interval
        if output > 1
            printf("Inputs are verified.")
        end
    end

    if X >= (a // 2) #theorem 10.11.2, we assume X<a/2
        if output > 1
            printf("Inputs are verified.")
        end
    end
    if X >= (a + d) // 4 #theorem 10.11.3, we assume X<(a+d)/4
        if output > 1
            printf("Inputs are verified.")
        end
    end

    if a < 1 || a > (3 * d)
        pre = false
        if output > 0
            printf("Bad input, VHBM failed. a ∉ {1...$(3*d)}.")
            printEnd()
        end
        return false
    end

    if X < a // 3
        pre = false
        if output > 0
            printf("Bad input, VHBM failed. X is less than a/3 ($(formatFrac(a//3))).")
            printEnd()
        end
        return false
    end

    output > 1 && printf("END OF PREPROCESSING STAGE.")

    if pre && V == 3

        if output > 1 #start of proof
            printHeader("START OF PROOF:")
            printf(
                "Assume, by way of contradiction, there is a ($muffins, $students) procedure where the smallest piece is > $alphaS.",
                "",
                " We are looking for a contradiction.",
            )
            printLine()

            printfT(
                "Solving for students and shares",
                "Every student has $(V) or $(V-1) shares.",
                "",
                "Let s₂ and s₃ equal the number of 2 and 3-students.",
                "",
                "s₂ = $(3*d*k+a-2*d). There are $(6*d*k+2*a-4*d) 2-shares.",
                "",
                "s₃ = $(2*d). There are $(6*d) 3-shares.",
            )

            printHeader("DEFINING TERMS FOR ALGORITHM: ")
            printfT(
                "Variables",
                "We will use the following variables for the HBM algorithm (where m=muffins and s=students):",
                "",
                "d = m-s = $d",
                "",
                "k is the largest number >= 0 such that 3dk > s. k = $k.",
                "",
                "a = s - 3dk = $a",
            )

            printfT(
                "Note",
                "A buddy of share size x is equal to 1-x. A match of share size y is equal to $muffins/$students-y.",
            )
        end

        #variables for interval diagram
        x = (d * k + a + d - 2 * X) // (comden)
        y = (d * k + d + X) // comden
        alphabuddy = 1 - alpha
        #formatting
        den = denominator(alpha)
        xS, yS, aB = formatFrac(x, den), formatFrac(y, den), formatFrac(alphabuddy)

        bij1, bij2 = formatFrac((2 * d * k + a - d - X) // comden, den),
        formatFrac((2 * d * k - d + 2 * X) // comden)
        int1, int2 = formatFrac(((d * k + a - X) // comden), den),
        formatFrac(((d * k + 2X) // comden), den)
        int3, int4 = formatFrac((2 * d * k + a - 2 * X) // comden, den),
        formatFrac((2 * d * k + X) // comden, den)
        bound1S = formatFrac((d * k + a - X) // comden, den)
        int5 = formatFrac((d * k + a // 2) // comden, den)

        if output > 1
            printHeader("INTERVAL DIAGRAM:")
            printf("The following diagram depicts what we know so far: ")
            println(
                "\n",
                interval(
                    ["(", "$alphaS"],
                    [")[", "$xS"],
                    ["](", "$yS"],
                    [")", "$aB"],
                    labels = ["$(V)s_$V $V-shs", "0", "$(V-1)s_$(V-1) $(V-1)-shs"],
                ),
            )
            printf("The assumption that X ≥ $(formatFrac(a//3)) ensures that the interval of $V-shares and $(V-1) shares do not intersect.")
            printLine()

            printfT(
                "Buddy/Match Bijections",
                "We define the following (where B signifies buddying and M signifies matching):",
                "",
                "M₀ = [$xS, $yS]",
                "B₀ = B(M₀) = [$bij1, $bij2]",
                "",
                "∀ ≤ i ≤ $(k-1), [Mᵢ = M(Bᵢ₋₁) = [($(d*k+a)+$d(i+1)-2X)/$comden, ($(d*k)+$(d)(i+1)+X)/$comden)]",
                "",
                "∀ ≤ i ≤ $(k-1), [Bᵢ = B(Mᵢ) = [($(2*d*k+a-X)-$(d)(i+1))/$comden, ($(2*d*k+2*X)-$(d)(i+1))/$comden)]",
            ) #printing bijections

            printfT(
                "Note",
                "Since M₀ is empty, so is the Mᵢ's' and Bᵢ's.",
                "Bκ₋₁ = [$int1, $int2] is also empty. We want this to be in the 3-shares.",
                "",
                "Therefore, we want $int2 < $(xS). This is true.",
            )

            #formatting issues
            z = 2 * a + 2 * d
            printHeader("INTERVAL ANALYSIS CONT.:")
            printf("Let the number of shares in [$alphaS, $bound1S] be 𝑧.")
            printfT(
                "HBM Algorithm",
                "Algorithm proves through buddy-match bijections that 𝑧 = 2a+2d = $(z).",
                "",
                "For complete derivation of 𝑧, read HBM chapter of Muffins Book (page 176-177).",
            )
            printf("The following picture captures what we know. ")

            printf("3-shares:")
            println(
                "\n",
                interval(
                    ["(", "$alphaS"],
                    [")[", "$int1"],
                    ["](", "$int2"],
                    [")[", "$xS"],
                    ["]", "$yS"],
                    labels = ["$z 3-shs", "0", "$(6*d-z) 3-shs", "0"],
                ),
            )
            println()

            printf("2-shares: ")
            println(
                "\n",
                interval(
                    ["(", "$yS"],
                    [")[", "$bij1"],
                    ["]", "$bij2"],
                    labels = ["$(6*d*k-10*d+2*a) 3-shs", "0"],
                ),
            )
            println(
                "\n",
                interval(
                    ["(", "$bij2"],
                    [")[", "$int3"],
                    ["](", "$int4"],
                    [")", "$aB"],
                    labels = ["$(6*d-z) 2-shs", "0", "$z 2-shs"],
                ),
            )

            printLine()

            casekeyS, casenum, endsum = formatFrac(k // 2), 1, 1
            midpointS = formatFrac(endsum * 1 // 2)
            parity, ivalue, method = "even", "k/2", "buddying"

            if k % 2 != 0
                casekeyS, casenum, endsum =
                    formatFrac((k + 1) // 2), 2, (3 * d * k + a + d) // comden
                midpointS = formatFrac(endsum * 1 // 2)
                parity, ivalue, method = "odd", "(k+1)/2", "matching"
            end

            sym1 = formatFrac((d * k + toFrac(casekeyS) * d + X) // comden, den)
            sym2 = formatFrac((d * k + a + d * toFrac(casekeyS) - X) // comden, den)

            printHeader("CASES:")
            printfT(
                "Case $casenum",
                "k ($k) is $parity. Let i = $ivalue = $casekeyS.",
                "",
                "Mᵢ = [$sym1, $sym2].",
                "",
                "The sum of the endpoints is $endsum, so the midpoint is $midpointS. So, Mκ/2 is symmetric by $method.",
            )

            printf("The following is what we know about the 3-shares: ")
            println(
                "\n",
                interval(
                    ["(", "$alphaS"],
                    ["|", "$int5"],
                    [")[", "$int1"],
                    ["](", "$int2"],
                    [")", "$xS"],
                    labels = ["$(a+d) 3-shs", "$(a+d) 3-shs", "0", "$(4*d-2*a) 3-shs"],
                ),
            )
            printLine()
            printfT(
                "Defining intervals",
                "We define the following intervals: ",
                "",
                "|J₁| = [$alphaS, $int5], |J₂| = [$int5, $int1], |J₃| = [$int2, $xS]",
                "",
                "J₁ = J₂ = $(a+d)      J₃ = $(4*d-2*a)",
            )

        end

        # Determine possible share combinations
        shares = combs(3, 3)
        workingcomb = []

        for (a, b, c) in comboTup(3, 3)
            NEsum = a * toFrac(int5) + b * toFrac(int1) + c * x
            TMsum = a * alpha + b * toFrac(int5) + c * toFrac(int2)
            if NEsum > muffins // students > TMsum && [a, b, c] in shares
                append!(workingcomb, [(a, b, c)])
            end
        end

        len = length(workingcomb)

        if len > 0 && (2, 0, 1) in workingcomb
            if output > 1
                sharecombs = [
                    "$a shares from J1, $b shares in J2, and $c shares in J3"
                    for (a, b, c) in workingcomb
                ]
                printHeader("POSSIBLE SHARES: ")
                printfT(
                    "Determining share combos",
                    "Given all possible share combos, the only shares from J1, J2, and J3 that are not forbidden are: ",
                    "",
                    sharecombs...,
                    "",
                )
            end

            output > 1 && printHeader("EQUATION ANALYSIS:")
            if len == 1
                if output > 1
                    printfT(
                        "One-variable equation",
                        "Since there is only 1 possible share combination, we look for a 2d solution in following set of equations:",
                        "",
                        "Equation 1 based on |J₁|+ |J₂| = a+d → 2y₁₁₃ = a+d = $(a+d)",
                        "",
                        "Equation 2 based on |J₃| = 4d-2a → y₁₁₃ = 4d-2a = $(4*d-2*a)",
                        "",
                        "Equation 3 based on s₃ = 2d → y₁₁₃ = 2d = $(2*d)",
                        "",
                        "Since y₁₁₃ has a solution from {0...$(2*d)}, VHBM verifies alpha",
                    )
                end
                return true

            elseif len == 2
                sol = false
                number = "Two"
                if (1, 1, 1) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ = a + d = $(a+d)",
                    "y₁₂₃ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃"
                elseif (0, 3, 0) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ = a + d = $(a+d)",
                    "3y₂₂₂ = a + d = $(a+d)",
                    "y₁₁₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₂₂₂ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₂₂₂"
                end

            elseif len == 3
                sol = false
                number = "Three"
                if (1, 1, 1) in workingcomb && (1, 0, 2) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ + y₁₃₃ = a + d = $(a+d)",
                    "y₁₂₃ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₁₃₃ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₁₃₃"
                elseif (1, 1, 1) in workingcomb && (0, 3, 0) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ = a + d = $(a+d)",
                    "y₁₂₃ + 3y₂₂₂ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₂₂₂ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₂₂₂"
                elseif (1, 2, 0) in workingcomb && (0, 3, 0) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₂ = a + d = $(a+d)",
                    "2y₁₂₂ + 3y₂₂₂ = a + d = $(a+d)",
                    "y₁₁₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₂ + y₂₂₂ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₂, y₂₂₂"
                elseif (1, 1, 1) in workingcomb && (0, 2, 1) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ = a + d = $(a+d)",
                    "y₁₂₃ + y₂₂₃ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ + y₂₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₂₂₃ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₂₂₃"
                end

            elseif len == 4
                sol = false
                number = "Four"
                if (2, 1, 0) in workingcomb &&
                   (1, 2, 0) in workingcomb &&
                   (0, 3, 0) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + 2y₁₁₂ + y₁₂₂ = a + d = $(a+d)",
                    "3y₂₂₂ = a + d = $(a+d)",
                    "no |J₃| equation",
                    "y₁₁₃ + y₁₁₂ + y₁₂₂ y₂₂₂ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₁₂, y₁₂₂, y₂₂₂"
                elseif (1, 1, 1) in workingcomb &&
                       (1, 2, 0) in workingcomb &&
                       (0, 3, 0) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ + y₁₂₂ = a + d = $(a+d)",
                    "3y₂₂₂ + 2y₁₂₂ + y₁₂₃ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₁₂₂ + y₂₂₂ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₁₂₂, y₂₂₂"
                elseif (1, 1, 1) in workingcomb &&
                       (1, 0, 2) in workingcomb &&
                       (0, 3, 0) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ + y₁₃₃ = a + d = $(a+d)",
                    "3y₂₂₂ + y₁₂₃ = a + d = $(a+d)",
                    "2y₁₃₃ + y₁₂₃= 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₁₃₃ + y₂₂₂ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₁₃₃, y₂₂₂"
                elseif (1, 1, 1) in workingcomb &&
                       (1, 0, 2) in workingcomb &&
                       (0, 2, 1) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ + y₁₃₃ = a + d = $(a+d)",
                    "y₁₂₃ + 2y₂₂₃ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ + 2y₁₃₃ + y₂₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₁₃₃ + y₂₂₃ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₁₃₃, y₂₂₃"
                elseif (1, 1, 1) in workingcomb &&
                       (0, 3, 0) in workingcomb &&
                       (0, 2, 1) in workingcomb
                    sol = true
                    eq1, eq2, eq3, eq4 = "2y₁₁₃ + y₁₂₃ = a + d = $(a+d)",
                    "3y₂₂₂ + y₁₂₃ + 2y₂₂₃ = a + d = $(a+d)",
                    "y₁₁₃ + y₁₂₃ + y₂₂₃ = 4d - 2a = $(4*d-2*a)",
                    "y₁₁₃ + y₁₂₃ + y₂₂₂ + y₂₂₃ = 2d = $(2*d)"
                    variables = "y₁₁₃, y₁₂₃, y₂₂₂, y₂₂₃"
                end
            else
                printfT(
                    "Equations",
                    "The system of equations formed by student combiations was unsolvable. VHBM failed.",
                )
                printEnd()
                return false
            end

        else
            if output > 1
                printf("There are no possible share combinations, vhbm failed")
            end
            return false
        end

        #conclusion
        if sol && len != 1
            if output > 1
                printfT(
                    "$number-variable equations",
                    "We look for a 2d solution in following set of equations:",
                    "",
                    "Equation 1 based on |J₁|+ |J₂| = a+d → ",
                    "",
                    "$eq1,  $eq2",
                    "",
                    "Equation 2 based on |J₃| = 4d-2a → ",
                    "",
                    "$eq3",
                    "",
                    "Equation 3 based on s₃ = 2d → ",
                    "",
                    "$eq4",
                    "",
                    "Since $variables has a solution from {0...$(2*d)}, VHBM verifies alpha",
                )
            end
            return true #verified for 2 variable solutions
        else
            if output > 1
                printf("No 2d solution found with derived 2-variable share combinations, vhbm failed.")
                printEnd()
            end
            return false
        end
    else
        if output > 0
            printf("Vhbm failed, hbm method needs students to have $V or $(V-1) shares.")
            printEnd()
        end
    end
end

# Using COND function and a,d,k values to derive candidates for X
function Xsol(a, d, k)
    Xsolutions = Array{Rational}(undef, 0)

    if d >= 1 && k >= 1 && a >= 1 && a <= (2 * d) #11.5.1: y₁₁₃ permitted, all other variables forbidden

        X = max((a + 2 * d) // 6, (2 * a - d) // 3) #corollary 11.7
        if COND(a, d, X)
            append!(Xsolutions, X)
        end

        if a != 1 || d != 1 #11.5.2: y₁₁₃, y₁₂₃ permitted, all others forb
            X = max((a + d) // 5, (2 * a - d) // 3, d // 2) #corollary 11.8
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a != (7 * d // 5) #11.5.3: y₁₁₃, y₂₂₂ permitted, all others forb
            X = max((3 * a - 2 * d) // 4, (a + 2 * d) // 6)#corollary 11.9
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a < (7 * d // 5) || a > (5 * d // 3) #11.5.4: y₁₁₃, y₁₂₂, y₂₂₂ permitted, all others forb
            X = max(a - d, (a + 2 * d) // 6) #corollary 11.10
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a != 1 || d != 1 #11.5.5: y₁₁₃, y₁₂₃, y₁₃₃ permitted, all others forb
            X = max((2 * a - d) // 3, d // 2) #corollary 11.11
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a < d || a > (7 * d) // 5 #11.5.6: y₁₁₃, y₁₂₃, y₂₂₂ permitted, all others forb
            X = max((3 * a - 2 * d) // 4, (a + d) // 5, d // 2) #corollary 11.12
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a != 1 || d != 1 #11.5.7: y₁₁₃, y₁₂₃, y₂₂₃ permitted, all others forb
            X = max((2 * a - d) // 3, (a + d) // 5) #corollarly 11.13
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a < (7 * d) // 5 #11.5.8: y₁₁₂, y₁₁₃, y₁₂₂, y₂₂₂ permitted, all others forb
            X = (a + 2 * d) // 6 #corollarly 11.14
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if d > a || a > (7 * d) // 3 #11.5.9: y₁₂₂, y₁₁₃, y₁₂₃, y₂₂₂ permitted, all others forb
            X = max(a - d, (a + d) // 5, d // 2) #corollarly 11.15
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a < d // 2 || a > (7 * d) // 5 #11.5.10: y₁₃₃, y₁₁₃, y₁₂₃, y₂₂₂ permitted, all others forb
            X = max((3 * a - 2 * d) // 4, d // 2) #corollarly 11.16
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a > d || a < d // 3 #11.5.11: y₁₃₃, y₁₁₃, y₁₂₃, y₂₂₃ permitted, all others forb
            X = max((2 * a - d) // 3, (a + 2 * d) // 8) #corollarly 11.17
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end

        if a < d || a > (7 * d) // 5 #11.5.12: y₁₁₃, y₁₂₃, y₂₂₂, y₂₂₃ permitted, all others forb
            X = max((a + d) // 5, (3 * a - 2 * d) // 4) #corollarly 11.17
            if COND(a, d, X)
                append!(Xsolutions, X)
            end
        end
    end

    return Xsolutions

    #else
    #A = Array{Int64}(undef, 0)
    #return A
    #end
end

# Helper function for Xsol -- determines if generated X meets specified condition
function COND(a, d, X)
    if a // 3 <= X && X < min(a // 2, (a + d) // 4)
        return true
    else
        return false
    end
end

# Helper function for VHBM -- returns all permutations from 0 to targ that sum to targ of size n
function combs(targ::Any, n::Any)
    openvec = Vector{Int64}(undef, 0)
    if n == 0 #if length is 0, return empty vector
        return openvec
    else

        sol = Vector{Int64}(undef, 0)
        fullsol = Vector{Vector{Int64}}(undef, 0)

        for i = 0:targ
            push!(sol, i)
        end

        fullsol = reverse.(digits.(0:targ^n-1, base = targ, pad = n))
        len = length(fullsol)

        finalsol = Array{Array{Int64,1}}(undef, 0)

        for i = 1:len
            element = fullsol[i]
            if sum(element) == targ
                push!(finalsol, element)
            end
        end
        othersolutions = Array{Array{Int64}}(undef, 0)

        othersolutions = targperm(targ, n)

        splitting(x, n) = [x[i:min(i + n - 1, length(x))] for i = 1:n:length(x)] #function will split targperm output into increments of length n

        append!(finalsol, splitting(othersolutions, n))
    end
end

# Helper function for combs -- returns permuations of target number and 0
function targperm(targ, n)
    littlesol = Array{Int64,1}(undef, 0)
    zero = Array{Int64,1}(undef, 0)
    bigsol = Array{Array{Int64,1}}(undef, 0)
    answer = Array{Int64,1}(undef, 0)

    for i = 1:n
        push!(zero, 0)
        push!(bigsol, zero) #generates n arrays consisting of zeros
    end

    #combinations of 0 and target number
    index = 1
    for i = 1:n
        littlesol = bigsol[i]
        for j = 1:n
            if index == j
                littlesol[j] = targ
                append!(answer, littlesol[j])
            else
                littlesol[j] = 0
                append!(answer, littlesol[j])
            end
        end
        index += 1

    end

    return answer
end

end


text:
# This implements src/snes/examples/tutorials/ex2.c from PETSc using the PETSc.jl package, using SNES
#
# This solves the equations sequentially
#
# Newton method to solve u'' + u^{2} = f, sequentially.

using PETSc, MPI, LinearAlgebra, SparseArrays, UnicodePlots

if ~MPI.Initialized()
    MPI.Init()
end

petsclib = PETSc.petsclibs[1]
PETSc.initialize(petsclib)

"""
    Computes initial guess
"""
function FormInitialGuess!(x)
    for i=1:length(x)
        x[i] = 0.50;
    end
end

"""
    Computes rhs forcing function
"""
function SetInitialArrays(n)
    h =  1.0/(n-1.0)
    F = zeros(n);
    xp = 0.0;
    for i=1:n
        v    = 6.0*xp + (xp+1.e-12)^6.0;
        F[i] = v;
        xp   = xp+h;
    end

    return F
end

"""
    Computes the residual f, given solution vector x
"""
function FormResidual!(cf,cx, args...)
    if typeof(cx) <: Ptr{Nothing}
        # When this routine is called from PETSc, cx is a pointer to a global vector
        # That's why we have to transfer it first to
        x   =   PETSc.unsafe_localarray(PETSc.scalartype(petsclib),cx)
    else
        x   = cx;
    end
    if typeof(cf) <: Ptr{Nothing}
        f   =   PETSc.unsafe_localarray(PETSc.scalartype(petsclib),cf)
    else
        f   = cf;
    end
    n       =   length(x);
    xp      =   LinRange(0.0,1.0, n);
    F       =   6.0.*xp .+ (xp .+1.e-12).^6.0;      # define source term function

    dx      =   1.0/(n-1.0);
    f[1]    =   x[1] - 0.0;
    for i=2:n-1
        f[i] = (x[i-1] - 2.0*x[i] + x[i+1])/dx^2 + x[i]*x[i] - F[i]
    end
    f[n]    =   x[n] - 1.0;
    Base.finalize(x)
    Base.finalize(f)

end

"""
    Computes the jacobian, given solution vector x
"""
function FormJacobian!(cx, args...)

    if typeof(cx) <: Ptr{Nothing}
        x   =   PETSc.unsafe_localarray(PETSc.scalartype(petsclib),cx)
    else
        x   =   cx;
    end

    J   =   args[1];        # preconditioner = args[2], in case we want it to be different from J
    n   =   length(x);
    dx  =   1.0/(n-1.0);

    # interior points (hand-coded jacobian)
    for i=2:n-1
        J[i,i-1] = 1.0/dx^2;
        J[i,i  ] = -2.0/dx^2 + 2.0*x[i];
        J[i,i+1] = 1.0/dx^2;
    end

    # boundary points
    J[1,1] = 1.0;
    J[n,n] = 1.0;

    if typeof(J) <: PETSc.AbstractMat
        PETSc.assemble(J);  # finalize assembly

    end

    Base.finalize(x)
end

# ==========================================
# Main code

# Compute initial solution
n   =   21;
F   =   SetInitialArrays(n);
x   =   zeros(n);

FormInitialGuess!(x);

# Compute initial jacobian using a julia structure to obtain the nonzero structure
# Note that we can also obtain this structure in a different manner
Jstruct  = zeros(n,n);
FormJacobian!(x, Jstruct);                              # jacobian in julia form
Jsp      =   sparse(Float64.(abs.(Jstruct) .> 0))       # sparse julia, with 1.0 in nonzero spots
PJ       =   PETSc.MatSeqAIJ(Jsp);                      # transfer to

# Setup snes
x_s = PETSc.VecSeq(x);                  # solution vector
res = PETSc.VecSeq(F);     # residual vector

S = PETSc.SNES{Float64}(PETSc.petsclibs[1],MPI.COMM_WORLD;
        snes_rtol=1e-12,
        snes_monitor=false,
        snes_converged_reason=false);
PETSc.setfunction!(S, FormResidual!, res)
PETSc.setjacobian!(S, FormJacobian!, PJ, PJ)

# solve
PETSc.solve!(x_s, S);

# Extract & plot solution
x_sol = x_s.array;                  # convert solution to julia format
FormResidual!(res.array,x_sol)      # just for checking, compute residual

@show norm(res.array)

PETSc.finalize(petsclib)

# plot solution in REPL
lineplot(LinRange(0,1,n),x_sol,xlabel="width",ylabel="solution")

ext: jl


text:
@testset "Notebook Function Tests" begin

  @test isdefined(Julz, :notebook) == true

end


text:
# include(joinpath("toolkit","runtests.jl"))
# include(joinpath("geometric_center", "runtests.jl"))
# include(joinpath("distance", "runtests.jl"))
# include(joinpath("principal_axes","runtests.jl"))
include(joinpath("angle","runtests.jl"))


text:
lines = readlines("inputs/6-10/09.txt")
nums = [[parse(Int, x) for x in line] for line in lines]

# gets position (r,c) of 2D vec grid, or 9 if out of bounds
function getat(r,c,grid,fallback=9)
  h = length(grid)
  w = length(grid[1])
  if r < 1 || c < 1 || r > h || c > w
    return fallback
  end
  return nums[r][c]
end

function part1(nums)
  t = 0
  for r in 1:length(nums), c in 1:length(nums[1])
    val = nums[r][c]
    if val < getat(r-1,c, nums) &&  val < getat(r+1,c, nums) &&  val < getat(r,c-1, nums) &&  val < getat(r,c+1, nums)
      t += val + 1
    end
  end
  t
end

# unique identifier for an x,y pair
f(x, y) = (100 * x) + y

function part2(nums)
  g = SimpleGraph(10000)

  for r in 1:length(nums)
    for c in 1:length(nums[1])
      if nums[r][c] == 9
        continue
      end
      for (dc, dr) in [(0, -1), (0, 1), (-1, 0), (1, 0)]
        # gets nums[r + dr][c + dc] or 9 if out of bounds
        if getat(r + dr, c + dc, nums, 9)  != 9
          add_edge!(g, f(r, c), f(r + dr, c + dc))
        end
      end
    end
  end

  lens = [length(cc) for cc in connected_components(g)] |> sort
  lens[end] * lens[end-1] * lens[end-2]
end


text:
# This file contains functions that create properties for the desired least-square spline

"""
```
default_slm_kwargs!
```

adds the default keyword arguments for an SLM object.
"""
function default_slm_kwargs!(kwargs::Dict, T::Type)

    # degree
    if !haskey(kwargs, :degree)
        kwargs[:degree] = 3
    end

    # scaling
    if !haskey(kwargs, :scaling)
        kwargs[:scaling] = true
    end

    # knots
    if !haskey(kwargs, :knots)
        kwargs[:knots] = 6
    end

    # C2 function
    if !haskey(kwargs, :C2)
        kwargs[:C2] = true
    end

    # Regularization parameter
    if !haskey(kwargs, :λ)
        kwargs[:λ] = convert(T, 1e-4)
    end

    # End point values
    if !haskey(kwargs, :left_value)
        kwargs[:left_value] = NaN
    end

    if !haskey(kwargs, :right_value)
        kwargs[:right_value] = NaN
    end

    # Global min and max
    if !haskey(kwargs, :min_value)
        kwargs[:min_value] = NaN
    end

    if !haskey(kwargs, :max_value)
        kwargs[:max_value] = NaN
    end

    if !haskey(kwargs, :min_max_sample_points)
        kwargs[:min_max_sample_points] = convert(Vector{T}, [.017037, .066987, .14645, .25, .37059,
                                                             .5, .62941, .75, .85355, .93301, .98296])
    end

    # Monotonicity
    if !haskey(kwargs, :increasing)
        kwargs[:increasing] = false
    end

    if !haskey(kwargs, :decreasing)
        kwargs[:decreasing] = false
    end

    if !haskey(kwargs, :increasing_intervals)
        kwargs[:increasing_intervals] = Matrix{T}(undef, 0, 0)
    end

    if !haskey(kwargs, :decreasing_intervals)
        kwargs[:decreasing_intervals] = Matrix{T}(undef, 0, 0)
    end

    # Curvature
    if !haskey(kwargs, :concave_up)
        kwargs[:concave_up] = false
    end

    if !haskey(kwargs, :concave_down)
        kwargs[:concave_down] = false
    end

    if !haskey(kwargs, :concave_up_intervals)
        kwargs[:concave_up_intervals] = Matrix{T}(undef, 0, 0)
    end

    if !haskey(kwargs, :concave_down_intervals)
        kwargs[:concave_down_intervals] = Matrix{T}(undef, 0, 0)
    end

    # Initial guess for spline's coefficients
    if !haskey(kwargs, :init)
        kwargs[:init] = Vector{T}(undef, 0)
    end

    # Use sparse matrices when constructing the
    # objective and constraints
    if !haskey(kwargs, :use_sparse)
        kwargs[:use_sparse] = false
    end

    # Use LLSModel rather than ADNLSModel (from NLPModels) to
    # define the constrained linear least squares problem.
    # This usually creates speed gains
    if !haskey(kwargs, :use_lls)
        kwargs[:use_lls] = true
    end

    return kwargs
end

"""
```
scale_problem!(x::AbstractVector, y::AbstractVector, kwargs)
```

scales y to minimize numerical issues when constructing splines.
"""
function scale_problem!(x::AbstractVector, y::AbstractVector, kwargs::Dict)
    scaling = haskey(kwargs, :scaling) ? kwargs[:scaling] : false
    kwargs[:scaling] = scaling # in case scaling didn't exist already

    if scaling
        # scale y so that the minimum value is 1/phi, and the maximum value phi.
        # where phi is the golden ratio, so phi = (sqrt(5) + 1)/2 = 1.6180...
        # Note that phi - 1/phi = 1, so the new range of yhat is 1. (Note that
        # this interval was carefully chosen to bring y as close to 1 as
        # possible, with an interval length of 1.)
        #
        # The transformation is:
        #   ŷ = y * y_scale + y_shift
        ϕ⁻¹ = (sqrt(5.) - 1.) / 2.;

        # Shift and scale are determined from min and max of y
        ymin = minimum(y)
        ymax = maximum(y)

        y_scale = 1. / (ymax - ymin)

        if isinf(y_scale)
            # If data passed is constant, then no scaling needed
            y_scale = 1.
        end

        y_shift = ϕ⁻¹ - y_scale * ymin

        ŷ = y * y_scale .+ y_shift

        # finally, shift/scale each part of the prescription as necessary.
        # derivative information need only be scaled of course. And since
        # y_scale will always be positive, monotonicity signs
        # and curvature signs will remain unchanged.

        # Those constraints that will be left untouched:
        #  C2, concave_down, concave_up
        #  decreasing, degree, increasing, knots
        #
        # All other constraints will potentially be affected.

        # left_value
        if haskey(kwargs, :left_value)
            kwargs[:left_value] *= y_scale
            kwargs[:left_value] += y_shift
        end

        # right_value
        if haskey(kwargs, :right_value)
            kwargs[:right_value] *= y_scale
            kwargs[:right_value] += y_shift
        end

        # left_value
        if haskey(kwargs, :min_value)
            kwargs[:min_value] *= y_scale
            kwargs[:min_value] += y_shift
        end

        # max_value
        if haskey(kwargs, :max_value)
            kwargs[:max_value] *= y_scale
            kwargs[:max_value] += y_shift
        end
    else
        y_shift = 0.
        y_scale = 1.
        ŷ = copy(y)  # Copy here to avoid accidentally changing original y values
    end

    kwargs[:y_scale] = y_scale
    kwargs[:y_shift] = y_shift

    return ŷ
end

"""
```
construct_design_matrix(x::AbstractVector{T}, knots::AbstractVector{T},
    dknots::AbstractVector{T}, x_bins::AbstractVector{Int}, nₓ::Int, nk::Int, nc::Int;
    use_sparse::Bool = false) where {T <: Real}
```

constructs the design matrix used to create an SLM.
"""
function construct_design_matrix(x::AbstractVector{T}, knots::AbstractVector{T},
                                 dknots::AbstractVector{T}, x_bins::AbstractVector{Int},
                                 nₓ::Int, nk::Int, nc::Int; use_sparse = false) where {T <: Real}

    # Create values used in the design matrix
    t  = (x - knots[x_bins]) ./ dknots[x_bins]
    t² = t.^2
    t³ = t.^3
    s² = (1. .- t) .^ 2
    s³ = (1. .- t) .^ 3

    vals = [3. .* s² .- 2. .* s³;
            3. .* t² .- 2. .* t³;
            (s² - s³) .* dknots[x_bins];
            (t³ - t²) .* dknots[x_bins]]

    # Coefficients are stored in two blocks,
    # first nk function values, then nk derivatives
    # Must use vector of CartesianIndex for the subscripts (unlike in Matlab)
    Mdes = if use_sparse
        sparse_accumarray([CartesianIndex(i, j) for (i, j) in
                    zip(repeat(collect(1:nₓ), 4),
                        [x_bins; x_bins .+ 1; nk .+ x_bins; (nk + 1) .+ x_bins])],
                   vals, (nₓ, nc))
    else
        accumarray([CartesianIndex(i, j) for (i, j) in
                    zip(repeat(collect(1:nₓ), 4),
                        [x_bins; x_bins .+ 1; nk .+ x_bins; (nk + 1) .+ x_bins])],
                   vals, (nₓ, nc))
    end

    return Mdes
end

"""
```
construct_regularizer(dknots::AbstractVector{T}, nk::Int;
    use_sparse::Bool = false) where {T <: Real}
```

constructs the regularizer equations used to fit the least-squares spline
when constructing an SLM.
"""
function construct_regularizer(dknots::AbstractVector{T}, nk::Int;
                               use_sparse::Bool = false) where {T <: Real}

    # We are integrating the piecewise linear f''(x)
    # as a quadratic form in terms of the (unknown) second
    # derivatives at the knots
    if use_sparse
        sf = spzeros(T, nk, nk)
        sd = spzeros(T, nk, nk)
    else
        sf = zeros(T, nk, nk)
        sd = zeros(T, nk, nk)
    end
    Mreg = zeros(T, nk, nk) # need to take a Cholesky of this matrix

    Mreg[1, 1] = dknots[1] / 3.
    Mreg[1, 2] = dknots[1] / 6.
    Mreg[nk, nk - 1] = dknots[end] / 6.
    Mreg[nk, nk] = dknots[end] / 3.
    for i in 2:(nk - 1)
        Mreg[i, i - 1] = dknots[i - 1] / 6.
        Mreg[i, i] = (dknots[i - 1] + dknots[i]) / 3.
        Mreg[i, i + 1] = dknots[i] / 6.
    end
    Mreg = use_sparse ? sparse(cholesky(Mreg).U) : cholesky(Mreg).U # Matrix square root, cholesky is easy way to do this

    # Write second derivativeas as a function of
    # the function values and first derivatives
     for i in 1:(nk - 1)
        sf[i, i] = -(6. / dknots[i]^2)
        sf[i, i + 1] = 6. / dknots[i]^2
        sd[i, i] = -4. / dknots[i]
        sd[i, i + 1] = -2. / dknots[i]
    end
    sf[nk, nk - 1] = 6. / dknots[end]^2
    sf[nk, nk] = -6. / dknots[end]^2
    sd[nk, nk - 1] = 2 / dknots[end]
    sd[nk, nk] = 4. / dknots[end]
    Mreg = Mreg * hcat(sf, sd)

    # Scale the regularizer before applied to the regularizing parameter
    Mreg ./= opnorm(Mreg, 1)

    return Mreg
end

"""
```
C2_matrix(nk::Int, nc::Int, dknots::AbstractVector{T},
    Meq::AbstractMatrix{T}, rhseq::AbstractVector{T},
    use_sparse::Bool = false) where {T <: Real}
```

constructs the matrix of equations that enforce twice-continuous differentiability.
"""
function C2_matrix(nk::Int, nc::Int, dknots::AbstractVector{T}, Meq::AbstractMatrix{T}, rhseq::AbstractVector{T};
                   use_sparse::Bool = false) where {T <: Real}
    MC2 = use_sparse ? spzeros(T, nk - 2, nc) : zeros(T, nk - 2, nc)
    for i in 1:(nk - 2)
        MC2[i, i] = 6. / dknots[i]^2
        MC2[i, i + 1] = -6. / dknots[i]^2 + 6. / dknots[i + 1]^2
        MC2[i, i + 2] = -6. / dknots[i + 1]^2

        MC2[i, nk + i] = 2. / dknots[i]
        MC2[i, nk + i + 1] = 4. / dknots[i] + 4. / dknots[i + 1]
        MC2[i, nk + i + 2] =  2. / dknots[i]
    end

    return vcat(Meq, MC2), vcat(rhseq, zeros(T, nk - 2))
end

"""
```
function monotone_increasing!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
    nk::Int) where {T <: Real}
```

adds the properties of a monotone increasing spline.
"""
function monotone_increasing!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                              nk::Int)
    push!(monotone_settings, (knotlist = (1, nk), increasing = true))
    return nk - 1
end

"""
```
function increasing_intervals_info!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                                   knots::AbstractVector{T}, increasing_intervals::AbstractMatrix{T}, nk::Int = length(knots)) where {T <: Real}
```

adds information about interval(s) over which the spline should be increasing
"""
function increasing_intervals_info!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                                   knots::AbstractVector{T}, increasing_intervals::AbstractMatrix{T}, nk::Int = length(knots)) where {T <: Real}
    out = 0
    for i in 1:size(increasing_intervals, 1)
        knotlist = (max(1, findlast(knots .<= increasing_intervals[i, 1])),
                    min(nk, 1 + findlast(knots .< increasing_intervals[i, 2])))
        push!(monotone_settings, (knotlist = knotlist, increasing = true))
        out += knotlist[2] - knotlist[1]
    end

    return out
end

"""
```
function monotone_decreasing!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
    nk::Int) where {T <: Real}
```

adds the properties of a monotone decreasing spline.
"""
function monotone_decreasing!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                              nk::Int) where {T <: Real}
    push!(monotone_settings, (knotlist = (1, nk), increasing = false))
    return nk - 1
end

"""
```
function decreasing_intervals_info!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                                   knots::AbstractVector{T}, decreasing_intervals::AbstractMatrix{T}, nk::Int = length(knots)) where {T <: Real}
```

adds information about interval(s) over which the spline should be decreasing
"""
function decreasing_intervals_info!(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                                   knots::AbstractVector{T}, decreasing_intervals::AbstractMatrix{T}, nk::Int = length(knots)) where {T <: Real}
    out = 0
    for i in 1:size(decreasing_intervals, 1)
        knotlist = (max(1, findlast(knots .<= decreasing_intervals[i, 1])),
                    min(nk, 1 + findlast(knots .< decreasing_intervals[i, 2])))
        push!(monotone_settings, (knotlist = knotlist, increasing = false))
        out += knotlist[2] - knotlist[1]
    end

    return out
end

"""
```
function construct_monotoncitiy_matrix(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                                      nc::Int, nk::Int, dknots::AbstractVector{T}, total_intervals::Int,
                                      Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
    use_sparse::Bool = false) where {T <: Real}
```
The function must be monotone between
knots j and j + 1. The direction over
that interval is specified. The constraint
system used comes from Fritsch & Carlson, see here:

http://en.wikipedia.org/wiki/Monotone_cubic_interpolation

Define delta = (y(i+1) - y(i))/(x(i+1) - x(i))
Thus delta is the secant slope of the curve across
a knot interval. Further, define alpha and beta as
the ratio of the derivative at each end of an
interval to the secant slope.

 alpha = d(i)/delta
 beta = d(i+1)/delta

Then we have an elliptically bounded region in the
first quadrant that defines the set of monotone cubic
segments. We cannot define that elliptical region
using a set of linear constraints. However, by use
of a system of 7 linear constraints, we can form a
set of sufficient conditions such that the curve
will be monotone. There will be some few cubic
segments that are actually monotone, yet lie outside
of the linear system formed. This is acceptable,
as our linear approximation here is a sufficient
one for monotonicity, although not a necessary one.
It merely says that the spline may be slightly over
constrained, i.e., slightly less flexible than is
absolutely necessary. (So?)

The 7 constraints applied for an increasing function
are (in a form that lsqlin will like):

   -delta          <= 0
   -alpha          <= 0
   -beta           <= 0
   -alpha + beta   <= 3
    alpha - beta   <= 3
    alpha + 2*beta <= 9
    2*alpha + beta   <= 9

Multiply these inequalities by (y(i+1) - y(i)) to
put them into a linear form.
"""
function construct_monotonicity_matrix(monotone_settings::Vector{NamedTuple{(:knotlist, :increasing), Tuple{Tuple{Int, Int}, Bool}}},
                                      nc::Int, nk::Int, dknots::AbstractVector{T}, total_intervals::Int,
                                      Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
                                       use_sparse::Bool = false) where {T <: Real}
    L = length(monotone_settings)
    M = use_sparse ? spzeros(T, 7 * total_intervals, nc) : zeros(T, 7 * total_intervals, nc)
    n = 0

    for i in 1:L
        for j in monotone_settings[i][:knotlist][1]:(monotone_settings[i][:knotlist][2] - 1)
            if monotone_settings[i][:increasing]
                M[n + 1, j]          = 1.
                M[n + 1, j + 1]      = -1
                M[n + 2, nk + j]     = -1.
                M[n + 3, nk + j + 1] = -1.

                M[n + 4, j]          = 3.
                M[n + 4, j + 1]      = -3.
                M[n + 4, nk + j]     = -1. * dknots[j]
                M[n + 4, nk + j + 1] = 1. * dknots[j]

                M[n + 5, j]          = 3.
                M[n + 5, j + 1]      = -3.
                M[n + 5, nk + j]     = 1. * dknots[j]
                M[n + 5, nk + j + 1] = -1. * dknots[j]

                M[n + 6, j]          = 9.
                M[n + 6, j + 1]      = -9.
                M[n + 6, nk + j]     = 1. * dknots[j]
                M[n + 6, nk + j + 1] = 2. * dknots[j]

                M[n + 7, j]          = 9.
                M[n + 7, j + 1]      = -9.
                M[n + 7, nk + j]     = 2. * dknots[j]
                M[n + 7, nk + j + 1] = 1. * dknots[j]
            else
                M[n + 1, j]          = -1.
                M[n + 1, j + 1]      = 1
                M[n + 2, nk + j]     = 1.
                M[n + 3, nk + j + 1] = 1.

                M[n + 4, j]          = -3.
                M[n + 4, j + 1]      = 3.
                M[n + 4, nk + j]     = 1. * dknots[j]
                M[n + 4, nk + j + 1] = -1. * dknots[j]

                M[n + 5, j]          = -3.
                M[n + 5, j + 1]      = 3.
                M[n + 5, nk + j]     = -1. * dknots[j]
                M[n + 5, nk + j + 1] = 1. * dknots[j]

                M[n + 6, j]          = -9.
                M[n + 6, j + 1]      = 9.
                M[n + 6, nk + j]     = -1. * dknots[j]
                M[n + 6, nk + j + 1] = -2. * dknots[j]

                M[n + 7, j]          = -9.
                M[n + 7, j + 1]      = 9.
                M[n + 7, nk + j]     = -2. * dknots[j]
                M[n + 7, nk + j + 1] = -1. * dknots[j]
            end
            n += 7
        end
    end

    return  vcat(Mineq, M), vcat(rhsineq, zeros(T, size(M, 1)))
end

"""
```
set_left_value(left_value::T, nc:Int, M::AbstractMatrix{T}, rhseq::AbstractVector{T};
    use_sparse::Bool = false) where {T <: Real}
```
"""
function set_left_value(left_value::T, nc::Int, M::AbstractMatrix{T}, rhseq::AbstractVector{T};
                        use_sparse::Bool = false) where {T <: Real}
    M_add = use_sparse ? spzeros(T, 1, nc) : zeros(T, 1, nc)
    M_add[1] = 1.
    return vcat(M, M_add), vcat(rhseq, left_value)
end

"""
```
set_right_value(right_value::T, nc::Int, nk::Int, M::AbstractMatrix{T}, rhseq::AbstractVector{T};
    use_sparse::Bool = false) where {T <: Real}
```
"""
function set_right_value(right_value::T, nc::Int, nk::Int, M::AbstractMatrix{T}, rhseq::AbstractVector{T};
                         use_sparse::Bool = false) where {T <: Real}
    M_add = use_sparse ? spzeros(T, 1, nc) : zeros(T, 1, nc)
    M_add[nk] = 1.
    return vcat(M, M_add), vcat(rhseq, right_value)
end

"""
```
set_min_value(min_value::T, nk::Int, nc::Int, dknots::AbstractVector{T}, Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
    sample_points::AbstractVector{T} = [.017037, .066987, .14645, .25, .37059,
                                        .5, .62941, .75, .85355, .93301, .98296],
    use_sparse::Bool = false) where {T <: Real}
```

creates the equations enforcing a maximum value for the spline.
"""
function set_min_value(min_value::T, nk::Int, nc::Int, dknots::AbstractVector{T}, Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
                       sample_points::AbstractVector{T} = [.017037, .066987, .14645, .25, .37059,
                                                           .5, .62941, .75, .85355, .93301, .98296],
                       use_sparse::Bool = false) where {T <: Real}

    # The default sample points are Chebyshev nodes
    nsamp = length(sample_points)
    ntot = nk + (nk - 1) * nsamp
    Mmin = use_sparse ? spzeros(T, ntot, nc) : zeros(T, ntot, nc)

    # Constrain values at knots
    for i in 1:nk # This is the same as Min[1:nk, 1:nk] = -Matrix{T}(I, nk, nk)
        Mmin[i, i] = -1.
    end

    # Intermediate sample points
    t² = sample_points .^ 2
    t³ = sample_points .^ 3
    s² = (1. .- sample_points) .^ 2
    s³ = (1. .- sample_points) .^ 3
    vals = hcat(3. .* s² .- 2. .* s³,
            3. .* t² .- 2. .* t³,
             s² - s³, t³ .- t²)

    for j in 1:(nk - 1)
        Mmin[(1:nsamp) .+ ((j - 1) * nsamp + nk), j .+ [0, 1, nk, nk + 1]] =
            -vals * Diagonal([1., 1., dknots[j], dknots[j]])
    end

    return vcat(Mineq, Mmin), vcat(rhsineq, fill(-min_value, ntot))
end

"""
```
set_max_value(max_value::T, nk::Int, nc::Int, dknots::AbstractVector{T}, Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
    sample_points::AbstractVector{T} = [.017037, .066987, .14645, .25, .37059,
                                        .5, .62941, .75, .85355, .93301, .98296],
    use_sparse:Bool = false) where {T <: Real}
```

creates the equations enforcing a maximum value for the spline.
"""
function set_max_value(max_value::T, nk::Int, nc::Int, dknots::AbstractVector{T}, Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
                       sample_points::AbstractVector{T} = [.017037, .066987, .14645, .25, .37059,
                                                           .5, .62941, .75, .85355, .93301, .98296],
                       use_sparse::Bool = false) where {T <: Real}

    # The default sample points are Chebyshev nodes
    nsamp = length(sample_points)
    ntot = nk + (nk - 1) * nsamp
    Mmax = use_sparse ? spzeros(T, ntot, nc) : zeros(T, ntot, nc)

    # Constrain values at knots
    for i in 1:nk # This is the same as Mmax[1:nk, 1:nk] = Matrix{T}(I, nk, nk)
        Mmax[i, i] = 1.
    end

    # Intermediate sample points
    t² = sample_points .^ 2
    t³ = sample_points .^ 3
    s² = (1. .- sample_points) .^ 2
    s³ = (1. .- sample_points) .^ 3
    vals = hcat(3. .* s² .- 2. .* s³,
            3. .* t² .- 2. .* t³,
            -s³ + s², t³ .- t²)

    # Create equations
    for j in 1:(nk - 1)
        Mmax[(1:nsamp) .+ ((j - 1) * nsamp + nk), j .+ [0, 1, nk, nk + 1]] =
            vals * Diagonal([1, 1, dknots[j], dknots[j]])
    end

    return vcat(Mineq, Mmax), vcat(rhsineq, fill(max_value, ntot))
end

"""
```
function concave_up_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Tuple{T, T}}}}) where {T <: Real}
```

adds the properties of an everywhere concave up spline.
"""
function concave_up_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}}) where {T <: Real}
    push!(curvature_settings, (concave_up = true, range = Vector{T}(undef, 0)))
end

"""
```
function concave_up_intervals_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}},
                                    concave_up_intervals::AbstractMatrix{T}) where {T <: Real}
```

adds information about interval(s) over which the spline should be concave up.
"""
function concave_up_intervals_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}},
                                    concave_up_intervals::AbstractMatrix{T}) where {T <: Real}
    for i in 1:size(concave_up_intervals, 1)
        push!(curvature_settings, (concave_up = true, range = sort(concave_up_intervals[i, :])))
    end
end

"""
```
function concave_down_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Tuple{T, T}}}}) where {T <: Real}
```

adds the properties of an everywhere concave down spline.
"""
function concave_down_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}}) where {T <: Real}
    push!(curvature_settings, (concave_up = false, range = Vector{T}(undef, 0)))
end

"""
```
function concave_down_intervals_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}},
                                      concave_up_intervals::AbstractMatrix{T}) where {T <: Real}
```

adds information about interval(s) over which the spline should be concave up.
"""
function concave_down_intervals_info!(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}},
                                      concave_down_intervals::AbstractMatrix{T}) where {T <: Real}
    for i in 1:size(concave_down_intervals, 1)
        push!(curvature_settings, (concave_up = false, range = sort(concave_down_intervals[i, :])))
    end
end

"""
```
function construct_curvature_matrix(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}},
                                    nc::Int, nk::Int, knots::AbstractVector{T}, dknots::AbstractVector{T},
                                    Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
                                    use_sparse::Bool = false) where {T <: Real}
```

creates the inequalities enforcing curvature for the spline
"""
function construct_curvature_matrix(curvature_settings::Vector{NamedTuple{(:concave_up, :range), Tuple{Bool, Vector{T}}}},
                                    nc::Int, nk::Int, knots::AbstractVector{T}, dknots::AbstractVector{T},
                                    Mineq::AbstractMatrix{T}, rhsineq::AbstractVector{T};
                                    use_sparse::Bool = false) where {T <: Real}

    L = length(curvature_settings)
    if L > 1
        # Then not entirely concave up or down so need to handle construction of M
        dim1 = nk - 1
        for i in 1:L
            end_range = map(a -> max(min(a, knots[end]), knots[1]), curvature_settings[i][:range])
            dim1 += length(bin_sort(end_range, knots))
        end
        M = use_sparse ? spzeros(T, dim1, nc) : zeros(T, dim1, nc)
    else
        M = use_sparse ? spzeros(T, nk, nc) : zeros(T, nk, nc)
    end
    n = 0

    for i in 1:L
        if isempty(curvature_settings[i][:range])
            # Enter domain specified to be curved in one direction
            if curvature_settings[i][:concave_up]
                for j in 1:(nk - 1)
                    n += 1
                    M[n, j]     = 6. / dknots[j] ^ 2
                    M[n, j + 1] = -6. / dknots[j] ^ 2

                    M[n, nk + j]     = 4. / dknots[j]
                    M[n, nk + j + 1] = 2. / dknots[j]
                end
            else
                for j in 1:(nk - 1)
                    n += 1
                    M[n, j]     = -6. / dknots[j] ^ 2
                    M[n, j + 1] = 6. / dknots[j] ^ 2

                    M[n, nk + j]     = -4. / dknots[j]
                    M[n, nk + j + 1] = -2. / dknots[j]
                end
            end

            n += 1
            if curvature_settings[i][:concave_up]
                M[n, nk - 1] = -6 / dknots[end] ^ 2
                M[n, nk]     = 6 / dknots[end] ^ 2

                M[n, 2 * nk - 1] = -2 / dknots[end]
                M[n, 2 * nk]     = -4 / dknots[end]
            else
                M[n, nk - 1] = 6 / dknots[end] ^ 2
                M[n, nk]     = -6 / dknots[end] ^ 2

                M[n, 2 * nk - 1] = 2 / dknots[end]
                M[n, 2 * nk]     = 4 / dknots[end]
            end
        else
            # Only enforce curvature between the given range limits, knot by knot first
            # and then at the endpoints of the range
            curv_setting = curvature_settings[i]
            if curv_setting[:concave_up]
                for j in 1:(nk - 1)
                    if (knots[j] < curv_setting[:range][2]) &&
                        (knots[j] >= curv_setting[:range][1])

                        n += 1
                        M[n, j]     = 6. / dknots[j] ^ 2
                        M[n, j + 1] = -6. / dknots[j] ^ 2

                        M[n, nk + j]     = 4. / dknots[j]
                        M[n, nk + j + 1] = 2. / dknots[j]
                    end
                end
            else
                for j in 1:(nk - 1)
                    if (knots[j] < curv_setting[:range][2]) &&
                        (knots[j] >= curv_setting[:range][1])

                        n += 1
                        M[n, j]     = -6. / dknots[j] ^ 2
                        M[n, j + 1] = 6. / dknots[j] ^ 2

                        M[n, nk + j]     = -4. / dknots[j]
                        M[n, nk + j + 1] = -2. / dknots[j]
                    end
                end
            end

            end_range = map(a -> max(min(a, knots[end]), knots[1]), curv_setting[:range])
            ind = bin_sort(end_range, knots)
            t = (end_range - knots[ind]) ./ dknots[ind]
            s = 1. .- t

            if curv_setting[:concave_up]
                for j in 1:length(ind)
                    M[n + j, ind[j]]          = -(6. - 12. * s[j]) / dknots[ind[j]] ^ 2
                    M[n + j, ind[j] + 1]      = -(6. - 12. * t[j]) / dknots[ind[j]] ^ 2
                    M[n + j, ind[j] + nk]     = -(2. - 6. * s[j]) / dknots[ind[j]]
                    M[n + j, ind[j] + nk + 1] = -(6. * t[j] - 2.) / dknots[ind[j]]
                end
            else
                for j in 1:length(ind)
                    M[n + j, ind[j]]          = (6. - 12. * s[j]) / dknots[ind[j]] ^ 2
                    M[n + j, ind[j] + 1]      = (6. - 12. * t[j]) / dknots[ind[j]] ^ 2
                    M[n + j, ind[j] + nk]     = (2. - 6. * s[j]) / dknots[ind[j]]
                    M[n + j, ind[j] + nk + 1] = (6. * t[j] - 2.) / dknots[ind[j]]
                end
            end

            n += length(ind)
        end
    end

    return  vcat(Mineq, M), vcat(rhsineq, zeros(T, size(M, 1)))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# test default packages

cd(joinpath(JULIA_HOME,"..","share","julia","site","v$(VERSION.major).$(VERSION.minor)")) do
    Pkg.Entry.test(convert(Vector{AbstractString}, readdir()))
end

max_stars_repo_name: tpapp/julia

max_issues_repo_name: anowacki/julia

max_forks_repo_name: anowacki/julia


text:
using Test

println("\n=== Starting tests ===")
println("NOTE: SPD generation failures may cascade into image solve.\n")

@time include("spr_riav_comparison.jl")
@time include("imagesolve.jl")


text:
# # Future climate data

using SimpleSDMLayers
using Plots
using Statistics

# **Justification for this use case:** we will often want to forecast the range
# of species under a variety of climate change scenarios. For this reason,
# `SimpleSDMLayers` offers access to CMIP5 and CMIP6 scenarios or models, to be
# used in these situations.

# For some data providers and datasets, `SimpleSDMLayers` offers access to future
# climate data. Future climates are usually specified by a model, and a
# scenario. For example, WorldClim 2.1 offers the full suite of BioClim variable
# under four SSPs and a number of CMIP6 models.

# We can use this to look at, for example, the temperature difference between the
# current and future climate. To illustrate this, we will do a simple example
# where we contrast the "historical" climate (*i.e.* what is assumed to be the
# current data) to the projected  data under SSP585 in the 2041-2060 period.

# We will start by getting the contemporary data:

boundaries = (left=-169.0, right=-50.0, bottom=24.0, top=71.0)
baseline = SimpleSDMPredictor(WorldClim, BioClim, 1; boundaries...)
contour(baseline; c=:cork, frame=:box, fill=true, clim=(-30, 30), levels=6)

# To get a future dataset, we need to specify the model - models are
# combinations of a CMIP version, and either a RCP or SSP:

instances(CMIP6)

# And we need to check the names of the SSP we want to use:

instances(SharedSocioeconomicPathway)

# We can now get our future temperature layer (and plot it), for a valid
# model/scenario combination, and plot it. We will pick an extreme scenario at a
# long (yet frighteningly short) time in the future. In order to facilitate the
# comparison with the previous plot, we will re-use the same color limites,
# where blue is negatuve temperatures.

future = SimpleSDMPredictor(WorldClim, BioClim, CanESM5, SSP585, 1; year="2061-2080", boundaries...)
contour(future; c=:cork, frame=:box, fill=true, clim=(-30, 30), levels=6)

# Note that the call to get the future data is almost the same as the historical
# one - the exception is the addition of the model and scenario, and the
# specification of the years. With this layer, we can now measure the difference
# in mean annual temperature, because layers can be substracted. Note that we
# are switching the scale: the difference between the two layers here is always
# positive.

plot(future - baseline, frame=:box, c=:lajolla)

# We might actually be interested in averaging multiple models. Because we know
# the variety of models worldclim has (`instances(CMIP6)`), we can do this
# fairly easily. One of the model has no predictions for SSP585 (which we would
# learn in the form of an error message), so we will filter it out directly.

ensemble = [
    SimpleSDMPredictor(
        WorldClim, BioClim, model, SSP585, 1;
        year="2061-2080", boundaries...
    ) for model in instances(CMIP6) if model != GFDLESM4
];

# We can create a layer of differences from each scenario to the baseline:

differences = [component - baseline for component in ensemble];

# This can be plotted as a grid of differences, using the same colorscale as in
# the previous plot:

plot(plot.(differences, c=:lajolla, grid=:none, axes=false, frame=:none, leg=false)...)

# Finally, we can plot the average of the expected conditions under the scenario we considered:

contour(mosaic(mean, ensemble); c=:cork, frame=:box, fill=true, clim=(-30, 30), levels=6)


text:
module SynapseClientTests

using SynapseClient, PyCall
using Test
using Dates
using UUIDs
import SynapseClient: utils,
                      AbstractEntity, AbstractSynapse, AbstractSynapseDict,
                      Activity, Folder, File, Project, Evaluation, Submission, DictObject, Team

ValueError = pybuiltin(:ValueError)
PyKeyError = pybuiltin(:KeyError)

_raise_for_status = SynapseClient.synapseclient.exceptions._raise_for_status
SynapseMalformedEntityError = SynapseClient.synapseclient.exceptions.SynapseMalformedEntityError
SynapseHTTPError = SynapseClient.synapseclient.exceptions.SynapseHTTPError

PyDictObject = SynapseClient.synapseclient.dict_object.DictObject
PyFile = SynapseClient.synapseclient.File

macro catchpyerror(expr)
	esc(quote
		begin
			err = nothing
			try
				$expr
			catch e
				typeof(e) == PyCall.PyError || rethrow(e)
				err=e.T
			end
			err
		end
	end)
end
macro test_pythrows( exception, expr )
	esc(:(@test @catchpyerror($expr) == $exception))
end

include("unit/run_unit_tests.jl")
include("integration/run_integration_tests.jl")

end


text:
for T in (:CLCL, :OPOP)
   @eval begin
       function (+){Q<:Qualia,C<:Clay}(a::Flex{$T,Q,C}, b::Flex{$T,Q,C})
           lo = (+)(a.lo, b.lo, RoundDown)
           hi = (+)(a.hi, b.hi, RoundUp)
           Flex{$T,Q,C}(lo,hi)
       end
   end
end

function (+){S<:Sculpt, T<:Sculpt, Q<:Qualia, C<:Clay}(a::Flex{S,Q,C}, b::Flex{T,Q,C})
    aLoIsOpen, aHiIsOpen = boundries(S)
    bLoIsOpen, bHiIsOpen = boundries(T)
    sculpting = boundries(aLoIsOpen|bLoIsOpen, aHiIsOpen|bHiIsOpen)

    lo = (+)(a.lo, b.lo, RoundDown)
    hi = (+)(a.hi, b.hi, RoundUp)

    Flex{sculpting, Q, C}(lo,hi)
end

(+){S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}, b::C) = (+)(a, Flex{S,Q,C}(b))
(+){S<:Sculpt,Q<:Qualia,C<:Clay}(a::C, b::Flex{S,Q,C}) = (+)(Flex{S,Q,C}(a), b)

(+){S<:Sculpt,Q<:Qualia,C<:Clay,T<:Real}(a::Flex{S,Q,C}, b::T) = (+)(a, Flex{S,Q,C}(convert(C,b)))
(+){S<:Sculpt,Q<:Qualia,C<:Clay,T<:Real}(a::T, b::Flex{S,Q,C}) = (+)(Flex{S,Q,C}(convert(C,a)), b)

(+){S<:Sculpt,T<:Sculpt,Q<:Qualia,R<:Qualia,C<:Clay,D<:Clay}(a::Flex{S,Q,C}, b::Flex{T,R,D}) = (+)(promote(a,b)...)


text:
"""
    SineOperation(A::AbstractTensor{rank}) where rank

Represents computing sin of A.
"""
struct SineOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor}
    freeindices::Tuple{}
    function SineOperation(A::AbstractTensor{rank}) where rank
        new{rank}(A.shape, (A,), ())
    end
end

"""
    Base.sin(s::AbstractTensor{0})

Returns a [`SineOperation`](@ref) represeting sine of s in radians.
"""
function Base.sin(s::AbstractTensor{0})
    return SineOperation(Tensor(s))
end

"""
    CosineOperation(A::AbstractTensor{rank}) where rank

Represents computing cos of A.
"""
struct CosineOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor}
    freeindices::Tuple{}
    function CosineOperation(A::AbstractTensor{rank}) where rank
        new{rank}(A.shape, (A,), ())
    end
end

"""
    Base.cos(s::AbstractTensor{0})

Returns a [`CosineOperation`](@ref) represeting cosine of s in radians.
"""
function Base.cos(s::AbstractTensor{0})
    return CosineOperation(Tensor(s))
end

"""
    TangentOperation(A::AbstractTensor{rank}) where rank

Represents computing tan of A.
"""
struct TangentOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor}
    freeindices::Tuple{}
    function TangentOperation(A::AbstractTensor{rank}) where rank
        new{rank}(A.shape, (A,),())
    end
end

"""
    Base.tan(s::AbstractTensor{0})

Returns a [`TangentOperation`](@ref) represeting tangent of s in radians.
"""
function Base.tan(s::AbstractTensor{0})
    return TangentOperation(Tensor(s))
end

"""
    ArcsineOperation(A::AbstractTensor{rank}) where rank

Represents computing arcsin of A.
"""
struct ArcsineOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor}
    freeindices::Tuple{}
    function ArcsineOperation(A::AbstractTensor{rank}) where rank
        if A isa ConstantTensor{T} where T<: Number && abs(A.value) > 1
            throw(DomainError(A, " asin is not defined for |x| > 1."))
        end
        new{rank}(A.shape, (A,),())
    end
end

"""
    Base.asin(s::AbstractTensor{0})

Returns a [`ArcsineOperation`](@ref) represeting arcsine of s.
"""
function Base.asin(s::AbstractTensor{0})
    return ArcsineOperation(Tensor(s))
end

"""
    CosineOperation(A::AbstractTensor{rank}) where rank

Represents computing cos of A.
"""
struct ArccosineOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor}
    freeindices::Tuple{}
    function ArccosineOperation(A::AbstractTensor{rank}) where rank
        if A isa ConstantTensor{T} where T<: Number && abs(A.value) > 1
            throw(DomainError(A, " acos is not defined for |x| > 1."))
        end
        new{rank}(A.shape, (A,),())
    end
end

"""
    Base.acos(s::AbstractTensor{0})

Returns a [`ArccosineOperation`](@ref) represeting arccosine of s.
"""
function Base.acos(s::AbstractTensor{0})
    return ArccosineOperation(Tensor(s))
end

"""
    ArctangentOperation(A::AbstractTensor{rank}) where rank

Represents computing arctan of A.
"""
struct ArctangentOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor}
    freeindices::Tuple{}
    function ArctangentOperation(A::AbstractTensor{rank}) where rank
        new{rank}(A.shape, (A,),())
    end
end

"""
    Base.atan(s::AbstractTensor{0})

Returns a [`ArctangentOperation`](@ref) represeting arctangent of s.
"""
function Base.atan(s::AbstractTensor{0})
    return ArctangentOperation(Tensor(s))
end

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
import Base: getindex, setindex!

immutable seq
    start::Cdouble
    stop::Cdouble
    step::Cdouble
end

function create_seq_object(r::Range)
    start = r.start - 1
    stop = r.stop - 1
    if :step in fieldnames(r)
        step = r.step
    else
        step = 1
    end
    s = seq(start, stop, step)
end

function create_seq_object(i::Int)
    start = i - 1
    stop = i - 1
    step = 1
    s = seq(start, stop, step)
end

function create_seq_object(::Colon)
    start = 1
    stop = 1
    step = 0
    s = seq(start, stop, step)
end

immutable index
    ptr::Ptr{Void}
end

function getindex{T}(a::AFArray{T}, idx::Union{Range,Int,Colon,AFArray}...)
    indexers = new_ptr()
    af_create_indexers(indexers)
    indexers = index(indexers[])
    set_up_index!(idx, indexers)
    out = new_ptr()
    n = ndims(a)
    l = length(idx)
    l <= n || throw(DimensionMismatch("Number of dimensions is lesser than number of indices"))
    af_index_gen(out, a, l, indexers)
    af_release_indexers(indexers)
    outarr = AFArray{T}(out[])
    if length(outarr) == 1
        return Array(outarr)[1]
    else
        return outarr
    end
end

function set_up_index!(idx::Tuple, indexers::index)
    for (i,thing) in enumerate(idx)
        t = typeof(thing)
        if t <: Range || t <: Int || t <: Colon
            s = create_seq_object(thing)
            af_set_seq_indexer(indexers, s, i-1, true)
        elseif t <: AFArray
            _shift = thing - 1
            af_set_array_indexer(indexers, _shift, i-1)
        end
    end
end

function setindex!{T,S}(lhs::AFArray{T}, rhs::AFArray{S}, idx::Union{Range,Int,Colon,AFArray}...)
    indexers = new_ptr()
    af_create_indexers(indexers)
    indexers = index(indexers[])
    set_up_index!(idx, indexers)
    n = ndims(lhs)
    l = length(idx)
    rhs = convert(AFArray{T}, rhs)
    l <= n || throw(DimensionMismatch("Number of dimensions is lesser than number of indices"))
    af_assign_gen(pointer_from_objref(lhs), lhs, l, indexers, rhs)
    af_release_indexers(indexers)
    rhs
end

function setindex!{T,S<:Real}(lhs::AFArray{T}, val::S, idx::Union{Range,Int,Colon,AFArray}...)
    sz = get_sizes(idx, lhs)
    rhs = constant(T(val), sz...)
    setindex!(lhs, rhs, idx...)
    val
end

function get_sizes(idx::Tuple, lhs::AFArray)
    s = Array(Int, length(idx))
    for i = 1:length(idx)
        if typeof(idx[i]) <: Range
            s[i] = length(idx[i])
        elseif idx[i] == Colon()
            s[i] = size(lhs,i)
        elseif typeof(idx[i]) <: Integer
            s[i] = 1
        elseif typeof(idx[i]) <: AFArray
            s[i] = length(idx[i])
        end
    end
    tuple(s...)
end

Base.LinearIndexing(::Type{AFArray}) = Base.LinearSlow()


text:
module Test_bouncingBall_with_measurement_analytic

using Measurements

const g  = 9.81
const e  = 0.7 ± 0.1
const h0 = 1.0 ± 0.2

#=
der(h) = v  ; h(0) = h0
der(v) = -g ; v(0) = 0

when h <= 0 then
	v = -e*v
end when

Analytic solution of differential equation with
h(t_ev) = h_ev, v(t_ev) = v_ev:

	v = -g*(t-t_ev) + v_ev
	h = -g/2*(t-t_ev)^2 + v_ev*(t-t_ev) + h_ev

First phase (t_ev = 0, v_ev=0, h_ev=h0)
	v = -g*t
	h = -g/2*t^2 + h0

First event:
	hev1 = 0 -> t_ev1 = sqrt(2*h0/g)
	            v_ev1_before = -g*t_ev1
				v_ev1_after  = -e*v_ev1_before

Second phase (t_ev = t_ev1, v_ev = v_ev1_after, h_ev = 0)
    v = -g*(t-t_ev1) + v_ev1_after
    h = -g/2*(t-t_ev1)^2 + v_ev1_after*(t-t_ev1)

Second event:
	h_ev2 = 0 = -g/2*(t-t_ev1)^2 + vev1_after*(t-tev_1)
	          = (t-t_ev1)*(-g/2*(t-t_ev1) + v_ev1_after)

	-> t_ev2-t_ev1 = 2*v_ev1_after/g

	v_ev2_before = -g*(t_ev2-t_ev1) + v_ev1_after
	             = -g*2*v_ev1_after/g + v_ev1_after
	   			 = -v_ev1_after
	v_ev2_after  = -e*v_ev2_before
=#

t_ev1        = sqrt(2*h0/g)
h_ev1        = -g/2*t_ev1^2 + h0
v_ev1_before = -g*t_ev1
v_ev1_after  = -e*v_ev1_before

t_ev2        = t_ev1 +2*v_ev1_after/g
h_ev2        = -g/2*(t_ev2-t_ev1)^2 + v_ev1_after*(t_ev2-t_ev1)
v_ev2_before = -v_ev1_after
v_ev2_after  = -e*v_ev2_before

@show t_ev1
@show h_ev1
@show v_ev1_before
@show v_ev1_after
@show t_ev2
@show h_ev2
@show v_ev2_before
@show v_ev2_after

# Alternative that computes t_ev1 with Roots.jl:

println("\n using Roots:")

using Roots
f(t)  = -g/2*t^2 + h0
t_ev1 = find_zero(f,(0.1±0.0, 0.5±0.0))
h_ev1 = f(t_ev1)
@show t_ev1, h_ev1
end


text:
"""
Flac stream encoder object.

All it contains is an opaque pointer.
The zero-argument constructor creates a new stream decoder and wraps the pointer.

The type is not an immutable because it uses a finalizer.
"""
type StreamEncoderPtr  # type not immutable so that finalizer can be applied
    v::Ptr{Void}
end

"""
Allows for passing the instance of the type in a `ccall`.
"""
Base.unsafe_convert(::Type{Ptr{Void}},en::StreamEncoderPtr) = en.v

function StreamEncoderPtr()
    en = StreamEncoderPtr(ccall((:FLAC__stream_encoder_new,libflac),Ptr{Void},()))
    finalizer(en,x->ccall((:FLAC__stream_encoder_delete,libflac),Void,(Ptr{Void},),x.v))
    en
end

for (nm,typ) in (("verify",:Bool),
                 ("compression_level",:Cuint), # a value between 0 (fastest, least compression) to 8 (most compression) - default 3
                 ("channels",:Cuint),          # default 2
                 ("bits_per_sample",:Cuint),   # default 16
                 ("sample_rate",:Cuint),       # default 44100
                 ("total_samples_estimate",:Cuint), # default 0
                 ("blocksize",:Cuint))         # default 0 - encoder estimates
    @eval begin
        function $(symbol(string("set_",nm)))(en::StreamEncoderPtr,val)
            ccall(($(string("FLAC__stream_encoder_set_",nm)),libflac),Bool,(Ptr{Void},$typ),en,val)
        end
        function $(symbol(string("get_",nm)))(en::StreamEncoderPtr)
            ccall(($(string("FLAC__stream_encoder_get_",nm)),libflac),$typ,(Ptr{Void},),en)
        end
    end
end

function set_metadata(en::StreamEncoderPtr,mdpv::Vector{Ptr{StreamMetaData}})
    ccall((:FLAC__stream_encoder_set_metadata,libflac),Cint,
          (Ptr{Void},Ptr{Ptr{StreamMetaData}},Cuint),en,mdpv,length(mpdv)) ||
        error("Call to stream_encoder_set_metadata failed")
    en
end

@enum(StreamEncoderState,
      EncoderOK,
      EncoderUninitialized,
      EncoderOggError,
      EncoderVerifyDecoderError,
      EncoderVerifyMismatchInAudioData,
      EncoderClientError,
      EncoderIOError,
      EncoderFramingError,
      EncoderMemoryAllocationError)

get_state(en::StreamEncoderPtr) =
    ccall((:FLAC__stream_encoder_get_state,libflac),StreamEncoderState,(Ptr{Void},),en)

@enum(StreamEncoderInitStatus,
      EncoderInitOK,
      EncoderInitEncoderError,
      EncoderInitUnsupportedContainer,
      EncoderInitInvalidCallbacks,
      EncoderInitInvalidNumberOfChannels,
      EncoderInitInvalidBitsPerSample,
      EncoderInitInvalidSampleRate,
      EncoderInitInvalidBlockSize,
      EncoderInitInvalidMaxLPCSize,
      EncoderInitInvalidQLPCoeffPrec,
      EncoderInitBlockSizeTooSmallForLPCOrder,
      EncoderInitNotStreamable,
      EncoderInitInvalidMetadata,
      EncoderInitAlreadyInitialized)

@enum(StreamEncoderReadStatus,
      ReadStatusContinue,
      ReadStatusEndOfStream,
      ReadStatusAbort,
      ReadStatusUnsupported)

@enum(StreamEncoderWriteStatus,
      WriteStatusOK,
      WriteStatusFatalError)

@enum(StreamEncoderSeekStatus,
      SeekStatusOK,
      SeekStatusError,
      SeekStatusUnsupported)

@enum(StreamEncoderTellStatus,
      TellStatusOK,
      TellStatusError,
      TellStatusUnsupported)

finish(en::StreamEncoderPtr) = ccall((:FLAC__stream_encoder_finish,libflac), Int32, (Ptr{Void},), en)

"standard progress callback function"
function pcallback(en::Ptr{Void}, bytesw::Int64,samplesw::Int64,
                   framesw::Cint, totalframesest::Cint, client::Ptr{Void})
    nothing
end

"""
### init_file

Initialize the `StreamEncoder` object `en` to write the file `fnm`.

Note that setting stream characteristics (`channels`, `bits_per_sample`, etc.)
must be done **before** initializing the encoder.
"""
function initfile!(en::StreamEncoderPtr, fnm::ByteString)
    ec = ccall((:FLAC__stream_encoder_init_file, libflac), StreamEncoderInitStatus,
               (Ptr{Void}, Ptr{UInt8}, Ptr{Void}, Ptr{Void}),
               en, fnm, pcallback_c, C_NULL)
    ec == EncoderInitOK || error("Error code $ec from stream_encoder_init_file")
    en
end

function process_interleaved(en::StreamEncoderPtr, buf::Vector{Int32})
    nsamp = div(length(buf), get_channels(en))
    ccall((:FLAC__stream_encoder_process_interleaved, libflac), Bool,
          (Ptr{Void}, Ptr{Int32}, UInt32), en, buf, nsamp) ||
         error("process_interleaved failed: encoder_state is $(get_state(en))")
    nothing
end

function Base.show(io::IO,en::StreamEncoderPtr)
    println("FLAC Stream Encoder")
    println(" channels: ", Int(get_channels(en)))
    println(" bits per sample: ", Int(get_bits_per_sample(en)))
    println(" sample rate: ", Int(get_sample_rate(en)))
    println(" blocksize: ", Int(get_blocksize(en)))
    println(" encoder state: ", get_state(en))
    println()
end

# Cheat save() implementation for people who want to input a 1d array instead of the proper 2d array
save{T<:Real}(f::File{format"FLAC"}, data::Array{T,1}, samplerate; kwargs...) = save(f, data'', samplerate; kwargs...)

function save{T<:Real}(f::File{format"FLAC"}, data::Array{T,2}, samplerate; bits_per_sample = 24, compression_level = 3)
    encoder = StreamEncoderPtr()

    # Set encoder parameters
    set_compression_level(encoder, compression_level)
    num_samples = prod(size(data))
    set_total_samples_estimate(encoder, num_samples)
    set_channels(encoder, size(data)[2])
    set_sample_rate(encoder, samplerate)
    set_bits_per_sample(encoder, bits_per_sample)

    # Open file, make sure encoder was properly initialized
    initfile!(encoder, f.filename)
    if get_state(encoder) != EncoderOK
        throw(InvalidStateException("Encoder state after init_file is $(get_state(en))"))
    end

    # Shove interleaved samples into the encoder by transposing, and convert to Int32
    data_t = round(Int32, data'*2^(bits_per_sample - 1))
    blocksize = get_blocksize(encoder)
    for idx in 1:div(num_samples, blocksize)
        block_idxs = ((idx - 1) * blocksize + 1):(idx * blocksize)
        process_interleaved(encoder, data_t[block_idxs])
    end
    process_interleaved(encoder, data_t[end - rem(num_samples, blocksize) + 1:end])
    finish(encoder)
    return nothing
end

# Calculate cfunction versions of all our callbacks once, at runtime, as is necessary with cfunction's
pcallback_c = Ptr{Void}(C_NULL)
function init_encoder_cfunctions()
    global pcallback_c = cfunction(pcallback, Void, (Ptr{Void}, Int64, Int64, Cint, Cint, Ptr{Void}))
end


text:
# TIPS:
# 1
const DEFAULT = 42

global x = 42
function f(args)
	#function body
end
y = f(x::Int + 1)

# 3
function myFunc(a::T, c::Int) where T
  # code
end

# 9
function with_keyword(x; name::String = "Smith")
   # ...
end

# 16
mutable struct Person
	name::String
	height::Float64
	weight::Float64
end
# instead of:
# type Person
# 	name
# 	height
# 	weight
# end

# 18
# Use
file = open("file2.txt", "w")
a = "Julia"
b = "Rust"
write(file, a, " ", b)
# instead of:
write(file, "$a $b")

# TOOLS:
# @time

# linter:
# add Lint
# using Lint
# lintfile("performance.jl")


text:
function part1(input)
    state = Int8[]
    risk = 0
    first_row = true
    i = 0
    while !eof(input)
        c = read(input, UInt8)
        if c == UInt8('\n')
            first_row = false
            i = 0
            continue
        end
        i += 1
        height = Int8(c - UInt8('0'))
        if first_row
            push!(state, height)
            if i > 1
                if abs(state[i - 1]) >= height
                    state[i - 1] = -abs(state[i - 1])
                elseif abs(state[i - 1]) <= height
                    state[i] = -height
                end
            end
        else
            if abs(state[i]) < height
                if state[i] >= 0
                    risk += 1 + state[i]
                end
                state[i] = -height
            elseif abs(state[i]) == height
                state[i] = -height
            else
                state[i] = height
            end

            if i > 1
                if abs(state[i - 1]) >= height
                    state[i - 1] = -abs(state[i - 1])
                elseif abs(state[i - 1]) <= height
                    state[i] = -height
                end
            end
        end
    end
    for height in state
        if height >= 0
            risk += 1 + height
        end
    end
    return risk
end

mutable struct State
    next::Int
    prev::Int
    size::Int
    ref::Int
end

function part2(input)
    state = State[]
    largest_basins = (0, 0, 0)
    first_row = true
    i = 0
    while !eof(input)
        c = read(input, UInt8)
        if c == UInt8('\n')
            first_row = false
            i = 0
            continue
        end
        i += 1
        height = c - UInt8('0')
        if first_row
            if height == 9
                push!(state, State(0, 0, 0, 0))
            else
                if i == 1 || state[i - 1].ref == 0
                    push!(state, State(i, i, 1, i))
                else
                    push!(state, State(state[i - 1].next, i - 1, 1, state[i - 1].ref))
                    state[i - 1].next = i
                    state[state[i].next].prev = i
                end
            end
        else
            if height == 9
                if state[i].next == i
                    largest_basins = update_largest_basins(largest_basins, state[i].size)
                elseif state[i].ref != 0
                    n = state[i].next
                    p = state[i].prev
                    state[n].prev = state[i].prev
                    state[p].next = state[i].next
                    state[n].size += state[i].size
                end
                state[i].next = 0
                state[i].prev = 0
                state[i].size = 0
                state[i].ref = 0
            else
                if state[i].ref != 0
                    state[i].size += 1
                else
                    state[i].next = i
                    state[i].prev = i
                    state[i].size = 1
                    state[i].ref = i
                end
                if i > 1 && state[i - 1].ref != 0
                    merge!(state, i - 1, i)
                end
            end
        end
    end
    for i = 1:length(state)
        if state[i].next == i
            largest_basins = update_largest_basins(largest_basins, state[i].size)
        elseif state[i].ref != 0
            n = state[i].next
            p = state[i].prev
            state[n].prev = state[i].prev
            state[p].next = state[i].next
            state[n].size += state[i].size
        end
    end
    return prod(largest_basins)
end

function merge!(state, i, j)
    state[i].ref == state[j].ref && return
    k = state[j].next
    while true
        state[k].ref = state[i].ref
        k == j && break
        k = state[k].next
    end
    in = state[i].next
    jn = state[j].next
    state[i].next = jn
    state[j].next = in
    state[in].prev = j
    state[jn].prev = i
end

function update_largest_basins(b, n)
    n > b[1] && return (n, b[1], b[2])
    n > b[2] && return (b[1], n, b[2])
    n > b[3] && return (b[1], b[2], n)
    return b
end


text:
export generate, ShipNoise, Pinger

function generate(signal::AbstractArray, nsamples, fs)
  length(signal) >= nsamples && return signal[1:nsamples]
  padded(signal, (0,nsamples-signal))
end

Base.@kwdef struct Pinger
  spl = db2amp(180.0)
  frequency = 38000.0
  duration = 0.02
  start = 0.0
  interval = 1.0
  window = nothing
end

function generate(pinger::Pinger, nsamples, fs)
  x = zeros(nsamples)
  ping = cw(pinger.frequency, pinger.duration; fs=fs, window=pinger.window)
  t = pinger.start
  while t < nsamples/fs
    ndx = round(Int, t*fs)
    n = length(ping)
    ndx+n-1 > nsamples && n = 1+nsamples-ndx
    x[ndx:ndx+n-1] .= ping[1:n]
    t += pinger.interval
  end
  pinger.spl * x
end

Base.@kwdef struct ShipNoise
  spl = db2amp(172.0)
  spec = [(50.0, 1.0)]
end

function generate(ship::ShipNoise, nsamples, fs)
  x = rand(PinkGaussianNoise(), nsamples)
  y = copy(x)
  for s in ship.spec
    y += s[2]*cw(s[1], nsamples/fs; fs=fs) .* x
  end
  lpf = digitalfilter(Lowpass(3000.0, fs=fs), FIRWindow(hanning(15)))
  y = filtfilt(lpf, y)
  √ship.spl * y/rms(y)
end


text:
using ExperimentsDB
using Test

using Tables
using DataStructures

# For testing `datamerge`
struct Merger
    val::Int
end
ExperimentsDB.datamerge(a::Merger, b::Merger) = Merger(a.val + b.val)

function keycheck(db::ExperimentsDB.Database)
    ks = keys(db)
    for entry in db
        @test collect(keys(entry)) == ks
    end
end

@testset "ExperimentsDB.jl" begin
    ### todict
    x = ExperimentsDB.todict((a = 1, b = 2))
    @test isa(x, DataStructures.SortedDict{Symbol, Any})
    @test length(x) == 2
    @test x[:a] == 1
    @test x[:b] == 2

    # Test pass-through logic
    @test ExperimentsDB.todict(x) === x

    # Pass in a normal dict
    y = Dict("a" => 1, "b" => 2)
    x = ExperimentsDB.todict(y)
    @test isa(x, DataStructures.SortedDict{Symbol, Any})
    @test length(x) == 2
    @test x[:a] == 1
    @test x[:b] == 2

    y = Dict(:a => 1, :b => 2)
    x = ExperimentsDB.todict(y)
    @test isa(x, DataStructures.SortedDict{Symbol, Any})
    @test length(x) == 2
    @test x[:a] == 1
    @test x[:b] == 2

    x = ExperimentsDB.todict(DataStructures.SortedDict{String,Any}, y)
    @test isa(x, DataStructures.SortedDict{String, Any})
    @test length(x) == 2
    @test x["a"] == 1
    @test x["b"] == 2

    ### Entry
    data = Dict(
        "entry_1" => 1,
        "entry_2" => 2,
    )
    A = ExperimentsDB.Entry((a = 1, b = 2), data)
    @test collect(keys(A)) == [:a, :b]
    @test A[] == data
    @test A[:a] == 1
    @test A[:b] == 2

    # Tables API
    @test Tables.columnnames(A) == [:a, :b]
    @test Tables.getcolumn(A, :a) == 1
    @test Tables.getcolumn(A,  1) == 1
    @test Tables.getcolumn(A, :b) == 2
    @test Tables.getcolumn(A,  2) == 2

    params = ExperimentsDB.parameters(A)
    @test length(params) == 2
    @test params[:a] == 1
    @test params[:b] == 2

    # Test ==
    B = ExperimentsDB.Entry(Dict("a" => 1, "b" => 2))
    @test A == B

    @test issubset(ExperimentsDB.todict((a = 1,)), A)
    @test issubset(ExperimentsDB.todict((b = 2,)), A)
    @test issubset(ExperimentsDB.todict((a = 1, b = 2,)), A)
    @test !issubset(ExperimentsDB.todict((a = 2, b = 2,)), A)
    @test !issubset(ExperimentsDB.todict((c = nothing,)), A)
    @test in(ExperimentsDB.todict((a = 1, b = 2,)), A)
    @test !in(ExperimentsDB.todict((a = 1,)), A)

    # expand!
    nt = (a = 1, b = 2)
    x = ExperimentsDB.todict(nt)
    ExperimentsDB.expand!(x, [:a, :c])
    @test x[:a] == 1
    @test x[:b] == 2
    @test x[:c] == nothing

    entry = ExperimentsDB.Entry(nt)
    ExperimentsDB.expand!(entry, [:a, :c])
    @test ExperimentsDB.parameters(entry) == x

    #####
    ##### DataBase
    #####
    db = ExperimentsDB.Database()
    @test length(db) == 0
    @test keys(db) == Symbol[]

    # Show empty database
    show(devnull, db)

    data1 = Dict(
        "merger" => Merger(10),
        "int" => 10,
    )
    data2 = Dict(
        "merger" => Merger(20),
        "int" => 20,
    )
    data3 = Dict(
        "merger" => Merger(30),
        "int" => 30
    )

    db[(a = 1, c = 2)] = data1
    @test length(db) == 1
    @test keys(db) == [:a, :c]

    ### We now have several different cases to test.
    ###
    ### Case 1: Keys of new entry are a subset of the existing keys, but a match is not found.
    ### Case 2: Keys of new entry are a subset of the existing keys, but a match is found.
    ### Case 3: Not a subset

    ## Case 1
    db[(a = 1,)] = data2
    @test length(db) == 2
    @test keys(db) == [:a, :c]

    # The second entry should have had its keys expanded
    entry = db[2]
    @test collect(keys(entry)) == [:a, :c]
    @test entry[:a] == 1
    @test entry[:c] == nothing
    @test db[1][] == data1
    @test db[2][] == data2

    ## Case 2
    db[(a = 1, c = 2)] = data2
    @test length(db) == 2
    keycheck(db)

    # Check if data was merged properly.
    data = db[1][]
    # The integer should have been replaced.
    @test data["int"] == 20
    # The Merger should have been merged
    @test data["merger"] == Merger(10 + 20)

    ## Case 3
    db[(b = 3,)] = data3
    @test length(db) == 3
    @test keys(db) == [:a, :b, :c]
    keycheck(db)

    ### Done testing merging logic
    # Test showing DB
    show(devnull, db)

    # Now, do some testing for "getindex"
    db2 = db[[1,2,]]
    @test isa(db2, ExperimentsDB.Database)
    @test length(db2) == 2
    @test db[1] == db2[1]
    @test db[1][] === db2[1][]
    @test db[2] == db2[2]
    @test db[2][] === db2[2][]

    # Test filtering logic
    db2 = db[a = 1]
    @test isa(db2, ExperimentsDB.Database)
    @test length(db2) == 2
    keycheck(db2)
    @test all(isequal(1), db2[:a])

    db2 = db[abba = nothing]
    @test isa(db2, ExperimentsDB.Database)
    @test length(db2) == 0

    # Finally, test the "in" and "subset" logic
    @test in(ExperimentsDB.todict((a = 1, b = nothing, c = 2)), db)
    @test !in(ExperimentsDB.todict((a = 1, c = 2)), db)
    @test issubset(ExperimentsDB.todict((a = 1, c = 2)), db)
    @test !issubset(ExperimentsDB.todict((a = 1, d = 2)), db)

    # Findonly
    @test isa(ExperimentsDB.findonly(db, (a = 1, b = nothing, c = 2)), ExperimentsDB.Entry)
    @test isnothing(ExperimentsDB.findonly(db, (a = 1, b = nothing)))
end