text:
# The file contains test examples definitions

# nearest TODO: generate answers for these
# and also consider doing this in another way

function rootn(n; ground=QQ)
    R, xs = PolynomialRing(ground, ["x$i" for i in 1:n])
    ans = [
        sum(map(prod, Combinatorics.combinations(xs, i)))
        for i in 1:n
    ]
    ans[end] -= (-1)^(n - 1)
    ans
end

function henrion5(;ground=QQ)

    R, (f1,f2,f3,f4,f5,t) = PolynomialRing(ground, ["f1","f2","f3","f4","f5","t"])
    fs = [
        2*f1*f2*f3*f4*f5-9823275,
        21//5*f1*f2*f4*f5+16//5*f1*f3*f4*f5+9//5*f2*f3*f4*f5+24//5*f1*f2*f3*f5+5*f4*f3*f1*f2-4465125,
        14//5*f4*f5*f1+14//5*f4*f5*f2+8//5*f3*f4*f5+18//5*f1*f2*f5+24//5*f1*f3*f5+18//5*f2*f3*f5+4*f3*f1*f2+6*f1*f2*f4+6*f3*f4*f1+4*f2*f3*f4-441486,
        7//5*f4*f5+12//5*f5*f1+12//5*f5*f2+12//5*f5*f3+3*f1*f2+4*f3*f1+4*f4*f1+3*f2*f3+4*f4*f2+3*f3*f4-15498,
        6//5*f5+2*f4+2*f3+2*f2+2*f1-215,
        f1+2*f2+3*f3+4*f4+5*f5+6*t
    ]
end

function katsura6(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7) = PolynomialRing(ground, ["x$i" for i in 1:7])

    fs = [
        1*x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7-1,
        2*x4*x3+2*x5*x2+2*x6*x1+2*x7*x2-1*x6,
        1*x3^2+2*x4*x2+2*x5*x1+2*x6*x2+2*x7*x3-1*x5,
        2*x3*x2+2*x4*x1+2*x5*x2+2*x6*x3+2*x7*x4-1*x4,
        1*x2^2+2*x3*x1+2*x4*x2+2*x5*x3+2*x6*x4+2*x7*x5-1*x3,
        2*x2*x1+2*x3*x2+2*x4*x3+2*x5*x4+2*x6*x5+2*x7*x6-1*x2,
        1*x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2-1*x1
    ]
    fs
end

function katsura9(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7, x8, x9) = PolynomialRing(ground, ["x$i" for i in 1:9])

    fs = [
        x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7+2*x8+2*x9-1,
        x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2+2*x8^2+2*x9^2-x1,
        2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6+2*x6*x7+2*x7*x8+2*x8*x9-x2,
        x2^2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6+2*x5*x7+2*x6*x8+2*x7*x9-x3,
        2*x2*x3+2*x1*x4+2*x2*x5+2*x3*x6+2*x4*x7+2*x5*x8+2*x6*x9-x4,
        x3^2+2*x2*x4+2*x1*x5+2*x2*x6+2*x3*x7+2*x4*x8+2*x5*x9-x5,
        2*x3*x4+2*x2*x5+2*x1*x6+2*x2*x7+2*x3*x8+2*x4*x9-x6,
        x4^2+2*x3*x5+2*x2*x6+2*x1*x7+2*x2*x8+2*x3*x9-x7,
        2*x4*x5+2*x3*x6+2*x2*x7+2*x1*x8+2*x2*x9-x8
    ]
end

function katsura10(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) = PolynomialRing(ground, ["x$i" for i in 1:10])

    fs = [
        x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7+2*x8+2*x9+2*x10-1,
        x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2+2*x8^2+2*x9^2+2*x10^2-x1,
        2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6+2*x6*x7+2*x7*x8+2*x8*x9+2*x9*x10-x2,
        x2^2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6+2*x5*x7+2*x6*x8+2*x7*x9+2*x8*x10-x3,
        2*x2*x3+2*x1*x4+2*x2*x5+2*x3*x6+2*x4*x7+2*x5*x8+2*x6*x9+2*x7*x10-x4,
        x3^2+2*x2*x4+2*x1*x5+2*x2*x6+2*x3*x7+2*x4*x8+2*x5*x9+2*x6*x10-x5,
        2*x3*x4+2*x2*x5+2*x1*x6+2*x2*x7+2*x3*x8+2*x4*x9+2*x5*x10-x6,
        x4^2+2*x3*x5+2*x2*x6+2*x1*x7+2*x2*x8+2*x3*x9+2*x4*x10-x7,
        2*x4*x5+2*x3*x6+2*x2*x7+2*x1*x8+2*x2*x9+2*x3*x10-x8,
        x5^2+2*x4*x6+2*x3*x7+2*x2*x8+2*x1*x9+2*x2*x10-x9
    ]
end

function katsura11(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11) = PolynomialRing(ground, ["x$i" for i in 1:11])

    fs = [
        x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7+2*x8+2*x9+2*x10+2*x11-1,
        x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2+2*x8^2+2*x9^2+2*x10^2+2*x11^2-x1,
        2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6+2*x6*x7+2*x7*x8+2*x8*x9+2*x9*x10+2*x10*x11-x2,
        x2^2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6+2*x5*x7+2*x6*x8+2*x7*x9+2*x8*x10+2*x9*x11-x3,
        2*x2*x3+2*x1*x4+2*x2*x5+2*x3*x6+2*x4*x7+2*x5*x8+2*x6*x9+2*x7*x10+2*x8*x11-x4,
        x3^2+2*x2*x4+2*x1*x5+2*x2*x6+2*x3*x7+2*x4*x8+2*x5*x9+2*x6*x10+2*x7*x11-x5,
        2*x3*x4+2*x2*x5+2*x1*x6+2*x2*x7+2*x3*x8+2*x4*x9+2*x5*x10+2*x6*x11-x6,
        x4^2+2*x3*x5+2*x2*x6+2*x1*x7+2*x2*x8+2*x3*x9+2*x4*x10+2*x5*x11-x7,
        2*x4*x5+2*x3*x6+2*x2*x7+2*x1*x8+2*x2*x9+2*x3*x10+2*x4*x11-x8,
        x5^2+2*x4*x6+2*x3*x7+2*x2*x8+2*x1*x9+2*x2*x10+2*x3*x11-x9,
        2*x5*x6+2*x4*x7+2*x3*x8+2*x2*x9+2*x1*x10+2*x2*x11-x10
    ]
end

function katsura12(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12) = PolynomialRing(ground, ["x$i" for i in 1:12])

    fs = [
        x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7+2*x8+2*x9+2*x10+2*x11+2*x12-1,
        x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2+2*x8^2+2*x9^2+2*x10^2+2*x11^2+2*x12^2-x1,
        2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6+2*x6*x7+2*x7*x8+2*x8*x9+2*x9*x10+2*x10*x11+2*x11*x12-x2,
        x2^2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6+2*x5*x7+2*x6*x8+2*x7*x9+2*x8*x10+2*x9*x11+2*x10*x12-x3,
        2*x2*x3+2*x1*x4+2*x2*x5+2*x3*x6+2*x4*x7+2*x5*x8+2*x6*x9+2*x7*x10+2*x8*x11+2*x9*x12-x4,
        x3^2+2*x2*x4+2*x1*x5+2*x2*x6+2*x3*x7+2*x4*x8+2*x5*x9+2*x6*x10+2*x7*x11+2*x8*x12-x5,
        2*x3*x4+2*x2*x5+2*x1*x6+2*x2*x7+2*x3*x8+2*x4*x9+2*x5*x10+2*x6*x11+2*x7*x12-x6,
        x4^2+2*x3*x5+2*x2*x6+2*x1*x7+2*x2*x8+2*x3*x9+2*x4*x10+2*x5*x11+2*x6*x12-x7,
        2*x4*x5+2*x3*x6+2*x2*x7+2*x1*x8+2*x2*x9+2*x3*x10+2*x4*x11+2*x5*x12-x8,
        x5^2+2*x4*x6+2*x3*x7+2*x2*x8+2*x1*x9+2*x2*x10+2*x3*x11+2*x4*x12-x9,
        2*x5*x6+2*x4*x7+2*x3*x8+2*x2*x9+2*x1*x10+2*x2*x11+2*x3*x12-x10,
        x6^2+2*x5*x7+2*x4*x8+2*x3*x9+2*x2*x10+2*x1*x11+2*x2*x12-x11
    ]
end

function eco5(;ground=GF(2^31-1))
    R, (x1, x2, x3, x4, x5) = PolynomialRing(ground, ["x$i" for i in 1:5])

    fs = [
    (x1 + x1*x2 + x2*x3 + x3*x4)*x5 - 1,
     (x2 + x1*x3 + x2*x4)*x5 - 2,
             (x3 + x1*x4)*x5 - 3,
                       x4*x5 - 4,
           x1 + x2 + x3 + x4 + 1
    ]
end

function eco7(;ground=GF(2^31-1))
    R, (x1, x2, x3, x4, x5, x6, x7) = PolynomialRing(ground, ["x$i" for i in 1:7])

    fs = [
        (x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6)*x7 - 1,
     (x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6)*x7 - 2,
     (x3 + x1*x4 + x2*x5 + x3*x6)*x7 - 3,
     (x4 + x1*x5 + x2*x6)*x7 - 4,
     (x5 + x1*x6)*x7 - 5,
     x6*x7 - 6,
     x1 + x2 + x3 + x4 + x5 + x6 + 1
    ]
    fs
end

function eco10(;ground=GF(2^31-1))
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9) = PolynomialRing(ground, ["x$i" for i in 1:11])

    fs = [
        x0*x1*x9+x1*x2*x9+x2*x3*x9+x3*x4*x9+x4*x5*x9+x5*x6*x9+x6*x7*x9+x7*x8*x9+x0*x9-1,
        x0*x2*x9+x1*x3*x9+x2*x4*x9+x3*x5*x9+x4*x6*x9+x5*x7*x9+x6*x8*x9+x1*x9-2,
        x0*x3*x9+x1*x4*x9+x2*x5*x9+x3*x6*x9+x4*x7*x9+x5*x8*x9+x2*x9-3,
        x0*x4*x9+x1*x5*x9+x2*x6*x9+x3*x7*x9+x4*x8*x9+x3*x9-4,
        x0*x5*x9+x1*x6*x9+x2*x7*x9+x3*x8*x9+x4*x9-5,
        x0*x6*x9+x1*x7*x9+x2*x8*x9+x5*x9-6,
        x0*x7*x9+x1*x8*x9+x6*x9-7,
        x0*x8*x9+x7*x9-8,
        x8*x9-9,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+1
    ]
    fs
end

function eco11(;ground=GF(2^31-1))
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) = PolynomialRing(ground, ["x$i" for i in 1:11])

    fs = [
        x0*x1*x10+x1*x2*x10+x2*x3*x10+x3*x4*x10+x4*x5*x10+x5*x6*x10+x6*x7*x10+x7*x8*x10+x8*x9*x10+x0*x10-1,
        x0*x2*x10+x1*x3*x10+x2*x4*x10+x3*x5*x10+x4*x6*x10+x5*x7*x10+x6*x8*x10+x7*x9*x10+x1*x10-2,
        x0*x3*x10+x1*x4*x10+x2*x5*x10+x3*x6*x10+x4*x7*x10+x5*x8*x10+x6*x9*x10+x2*x10-3,
        x0*x4*x10+x1*x5*x10+x2*x6*x10+x3*x7*x10+x4*x8*x10+x5*x9*x10+x3*x10-4,
        x0*x5*x10+x1*x6*x10+x2*x7*x10+x3*x8*x10+x4*x9*x10+x4*x10-5,
        x0*x6*x10+x1*x7*x10+x2*x8*x10+x3*x9*x10+x5*x10-6,
        x0*x7*x10+x1*x8*x10+x2*x9*x10+x6*x10-7,
        x0*x8*x10+x1*x9*x10+x7*x10-8,
        x0*x9*x10+x8*x10-9,
        x9*x10-10,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+x9+1
    ]
    fs
end

function eco12(;ground=GF(2^31-1))
    R, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11) = PolynomialRing(ground, ["x$i" for i in 1:11])

    fs = [
        x0*x1*x11+x1*x2*x11+x2*x3*x11+x3*x4*x11+x4*x5*x11+x5*x6*x11+x6*x7*x11+x7*x8*x11+x8*x9*x11+x9*x10*x11+x0*x11-1,
        x0*x2*x11+x1*x3*x11+x2*x4*x11+x3*x5*x11+x4*x6*x11+x5*x7*x11+x6*x8*x11+x7*x9*x11+x8*x10*x11+x1*x11-2,
        x0*x3*x11+x1*x4*x11+x2*x5*x11+x3*x6*x11+x4*x7*x11+x5*x8*x11+x6*x9*x11+x7*x10*x11+x2*x11-3,
        x0*x4*x11+x1*x5*x11+x2*x6*x11+x3*x7*x11+x4*x8*x11+x5*x9*x11+x6*x10*x11+x3*x11-4,
        x0*x5*x11+x1*x6*x11+x2*x7*x11+x3*x8*x11+x4*x9*x11+x5*x10*x11+x4*x11-5,
        x0*x6*x11+x1*x7*x11+x2*x8*x11+x3*x9*x11+x4*x10*x11+x5*x11-6,
        x0*x7*x11+x1*x8*x11+x2*x9*x11+x3*x10*x11+x6*x11-7,
        x0*x8*x11+x1*x9*x11+x2*x10*x11+x7*x11-8,
        x0*x9*x11+x1*x10*x11+x8*x11-9,
        x0*x10*x11+x9*x11-10,
        x10*x11-11,
        x0+x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+1
    ]
    fs
end

function noon3(;ground=QQ)
    R, (x1, x2, x3) = PolynomialRing(ground, ["x$i" for i in 1:3])
    fs = [
    10x1*x2^2 + 10x1*x3^2 - 11x1 + 10,
    10x2*x1^2 + 10x2*x3^2 - 11x2 + 10,
    10x3*x1^2 + 10x3*x2^2 - 11x3 + 10,
    ]
    fs
end

function noon4(; ground=QQ)
    R, (x1, x2, x3, x4) = PolynomialRing(ground, ["x$i" for i in 1:4])

    fs = [
    10*x1^2*x4+10*x2^2*x4+10*x3^2*x4-11*x4+10,
    10*x1^2*x3+10*x2^2*x3+10*x3*x4^2-11*x3+10,
    10*x1*x2^2+10*x1*x3^2+10*x1*x4^2-11*x1+10,
    10*x1^2*x2+10*x2*x3^2+10*x2*x4^2-11*x2+10
    ]
    fs
end

function noon5(;ground=QQ)
    R, (x1, x2, x3, x4, x5) = PolynomialRing(ground, ["x$i" for i in 1:5])

    fs = [
    10*x1^2*x5+10*x2^2*x5+10*x3^2*x5+10*x4^2*x5-11*x5+10,
    10*x1^2*x4+10*x2^2*x4+10*x3^2*x4+10*x4*x5^2-11*x4+10,
    10*x1^2*x3+10*x2^2*x3+10*x3*x4^2+10*x3*x5^2-11*x3+10,
    10*x1*x2^2+10*x1*x3^2+10*x1*x4^2+10*x1*x5^2-11*x1+10,
    10*x1^2*x2+10*x2*x3^2+10*x2*x4^2+10*x2*x5^2-11*x2+10
    ]
    fs
end

function noon6(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6) = PolynomialRing(ground, ["x$i" for i in 1:6])

    fs = [
    10*x1^2*x6+10*x2^2*x6+10*x3^2*x6+10*x4^2*x6+10*x5^2*x6-11*x6+10,
    10*x1^2*x5+10*x2^2*x5+10*x3^2*x5+10*x4^2*x5+10*x5*x6^2-11*x5+10,
    10*x1^2*x4+10*x2^2*x4+10*x3^2*x4+10*x4*x5^2+10*x4*x6^2-11*x4+10,
    10*x1^2*x3+10*x2^2*x3+10*x3*x4^2+10*x3*x5^2+10*x3*x6^2-11*x3+10,
    10*x1*x2^2+10*x1*x3^2+10*x1*x4^2+10*x1*x5^2+10*x1*x6^2-11*x1+10,
    10*x1^2*x2+10*x2*x3^2+10*x2*x4^2+10*x2*x5^2+10*x2*x6^2-11*x2+10
    ]
    fs
end

function noonn(n; ground=QQ)
    without(x, k) = x[1:end .!= k]

    R, xs = PolynomialRing(ground, ["x$i" for i in 1:n])
    fs = zeros(R, n)
    for i in 1:n
        other = without(xs, i)
        fs[i] = xs[i] * (10*sum(other .^ 2) - 11) + 10
    end
    fs
end

function ku10(;ground=QQ)
    R, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) = PolynomialRing(ground, ["x$i" for i in 1:10])

    fs = [
        5*x1*x2+ 5*x1+ 3*x2+ 55,
        7*x2*x3+ 9*x2+ 9*x3+ 19,
        3*x3*x4+ 6*x3+ 5*x4-4,
        6*x4*x5+ 6*x4+ 7*x5+ 118,
        x5*x6+ 3*x5+ 9*x6+ 27,
        6*x6*x7+ 7*x6+x7+ 72,
        9*x7*x8+ 7*x7+x8+ 35,
        4*x8*x9+ 4*x8+ 6*x9+ 16,
        8*x9*x10+ 4*x9+ 3*x10-51,
        3*x1*x10-6*x1+x10+ 5
    ]
    fs
end

function kinema(;ground=QQ)
    R, (z1, z2, z3, z4, z5, z6, z7, z8, z9) = PolynomialRing(ground, ["z$i" for i in 1:9])

    fs = [
    z1^2 + z2^2 + z3^2 - 12*z1 - 68;
    z4^2 + z5^2 + z6^2 - 12*z5 - 68;
    z7^2 + z8^2 + z9^2 - 24*z8 - 12*z9 + 100;
    z1*z4 + z2*z5 + z3*z6 - 6*z1 - 6*z5 - 52;
    z1*z7 + z2*z8 + z3*z9 - 6*z1 - 12*z8 - 6*z9 + 64;
    z4*z7 + z5*z8 + z6*z9 - 6*z5 - 12*z8 - 6*z9 + 32;
    2*z2 + 2*z3 - z4 - z5 - 2*z6 - z7 - z9 + 18;
    z1 + z2 + 2*z3 + 2*z4 + 2*z6 - 2*z7 + z8 - z9 - 38;
    z1 + z3 - 2*z4 + z5 - z6 + 2*z7 - 2*z8 + 8;
    ]
end

function sparse5(; ground=QQ)
    R, (x1, x2, x3, x4, x5) = PolynomialRing(ground, ["x$i" for i in 1:5])

    fs = [
        x1^2*x2^2*x3^2*x4^2*x5^2 + 3*x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + 3*x2^2 + x3^2 + x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + x2^2 + 3*x3^2 + x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + x2^2 + x3^2 + 3*x4^2 + x5^2 + x1*x2*x3*x4*x5 + 5,
        x1^2*x2^2*x3^2*x4^2*x5^2 + x1^2 + x2^2 + x3^2 + x4^2 + 3*x5^2 + x1*x2*x3*x4*x5 + 5
    ]
end

function generate_set(nvariables, exps, nterms, npolys, csz, rng,
                            ground, ordering)
    R, _ = PolynomialRing(
        ground,
        ["x$i" for i in 1:nvariables],
        ordering=ordering
    )

    return filter!(!iszero, [
        map_coefficients(
            c -> ground(data(c) % csz),
            rand(rng, R, exps, nterms)
        )
        for _ in 1:rand(rng, npolys)
    ])
end

function generate_set(nvariables, exps, nterms, npolys, csz, rng,
                        ground::AbstractAlgebra.Rationals, ordering)

    semiground = GF(2^31 - 1)
    R, _ = PolynomialRing(
        semiground,
        ["x$i" for i in 1:nvariables],
        ordering=ordering
    )

    csz = BigInt(csz)
    zzbase = BigInt(9223372036854775837)
    randzz() = rand(-zzbase:zzbase)

    return filter!(!iszero, [
        map_coefficients(
            c -> ground(mod(randzz(), csz), mod(randzz(), csz) + 1),
            rand(rng, R, exps, nterms)
        )
        for _ in 1:rand(rng, npolys)
    ])
end


text:
"""
    lorenztriple(;uᵢ=rand(9),
                σ₁ = 10.0, σ₂ = 10.0, σ₃ = 10.0,
                ρ₁ = 28.0, ρ₂ = 28.0, ρ₃ = 28.0,
                β₁ = 8/3,  β₂ = 8/3,  β₃ = 8.3,
				c₁ = 1.0, c₂ = 1.0) -> ContinuousDynamicalSystem

Equations of motion for three coupled Lorenz attractors with a transitive
causality chain where X₁ → X₂ and X₂ → X₃. In total, the three 3D-subsystems
create a 9-dimensional dynamical system.

The strength of the forcing X₁ → X₂ is controlled by the parameter `c₁`, and
the forcing from X₂ → X₃ by `c₂`. The remaining parameters are the usual
parameters for the Lorenz system, where the subscript `i` refers to the
subsystem Xᵢ. This system was studied by Papana et al. (2013) for coupling
strengths `c = 0, 1, 3, 5`.

# References
Papana et al., Simulation Study of Direct Causality Measures in Multivariate Time Series. Entropy 2013, 15(7), 2635-2661; doi:10.3390/e15072635
"""
function eom_lorenz_triple(u, p, t)
    x₁, y₁, z₁, x₂, y₂, z₂, x₃, y₃, z₃ = (u...,)
    σ₁, σ₂, σ₃, ρ₁, ρ₂, ρ₃, β₁, β₂, β₃, c₁, c₂ = (p...,)

    # Subsystem 1
    dx₁ = σ₁*(y₁-x₁)
    dy₁ = ρ₁*x₁ - y₁ - x₁*z₁
    dz₁ = x₁*y₁ - β₁*z₁

    # Subsystem 2
    dx₂ = σ₂*(y₂-x₂) + c₁*(x₁ - x₂)
    dy₂ = ρ₂*x₂ - y₂ - x₂*z₂
    dz₂ = x₂*y₂ - β₂*z₂

    # Subsystem 3
    dx₃ = σ₃*(y₃-x₃) + c₂*(x₂ - x₃)
    dy₃ = ρ₃*x₃ - y₃ - x₃*z₃
    dz₃ = x₃*y₃ - β₃*z₃
    return SVector{9}(dx₁, dy₁, dz₁, dx₂, dy₂,dz₂, dx₃, dy₃, dz₃)
end

function lorenz_triple(uᵢ, σ₁, σ₂, σ₃, ρ₁, ρ₂, ρ₃, β₁, β₂, β₃, c₁, c₂)
    p = [σ₁, σ₂, σ₃, ρ₁, ρ₂, ρ₃, β₁, β₂, β₃, c₁, c₂]
    ContinuousDynamicalSystem(eom_lorenz_triple, uᵢ, p)
end

"""
    lorenztriple(;uᵢ=rand(9),
                σ₁ = 10.0, σ₂ = 10.0, σ₃ = 10.0,
                ρ₁ = 28.0, ρ₂ = 28.0, ρ₃ = 28.0,
                β₁ = 8/3,  β₂ = 8/3,  β₃ = 8.3,
                c₁ = 1.0, c₂ = 1.0) -> ContinuousDynamicalSystem

Initalise a dynamical system consisting of three coupled Lorenz attractors with
a transitive causality chain where X₁ → X₂ and X₂ → X₃. In total, the three
3D-subsystems create a 9-dimensional dynamical system.

The strength of the forcing X₁ → X₂ is controlled by the parameter `c₁`, and
the forcing from X₂ → X₃ by `c₂`. The remaining parameters are the usual
parameters for the Lorenz system, where the subscript `i` refers to the
subsystem Xᵢ. This system was studied by Papana et al. (2013) for coupling
strengths `c = 0, 1, 3, 5`.

# References
Papana et al., Simulation Study of Direct Causality Measures in Multivariate Time Series. Entropy 2013, 15(7), 2635-2661; doi:10.3390/e15072635
"""
lorenz_triple(;uᵢ=rand(9),
            σ₁ = 10.0, σ₂ = 10.0, σ₃ = 10.0,
            ρ₁ = 28.0, ρ₂ = 28.0, ρ₃ = 28.0,
            β₁ = 8/3,  β₂ = 8/3,  β₃ = 8.3,
            c₁ = 1.0, c₂ = 1.0) =
    lorenz_triple(uᵢ, σ₁, σ₂, σ₃, ρ₁, ρ₂, ρ₃, β₁, β₂, β₃, c₁, c₂)

ext: jl


text:
using Documenter, Images

makedocs(doctest=false)

deploydocs(deps=Deps.pip("pygments", "mkdocs", "mkdocs-material"),
           repo="github.com/timholy/Images.jl.git",
           julia="0.4",
           osname="linux")


text:
# -*- encoding: utf-8 -*-
#
# MIT License
#
# Copyright (c) 2016 Maurizio Tomasi
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
#
# # TERMINOLOGY
#
# The code in this file tries to follow the conventions laid in
# document LSPE-STRIP-SP-017 ("STRIP Telescope Coordinate Systems
# Specification") as much as possible. As the reader might not be
# acquainted with it, here is a short summary of the relevant
# terminology:
#
# +-----+-----------------------------+------------------------------------------+
# | SCS | Site Coordinate System      | Centered on the GPS station              |
# +-----+-----------------------------+------------------------------------------+
# | MCS | Mount Coordinate System     | Centered at the basis of the telescope   |
# +-------+---------------------------+------------------------------------------+
# | ACS | A zimuth Coordinate System  | Moveable, centered in the middle of      |
# |     |                             | ground motor. Positive angles go from    |
# |     |                             | North to *East* (beware!)                |
# +-----+-----------------------------+------------------------------------------+
# | ECS | Elevation Coordinate System | Moveable, the Y axis is aligned with     |
# |     |                             | the motor axis                           |
# +-----+-----------------------------+------------------------------------------+
# | BCS | Boresight Coordinate System | It should have been moveable, but the    |
# |     |                             | current design does not allow it (alas!) |
# +-----+-----------------------------+------------------------------------------+
# | RDP | Reference Detector Plane    | Centered in the focus of the telescope.  |
# |     |                             | All the beam vectors are specified in    |
# |     |                             | this frame.                              |
# +-----+-----------------------------+------------------------------------------+
#
# We always use a right-handed coordinate system:
#
#                 ^ Z axis (aligned with "up" direction, e.g., Zenith)
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#          Origin O----------------> Y axis
#                /
#               /
#              /
#             /
#            /
#           v  X axis (aligned with North, when applicable)
#

import Quaternions
import Healpix
import StaticArrays
import LinearAlgebra: ×, dot
import AstroLib
import Dates

export TENERIFE_LATITUDE_DEG, TENERIFE_LONGITUDE_DEG, TENERIFE_HEIGHT_M
export qrotation_x, qrotation_y, qrotation_z
export timetorotang, telescopetoground, groundtoearth
export genpointings, polarizationangle

const TENERIFE_LATITUDE_DEG = 28.3
const TENERIFE_LONGITUDE_DEG = -16.509722
const TENERIFE_HEIGHT_M = 2390

# These functions are faster than Quaternions.qrotation2

function qrotation_x(theta)
    Quaternions.Quaternion(cos(theta / 2), sin(theta / 2), 0.0, 0.0, true)
end

function qrotation_y(theta)
    Quaternions.Quaternion(cos(theta / 2), 0.0, sin(theta / 2), 0.0, true)
end

function qrotation_z(theta)
    Quaternions.Quaternion(cos(theta / 2), 0.0, 0.0, sin(theta / 2), true)
end

"""
    qrotation_x(theta)
    qrotation_y(theta)
    qrotation_z(theta)

Return a `Quaternions.Quaternion` object representing a rotation
around the ``e_x``, ``e_y``, or ``e_z`` axis by an angle `theta` (in
radians).
"""
qrotation_x, qrotation_y, qrotation_z

"""
    timetorotang(time, rpm)

Convert a time into a rotation angle, given the number of rotations per minute.
The time should be expressed in seconds. The return value is in radians.
`time` can either be a scalar or a vector.
"""
function timetorotang(time_s, rpm)
    if rpm == 0
        0.0
    else
        2 * π * time_s * (rpm / 60)
    end
end

@doc raw"""
    telescopetoground(wheelanglesfn, time_s)

Return a quaternion of type `Quaternion{Float64}` representing the
coordinate transform from the focal plane to the ground of the
telescope. The parameter `wheelanglesfn` must be a function which
takes as input a time, `time_s`, in seconds, and it must return a
3-tuple containing the angles of the following motors, measured in
**radians**:

1. The boresight motor (rotation around the ``z`` axis,
   counterclockwise)

2. The altitude motor (rotation around the ``y`` axis,
   counterclockwise)

3. The ground motor (rotation around the ``z`` axis, **clockwise**:
   N→E→S→W)

# Example

`````julia
telescopetoground(3600.0) do
    # Boresight motor keeps a constant angle equal to 0°
    # Altitude motor remains at 20° from the Zenith
    # Ground motor spins at 1 RPM
    (0.0, deg2rad(20.0), timetorotang(time_s, 1))
end
`````
"""
function telescopetoground(wheelanglesfn, time_s)
    (wheel1ang, wheel2ang, wheel3ang) = wheelanglesfn(time_s)

    qwheel1 = qrotation_z(wheel1ang)
    qwheel2 = qrotation_y(wheel2ang)

    # The minus sign here takes into account the fact that the azimuth
    # motor requires positive angles to turn North into East
    qwheel3 = qrotation_z(-wheel3ang)

    qwheel3 * (qwheel2 * qwheel1)
end

"""
    groundtoearth(groundq, time_s, latitude_deg; day_duration_s=86400.0)

Return a quaternion of type `Quaternion{Float64}` representing the
coordinate transformation from the ground of the telescope to the
Equatorial coordinate system. The parameter `groundq` must be a
quaternion describing the coordinate transformation from the focal
plane of the telescope to the ground. The parameter `time_s` must be a
time in seconds, and `latitude_deg` is the latitude (in degrees, N is
positive) of the location where the observation is made.

The keyword `day_duration_s` specifies the length of a sidereal day in
seconds.
"""
function groundtoearth(groundq, time_s, latitude_deg; day_duration_s=86400.0)
    locq = qrotation_y(deg2rad(90 - latitude_deg))
    earthq = qrotation_z(2π * time_s / day_duration_s)

    earthq * (locq * groundq)
end

"""
    vector2equatorial(dir, jd, latitude_deg, longitude_deg, height_m;
                      prec=true, nut=true, aber=true)

Transform the Healpix coordinates of a vector into equatorial coordinates.
The parameter `vector` is 3D vector, `jd` is the julian date. The
paramters `latitude_deg`, `longitude_deg` and `height_m` should contain the
latitude (in degrees, N is positive), the longitude (in degrees, counterclockwise
is positive) and the height (in meters) of the location where the observation is
made. The parameters `prec`, `nut` and `aber` allow to consider secondary effects
respectively of precession, nutation and aberration.
"""
function vector2equatorial(vector, jd, latitude_deg, longitude_deg, height_m,
                           prec, nut, aber, ref)
    (θ, ϕ) = Healpix.vec2ang(vector...)
    alt_rad = π/2 - θ
    az_rad = 2π - ϕ

    ra_deg, dec_deg, _ = AstroLib.hor2eq(rad2deg(alt_rad),
                                         rad2deg(az_rad),
                                         jd,
                                         latitude_deg,
                                         longitude_deg,
                                         height_m,
                                         precession=prec,
                                         nutate=nut,
                                         aberration=aber,
                                         refract=ref)
    (deg2rad(dec_deg), deg2rad(ra_deg))
end

"""
    polarizationangle(northdir, poldir)

Calculate the polarization angle projected in the sky in IAU
conventions.  The parameter `northdir` must be a versor that points
the North and `poldir` must be a versor that identify the polarization
direction projected in the sky.  The return value is in radians.
"""
function polarizationangle(northdir, poldir)
    cosψ = clamp(dot(northdir, poldir), -1, 1)
    crosspr = northdir × poldir
    sinψ = clamp(sqrt(dot(crosspr, crosspr)), -1, 1)
    ψ = atan(sinψ, cosψ)
    ψ
end

"""
    get_northdir(θ, ϕ)

Compute the North of a vector. The North for a vector v is -dv/dθ, as
θ is the colatitude and moves along the meridian.
"""
get_northdir(θ, ϕ) = StaticArrays.@SArray [-cos(θ) * cos(ϕ), -cos(θ) * sin(ϕ), sin(θ)]

"""
    quat_to_angles(boreaxis, polaxis, quat)

Transform the boresight direction `dir` and the polarization direction
`poldir` according to quaternion `quat`. Return the 3-tuple (θ, ϕ, ψ)
representing the colatitude, longitude, and polarization angle
(calculated northward).

This function is used internally by [`genpointings`](@ref).
"""
function quat_to_angles(boreaxis, polaxis, quat)
    rotmatr = Quaternions.rotationmatrix(quat)
    boresight = rotmatr * boreaxis
    poldir = rotmatr * polaxis

    (θ, ϕ) = Healpix.vec2ang(boresight...)

    northdir = get_northdir(θ, ϕ)
    (θ, ϕ, polarizationangle(northdir, poldir))
end

function genpointings(wheelanglesfn,
                      dir,
                      timerange_s;
                      polaxis=Float64[1.0, 0.0, 0.0],
                      latitude_deg=0.0,
                      ground=false,
                      day_duration_s=86400.0)

    if ground
        dirs = Array{Float64}(undef, length(timerange_s), 4)
        ψ = Array{Float64}(undef, length(timerange_s), 2)
    else
        dirs = Array{Float64}(undef, length(timerange_s), 2)
        ψ = Array{Float64}(undef, length(timerange_s))
    end

    for (idx, time_s) = enumerate(timerange_s)

        # This converts the RDP into the MCS (ground reference frame)
        groundq = telescopetoground(wheelanglesfn, time_s)
        # This converts the MCS into the celestial reference frame
        quat = groundtoearth(groundq, time_s, latitude_deg; day_duration_s=day_duration_s)

        θ, ϕ, curψ = quat_to_angles(dir, polaxis, quat)
        (dirs[idx, 1], dirs[idx, 2]) = (θ, ϕ)

        if ground
            # Re-run the transformation algorithm using the ground quaternion
            θ_ground, ϕ_ground, ψ_ground = quat_to_angles(dir, polaxis, groundq)

            (dirs[idx, 3], dirs[idx, 4]) = (θ_ground, ϕ_ground)
            (ψ[idx, 1], ψ[idx, 2]) = (curψ, ψ_ground)
        else
            ψ[idx] = curψ
        end
    end

    (dirs, ψ)
end

function genpointings(wheelanglesfn,
                      dir,
                      timerange_s,
                      t_start;
                      polaxis=Float64[1.0, 0.0, 0.0],
                      latitude_deg=0.0,
                      longitude_deg=0.0,
                      height_m=0,
                      precession=true,
                      nutation=true,
                      aberration=true,
                      refract=true)

    skydirs = Array{Float64}(undef, length(timerange_s), 2)
    skyψ = Array{Float64}(undef, length(timerange_s))

    for (idx, time_s) = enumerate(timerange_s)
        groundq = telescopetoground(wheelanglesfn, time_s)
        rotmatr = Quaternions.rotationmatrix(groundq)
        vector = rotmatr * dir

        jd = AstroLib.jdcnv(t_start + Dates.Nanosecond(round(Int64, time_s * 1e9)))
        Dec_rad, Ra_rad = vector2equatorial(vector,
                                            jd,
                                            latitude_deg,
                                            longitude_deg,
                                            height_m,
                                            precession,
                                            nutation,
                                            aberration,
                                            refract)

        poldir = rotmatr * polaxis
        northdir = get_northdir(π/2 - Dec_rad, Ra_rad)

        skydirs[idx, 1] = Dec_rad
        skydirs[idx, 2] = Ra_rad
        skyψ[idx] = polarizationangle(northdir, poldir)
    end

    (skydirs, skyψ)
end

@doc raw"""
    genpointings(wheelanglesfn, dir, timerange_s;
                 polaxis=Float64[1.0, 0.0, 0.0],
                 latitude_deg=0.0, ground=false)
    genpointings(wheelanglesfn, dir, timerange_s, t_start;
                 polaxis=Float64[1.0, 0.0, 0.0],
                 latitude_deg=0.0, longitude_deg=0.0, height_m=0.0)

Generate a set of pointings for a STRIP detector. The parameter
`wheelanglesfn` must be a function which takes as input a time in
seconds and returns a 3-tuple containing the angles (in radians) of
the three motors:

1. The boresight motor (rotation around the ``z`` axis,
   counterclockwise)

2. The altitude motor (rotation around the ``y`` axis,
   counterclockwise)

3. The ground motor (rotation around the ``z`` axis, **clockwise**:
   N→E→S→W)

The meaning of the parameters/keywords is the following:

- `dir` specifies the pointing direction of the mean (the boresight is
  [0, 0, 1]). It must be normalized.

- `timerange_s` is an enumerable type that specifies at which times
  (in seconds) pointings must be computed.

- `t_start` is a DateTime which tells the exact UTC date and time of the
  observation.

- `latitude_deg` is the latitude of the location where the observation
  is made (in degrees, North is positive).

- `ground` is a Boolean: if `true`, the function will return a 4-tuple
  containing the colatitude and longitude measured in Equatorial
  coordinates (columns 1 and 2) and in ground coordinates (columns 3
  and 4). If false, only the Equatorial coordinates are computed.

- `polaxis` is the polarization axis; it must be normalized.

- `precession` is a boolen parameter that allow to include the Earth's precession
  effect; default is true.

- `nutation` is a boolen parameter that allow to include the Earth's nutation
  effect; default is true.

- `aberration` is a boolen parameter that allow to include the stellar aberration
  effect; default is true.

- `refract` is a boolen parameter that allow to include the refraction correction
  ; default is true.

# Return values

In the first form, the function returns a 2-tuple containing the sky
directions (a N×2 array containing declination and right ascension, in
Equatorial coordinates) and the polarization angle for each time step.

In the second form, the function returns a 2-tuple (4-tuple)
containing the directions (a N×2 or Nx4 array containing the
colatitude and the longitude) and the polarization angles at each time
step.

# Examples

Here is an example of the first kind:

`````julia
dir, psi = genpointings([0, 0, 1], 0:0.1:1) do time_s
    # Boresight motor keeps a constant angle equal to 0°
    # Altitude motor remains at 20° from the Zenith
    # Ground motor spins at 1 RPM
    return (0.0, deg2rad(20.0), timetorotang(time_s, 1))
end
`````

And here is an example of the second kind. Note that, unlike the
previous example, we use a lambda function instead of the `do...end`
notation.

`````julia
import Dates

dirs, psi = genpointings(time_s -> (0, deg2rad(20),
                                    timetorotang(time_s, 1)),
                         [0, 0, 1],
                         0:0.1:1,
                         Dates.DateTime(2019, 01, 01, 0, 0, 0),
                         latitude_deg=10.0,
                         longitude_deg=20.0,
                         height_m = 1000) do time_s
`````
"""
genpointings


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

# Import QGIS-cropped layer as example
ref = geotiff(SimpleSDMPredictor, "../../rasters/wc1_cropped.tif")

## Step 3: Update match_longitude for range not LinRange ####

# Load things
layer = SimpleSDMPredictor(WorldClim, BioClim, 1)
layer_subset = layer[coords]
size(layer_subset) # should be 330, 570
left, right, bottom, top = coords
import SimpleSDMLayers: _match_latitude, _match_longitude

# Preliminary tests
lon = left
@time _match_longitude(layer, left; side=:left) # 15k alloc
ldiff = abs.(lon .- longitudes(layer))
@time any(x -> isapprox(x, stride(layer, 1)), ldiff)
@time any(x -> isequal(x, stride(layer, 1)), ldiff) # not much difference
@time findlast(x -> isapprox(x, stride(layer, 1)), ldiff)
@time findlast(x -> isequal(x, stride(layer, 1)), ldiff) # not much difference either
@time findfirst(x -> isequal(x, stride(layer, 1)), ldiff)
@time findall(x -> isequal(x, stride(layer, 1)), ldiff)
@time filter(x -> x == stride(layer, 1), ldiff)
@time lon in latitudes(layer)
@time lon in ldiff
@time any(isequal.(ldiff, stride(layer, 1)))
@time any(isapprox.(ldiff, stride(layer, 1)))
@time isequal.(ldiff, stride(layer,1))
# With isequal.
@time _match_longitude(layer, left; side=:left) # 13k
@time _match_longitude(layer, left; side=:none) # 29 alloc (no change)
# with lequal
@time _match_longitude(layer, left; side=:left) # 39 alloc
@time _match_longitude(layer, left; side=:none) # 40 alloc
@time _match_longitude(layer, left; side=:right) # 39 alloc
@time _match_longitude(layer, left+0.001; side=:right) # 41 alloc
# with lapprox
@time _match_longitude(layer, left; side=:left) # 39 alloc
@time _match_longitude(layer, left; side=:none) # 40 alloc
@time _match_longitude(layer, left; side=:right) # 39 alloc
@time _match_longitude(layer, left+0.001; side=:right) # 41 alloc
# Check with latitude
@time _match_latitude(layer, bottom; side=:none) # 40 alloc
@time _match_latitude(layer, bottom; side=:bottom) # 39 alloc
@time _match_latitude(layer, top; side=:top) # 39 alloc
@time _match_latitude(layer, top+0.001; side=:top) # 41 alloc
# Check getindex result
@time layer[coords] # 278 alloc!
# compared to 47.69 k with previous fix, and 222 alloc before previous fix


text:
module SpectralClustering
export EigenvectorEmbedder,
       embedding

abstract type EigenvectorEmbedder end

include("Utils/DataAccess.jl")
include("Utils/DataProcessing.jl")
include("Graph/Graphs.jl")
include("LandmarkSelection.jl")
include("Embedding.jl")
include("ApproximateEmbedding.jl")
include("MultiView.jl")
include("EigenvectorClustering.jl")
end


text:
module RegEdit

export RegBranch
export register

using AutoHashEquals
using LibGit2
using Pkg: Pkg, TOML, GitTools
using UUIDs

const DEFAULT_REGISTRY_URL = "https://github.com/JuliaRegistries/General"

include("builtin_pkgs.jl")
include("types.jl")
include("register.jl")
include("utils.jl")

end


text:
"""
    InitialConditions

Functions and types for generating initial conditions for simulations.
"""
module InitialConditions

include("QuantisedDiatomic/QuantisedDiatomic.jl")
export QuantisedDiatomic

include("ThermalMonteCarlo.jl")
export ThermalMonteCarlo

include("MetropolisHastings.jl")
export MetropolisHastings

end # module


text:
# we collect the list of unit test files here

include("TestHullWhiteModel.jl")
include("TestHybridModel.jl")
include("TestPayoffScripting.jl")
include("TestQuasiGaussianModel.jl")

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
function converging!(ip::Vector{Int64}, population::Matrix{Float64},
                        n_population::Int64, n_gene::Int64)::Matrix{Float64}
    n_children::Int8 = 10
    children::Matrix{Float64} = zeros(n_children, n_gene+1)
    for i in 1:n_children
        ip[3:end] = sample(
            collect(1:n_population), n_gene, replace=false
        )
        children[i,:] = crossover(population[ip,:], n_gene)
    end

    family::Matrix{Float64} = zeros(n_children+2, n_gene+1)
    @inbounds for i in 1:n_gene+1
        @simd for j in 1:n_children
            family[j,i] = children[j,i]
        end
        family[n_children+1,i] = population[ip[1],i]
        family[n_children+2,i] = population[ip[2],i]
    end

    family = sortslices(family, dims=1, by=x->x[end])

    ic2::Int8 = rand(2:n_children+2)
    @inbounds for i in 1:n_gene+1
        population[ip[1],i] = family[1,i]  # Best
        population[ip[2],i] = family[ic2,i]  # Random
    end

    if !isfinite(population[ip[2],end])
        population[ip[2],end] = objective(
            population[ip[2],1:n_gene]
        )
    end

    population = sortslices(population, dims=1, by=x->x[end])

    return population
end

function crossover(parents::Matrix{Float64}, n_gene::Int64)::Vector{Float64}
    child::Vector{Float64} = ENDX(parents, n_gene)
    for i in 1:n_gene
        @inbounds child[i] = clamp(child[i], 0.0, 1.0)
    end

    child[end] = Inf

    return child
end

# Extended Normal Distribution Xover
function ENDX(parents::Matrix{Float64}, n_gene::Int64)::Vector{Float64}
    child::Vector{Float64} = zeros(n_gene+1)

    ALPHA::Float64 = sqrt((1.0-2.0*0.35^2.0)) / 2.0
    BETA::Float64 = 0.35 / sqrt(n_gene-1)

    t1::Vector{Float64} = (parents[2,1:n_gene] - parents[1,1:n_gene])./2.0
    t2::Vector{Float64} = randn()*ALPHA*(parents[2,1:n_gene] - parents[1,1:n_gene])
    t3::Vector{Float64} = zeros(n_gene)
    centroid::Vector{Float64} = reshape(
        mean(parents[3:end,1:n_gene], dims=1), n_gene
    )

    for i in 1:n_gene
        t3 += randn()*BETA*(parents[i+2,1:n_gene] - centroid)
    end

    for i in 1:n_gene
        @inbounds child[i] = t1[i] + t2[i] + t3[i]
    end

    return child
end


text:
################################################################################
### Lattice via LatticePhysics
################################################################################

using LinearAlgebra, LatticePhysics, LatPhysUnitcellLibrary

function LatPhysUnitcellLibrary.getUnitcellSquare(
            unitcell_type  :: Type{U},
            implementation :: Val{17}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,

        # Bravais lattice vectors
        [[1.0, +1.0], [1.0, -1.0]],

        # Sites
        S[
            newSite(S, [0.0, 0.0], getDefaultLabelN(LS, 1)),
            newSite(S, [0.0, 1.0], getDefaultLabelN(LS, 2))
        ],

        # Bonds
        B[
            # NN, directed
            # bonds from ref plot, π/4 weight for spin up
            newBond(B, 1, 2, getDefaultLabelN(LB, 1), (0, 1)),
            newBond(B, 1, 2, getDefaultLabelN(LB, 1), (-1, 0)),
            newBond(B, 2, 1, getDefaultLabelN(LB, 1), (+1, -1)),
            newBond(B, 2, 1, getDefaultLabelN(LB, 1), (0, 0)),

            # NN reversal
            newBond(B, 2, 1, getDefaultLabelN(LB, 2), (0, -1)),
            newBond(B, 2, 1, getDefaultLabelN(LB, 2), (+1, 0)),
            newBond(B, 1, 2, getDefaultLabelN(LB, 2), (-1, +1)),
            newBond(B, 1, 2, getDefaultLabelN(LB, 2), (0, 0)),

            # NNN
            # positive weight (we need forward and backward facing bonds here too)
            newBond(B, 1, 1, getDefaultLabelN(LB, 3), (+1, 0)),
            newBond(B, 1, 1, getDefaultLabelN(LB, 3), (-1, 0)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 3), (0, +1)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 3), (0, -1)),
            # negative weight
            newBond(B, 1, 1, getDefaultLabelN(LB, 4), (0, +1)),
            newBond(B, 1, 1, getDefaultLabelN(LB, 4), (0, -1)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 4), (+1, 0)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 4), (-1, 0)),

            # Fifth nearest neighbors
            newBond(B, 1, 1, getDefaultLabelN(LB, 5), (2, 0)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 5), (2, 0)),
            newBond(B, 1, 1, getDefaultLabelN(LB, 5), (0, 2)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 5), (0, 2)),
            # backwards facing bonds
            newBond(B, 1, 1, getDefaultLabelN(LB, 5), (-2, 0)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 5), (-2, 0)),
            newBond(B, 1, 1, getDefaultLabelN(LB, 5), (0, -2)),
            newBond(B, 2, 2, getDefaultLabelN(LB, 5), (0, -2)),
        ]
    )
end

################################################################################
### Model implementation
################################################################################

using MonteCarlo
using MonteCarlo: StructArray, CMat64, CVec64
using MonteCarlo: AbstractLattice, @bm, HubbardConf, @turbo, BlockDiagonal, vmul!, conf

"""
t5 = 0  <=> F = 0.2
t5 = (1 - sqrt(2)) / 4 <=> F = 0.009
"""
MonteCarlo.@with_kw_noshow struct HBCModel{LT<:AbstractLattice} <: HubbardModel
    # parameters with defaults based on paper
    mu::Float64 = 0.0
    U::Float64 = 1.0
    @assert U >= 0. "U must be positive."
    t1::Float64 = 1.0
    t2::Float64 = 1.0 / sqrt(2.0)
    t5::Float64 = (1 - sqrt(2)) / 4

    # lattice
    l::LT

    # two fermion flavors (up, down)
    flv::Int = 2

    # temp storage to avoid allocations in propose_local and accept_local
    IG::CMat64  = StructArray(Matrix{ComplexF64}(undef, length(l), 2))
    IGR::CMat64 = StructArray(Matrix{ComplexF64}(undef, length(l), 2))
    R::Diagonal{ComplexF64, CVec64} = Diagonal(StructArray(Vector{ComplexF64}(undef, 2)))
end

# Constructors
HBCModel(params::Dict{Symbol}) = HBCModel(; params...)
HBCModel(params::NamedTuple) = HBCModel(; params...)
function HBCModel(lattice::AbstractLattice; kwargs...)
    HBCModel(l = lattice; kwargs...)
end
function HBCModel(L, dims; kwargs...)
    l = choose_lattice(HBCModel, dims, L)
    HBCModel(l = l; kwargs...)
end

# type information for DQMC stack
MonteCarlo.hoppingeltype(::Type{DQMC}, ::HBCModel) = ComplexF64
MonteCarlo.hopping_matrix_type(::Type{DQMC}, ::HBCModel) = BlockDiagonal{ComplexF64, 2, CMat64}
MonteCarlo.greenseltype(::Type{DQMC}, ::HBCModel) = ComplexF64
MonteCarlo.greens_matrix_type( ::Type{DQMC}, ::HBCModel) = BlockDiagonal{ComplexF64, 2, CMat64}

# cosmetics
import Base.summary
import Base.show
Base.summary(model::HBCModel) = "Hofmann Berg Hubbard model"
function Base.show(io::IO, model::HBCModel)
    print(io, "Hofmann Berg Hubbard model, $(length(model.l)) sites")
end
Base.show(io::IO, m::MIME"text/plain", model::HBCModel) = print(io, model)

# Convenience
@inline MonteCarlo.parameters(m::HBCModel) = (N = length(m.l), t1 = m.t1, t2 = m.t2, t5 = m.t5, U = -m.U, mu = m.mu)

# generate hopping matrix
function MonteCarlo.hopping_matrix(mc::DQMC, m::HBCModel{<: LatPhysLattice})
    # number of sites
    N = length(m.l)

    # spin up and spin down blocks of T
    tup = diagm(0 => fill(-ComplexF64(m.mu), N))
    tdown = diagm(0 => fill(-ComplexF64(m.mu), N))

    # positive and negative prefactors for t1, t2
    t1p = m.t1 * cis(+pi/4) # ϕ_ij^↑ = + π/4
    t1m = m.t1 * cis(-pi/4) # ϕ_ij^↓ = - π/4
    t2p = + m.t2
    t2m = - m.t2

    for b in bonds(m.l.lattice)
        # NN paper direction
        if b.label == 1
            tup[b.from, b.to]   = - t1p
            tdown[b.from, b.to] = - t1m

        # NN reverse direction
        elseif b.label == 2
            tup[b.from, b.to]   = - t1m
            tdown[b.from, b.to] = - t1p

        # NNN solid bonds
        elseif b.label == 3
            tup[b.from, b.to]   = - t2p
            tdown[b.from, b.to] = - t2p

        # NNN dashed bonds
        elseif b.label == 4
            tup[b.from, b.to]   = - t2m
            tdown[b.from, b.to] = - t2m

        # Fifth nearest neighbors
        else
            tup[b.from, b.to]   = - m.t5
            tdown[b.from, b.to] = - m.t5
        end
    end

    return BlockDiagonal(StructArray(tup), StructArray(tdown))
end

"""
Calculate the interaction matrix exponential `expV = exp(- power * delta_tau * V(slice))`
and store it in `result::Matrix`.

This is a performance critical method.
"""
@inline @bm function MonteCarlo.interaction_matrix_exp!(mc::DQMC, m::HBCModel,
            result::Diagonal, conf::HubbardConf, slice::Int, power::Float64=1.)
    dtau = mc.parameters.delta_tau
    lambda = acosh(exp(0.5 * m.U * dtau))
    N = length(lattice(m))

    # spin up block
    @inbounds for i in 1:N
        result.diag[i] = exp(sign(power) * lambda * conf[i, slice])
    end

    # spin down block
    @inbounds for i in 1:N
        result.diag[N+i] = exp(sign(power) * lambda * conf[i, slice])
    end
    nothing
end

@inline @bm function MonteCarlo.propose_local(
        mc::DQMC, model::HBModel, i::Int, slice::Int, conf::HubbardConf
    )
    N = length(model.l)
    G = mc.stack.greens
    Δτ = mc.parameters.delta_tau
    R = model.R

    α = acosh(exp(0.5Δτ * model.U))
    ΔE_boson = -2.0α * conf[i, slice]
    Δ = exp(ΔE_boson) - 1.0

    # Unrolled R = I + Δ * (I - G)
    # up-up term
    R.diag.re[1] = 1.0 + Δ * (1.0 - G.blocks[1].re[i, i])
    R.diag.im[1] = - Δ * G.blocks[1].im[i, i]
    # down-down term
    R.diag.re[2] = 1.0 + Δ * (1.0 - G.blocks[2].re[i, i])
    R.diag.im[2] = - Δ * G.blocks[2].im[i, i]

    # Calculate "determinant"
    detratio = ComplexF64(
        R.diag.re[1] * R.diag.re[2] - R.diag.im[1] * R.diag.im[2],
        R.diag.re[1] * R.diag.im[2] + R.diag.im[1] * R.diag.re[2]
    )

    return detratio, ΔE_boson, Δ
end

@inline @bm function MonteCarlo.accept_local!(
        mc::DQMC, model::HBModel, i::Int, slice::Int, conf::HubbardConf,
        detratio, ΔE_boson, Δ)

    @bm "accept_local (init)" begin
        N = length(model.l)
        G = mc.stack.greens
        IG = model.IG
        IGR = model.IGR
        R = model.R
    end

    # compute R⁻¹ Δ, using that R is Diagonal, Δ is Number
    # using Δ / (a + ib) = Δ / (a^2 + b^2) * (a - ib)
    @bm "accept_local (inversion)" begin
        f = Δ / (R.diag.re[1]^2 + R.diag.im[1]^2)
        R.diag.re[1] = +f * R.diag.re[1]
        R.diag.im[1] = -f * R.diag.im[1]
        f = Δ / (R.diag.re[2]^2 + R.diag.im[2]^2)
        R.diag.re[2] = +f * R.diag.re[2]
        R.diag.im[2] = -f * R.diag.im[2]
    end

    # Compute (I - G) R^-1 Δ
    # Note IG is reduced to non-zero entries. Full IG would be
    # (I-G)[:, i]        0
    #     0         (I-G)[:, i+N]
    # our IG is [(I-G)[:, i]  (I-G)[:, i+N]]
    @bm "accept_local (IG, R)" begin
        # Calculate IG = I - G (relevant entries only)
        @turbo for m in axes(IG, 1)
            IG.re[m, 1] = -G.blocks[1].re[m, i]
        end
        @turbo for m in axes(IG, 1)
            IG.re[m, 2] = -G.blocks[2].re[m, i]
        end
        @turbo for m in axes(IG, 1)
            IG.im[m, 1] = -G.blocks[1].im[m, i]
        end
        @turbo for m in axes(IG, 1)
            IG.im[m, 2] = -G.blocks[2].im[m, i]
        end
        IG.re[i, 1] += 1.0
        IG.re[i, 2] += 1.0

        # Calculate IGR = IG * R where R = R⁻¹ Δ from the
        # previous calculation (relevant entries only)
        # spin up-up block
        @turbo for m in axes(IG, 1)
            IGR.re[m, 1] = IG.re[m, 1] * R.diag.re[1]
        end
        @turbo for m in axes(IG, 1)
            IGR.re[m, 1] -= IG.im[m, 1] * R.diag.im[1]
        end
        @turbo for m in axes(IG, 1)
            IGR.im[m, 1] = IG.re[m, 1] * R.diag.im[1]
        end
        @turbo for m in axes(IG, 1)
            IGR.im[m, 1] += IG.im[m, 1] * R.diag.re[1]
        end

        # spin down-down block
        @turbo for m in axes(IG, 1)
            IGR.re[m, 2] = IG.re[m, 2] * R.diag.re[2]
        end
        @turbo for m in axes(IG, 1)
            IGR.re[m, 2] -= IG.im[m, 2] * R.diag.im[2]
        end
        @turbo for m in axes(IG, 1)
            IGR.im[m, 2] = IG.re[m, 2] * R.diag.im[2]
        end
        @turbo for m in axes(IG, 1)
            IGR.im[m, 2] += IG.im[m, 2] * R.diag.re[2]
        end
    end

    # Update G according to G = G - (I - G)[:, i:N:end] * R⁻¹ * Δ * G[i:N:end, :]
    # We already have IG = (I - G)[:, i:N:end] * R⁻¹ * Δ
    @bm "accept_local (finalize computation)" begin
        # get blocks to write less
        G1 = G.blocks[1]
        G2 = G.blocks[2]
        temp1 = mc.stack.greens_temp.blocks[1]
        temp2 = mc.stack.greens_temp.blocks[2]

        # compute temp = IG[:, i:N:end] * G[i:N:end, :]
        # spin up-up block
        @turbo for m in axes(G1, 1), n in axes(G1, 2)
            temp1.re[m, n] = IGR.re[m, 1] * G1.re[i, n]
        end
        @turbo for m in axes(G1, 1), n in axes(G1, 2)
            temp1.re[m, n] -= IGR.im[m, 1] * G1.im[i, n]
        end
        @turbo for m in axes(G1, 1), n in axes(G1, 2)
            temp1.im[m, n] = IGR.im[m, 1] * G1.re[i, n]
        end
        @turbo for m in axes(G1, 1), n in axes(G1, 2)
            temp1.im[m, n] += IGR.re[m, 1] * G1.im[i, n]
        end

        # spin down-down block
        @turbo for m in axes(G2, 1), n in axes(G2, 2)
            temp2.re[m, n] = IGR.re[m, 2] * G2.re[i, n]
        end
        @turbo for m in axes(G2, 1), n in axes(G2, 2)
            temp2.re[m, n] -= IGR.im[m, 2] * G2.im[i, n]
        end
        @turbo for m in axes(G2, 1), n in axes(G2, 2)
            temp2.im[m, n] = IGR.im[m, 2] * G2.re[i, n]
        end
        @turbo for m in axes(G2, 1), n in axes(G2, 2)
            temp2.im[m, n] += IGR.re[m, 2] * G2.im[i, n]
        end

        # Calculate G = G - temp
        # spin up-up block
        @turbo for m in axes(G1, 1), n in axes(G1, 2)
            G1.re[m, n] = G1.re[m, n] - temp1.re[m, n]
        end
        @turbo for m in axes(G1, 1), n in axes(G1, 2)
            G1.im[m, n] = G1.im[m, n] - temp1.im[m, n]
        end

        # spin down-down block
        @turbo for m in axes(G2, 1), n in axes(G2, 2)
            G2.re[m, n] = G2.re[m, n] - temp2.re[m, n]
        end
        @turbo for m in axes(G2, 1), n in axes(G2, 2)
            G2.im[m, n] = G2.im[m, n] - temp2.im[m, n]
        end

        # Update configuration
        conf[i, slice] *= -1
    end

    nothing
end

# this enables boson energy measurements and global updates
@inline function MonteCarlo.energy_boson(mc::DQMC, m::HBCModel, hsfield = conf(mc))
    dtau = mc.parameters.delta_tau
    lambda = acosh(exp(m.U * dtau/2))
    return lambda * sum(hsfield)
end

# to save only what is necessary and not rely on the model type remaining the same
function MonteCarlo.save_model(
        file::MonteCarlo.JLDFile, m::HBCModel,
        entryname::String="Model"
    )
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "HBCModel")

    write(file, entryname * "/mu", m.mu)
    write(file, entryname * "/U", m.U)
    write(file, entryname * "/t1", m.t1)
    write(file, entryname * "/t2", m.t2)
    write(file, entryname * "/t5", m.t5)
    MonteCarlo.save_lattice(file, m.l, entryname * "/l")
    write(file, entryname * "/flv", m.flv)

    nothing
end

# I used HBModel as a type name and tag before. Adding this will alow the old
# tag to still work
MonteCarlo._load(data, ::Val{:HBModel}) = MonteCarlo._load(data, Val(:HBCModel))
function MonteCarlo._load(data, ::Val{:HBCModel})
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load HBCModel version $(data["VERSION"])"))
    end

    l = MonteCarlo._load(data["l"], MonteCarlo.to_tag(data["l"]))
    HBCModel(
        mu = data["mu"],
        U = data["U"],
        t1 = data["t1"],
        t2 = data["t2"],
        t5 = data["t5"],
        l = l,
        flv = data["flv"]
    )
end

################################################################################
### Measurement kernels
################################################################################

MonteCarlo.checkflavors(model::HBCModel) = MonteCarlo.checkflavors(model, 2)

function MonteCarlo.intE_kernel(mc, model::HBCModel, G::GreensMatrix)
    # ⟨U (n↑ - 1/2)(n↓ - 1/2)⟩ = ...
    # = U [G↑↑ G↓↓ - G↓↑ G↑↓ - 0.5 G↑↑ - 0.5 G↓↓ + G↑↓ + 0.25]
    # = U [(G↑↑ - 1/2)(G↓↓ - 1/2) + G↑↓(1 + G↑↓)]
    # with up-up = down-down and up-down = 0
    - model.U * sum((diag(G.val) .- 0.5).^2)
end


text:
# # Barotropic gyre

using Oceananigans
using Oceananigans.Grids

using Oceananigans.Coriolis:
    HydrostaticSphericalCoriolis,
    VectorInvariantEnergyConserving,
    VectorInvariantEnstrophyConserving

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    HydrostaticFreeSurfaceModel,
    VectorInvariant,
    ExplicitFreeSurface,
    ImplicitFreeSurface

using Oceananigans.TurbulenceClosures: HorizontallyCurvilinearAnisotropicDiffusivity
using Oceananigans.Utils: prettytime, hours, day, days, years
using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval, IterationInterval

using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary, RasterDepthMask

using Statistics
using JLD2
using Printf

Nx = 60
Ny = 60

# A spherical domain
underlying_grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                        longitude = (-30, 30),
                                        latitude = (15, 75),
                                        z = (-4000, 0))

@inline raster_depth(i, j) = 30 < i < 35 && 42 < j < 48

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBoundary(raster_depth, mask_type=RasterDepthMask()))

free_surface = ImplicitFreeSurface(gravitational_acceleration=0.1)
# free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1)

coriolis = HydrostaticSphericalCoriolis(scheme = VectorInvariantEnstrophyConserving())

@show surface_wind_stress_parameters = (τ₀ = 1e-4,
                                        Lφ = grid.Ly,
                                        φ₀ = 15)

surface_wind_stress(λ, φ, t, p) = p.τ₀ * cos(2π * (φ - p.φ₀) / p.Lφ)

surface_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress,
                                               parameters = surface_wind_stress_parameters)

μ = 1 / 60days

@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]

u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag,
                                         discrete_form = true,
                                         parameters = μ)

v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag,
                                         discrete_form = true,
                                         parameters = μ)

u_bcs = FieldBoundaryConditions(top = surface_wind_stress_bc,
                                bottom = u_bottom_drag_bc)

v_bcs = FieldBoundaryConditions(bottom = v_bottom_drag_bc)

@show const νh₀ = 5e3 * (60 / grid.Nx)^2

@inline νh(λ, φ, z, t) = νh₀ * cos(π * φ / 180)

variable_horizontal_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh)
constant_horizontal_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh₀)

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    architecture = CPU(),
                                    momentum_advection = VectorInvariant(),
                                    free_surface = free_surface,
                                    coriolis = coriolis,
                                    boundary_conditions = (u=u_bcs, v=v_bcs),
                                    closure = constant_horizontal_diffusivity,
                                    #closure = variable_horizontal_diffusivity,
                                    tracers = nothing,
                                    buoyancy = nothing)

g = model.free_surface.gravitational_acceleration

gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

# Time-scale for gravity wave propagation across the smallest grid cell
wave_propagation_time_scale = min(grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ)) * deg2rad(grid.Δλ),
                                  grid.radius * deg2rad(grid.Δφ)) / gravity_wave_speed

mutable struct Progress
    interval_start_time :: Float64
end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9

    @info @sprintf("Time: %s, iteration: %d, max(u): %.2e m s⁻¹, wall time: %s",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   maximum(sim.model.velocities.u),
                   prettytime(wall_time))

    p.interval_start_time = time_ns()

    return nothing
end

simulation = Simulation(model,
                        Δt = 3600,
                        stop_time = 1years,
                        iteration_interval = 100,
                        progress = Progress(time_ns()))

output_fields = merge(model.velocities, (η=model.free_surface.η,))

output_prefix = "barotropic_gyre_Nx$(grid.Nx)_Ny$(grid.Ny)"

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = TimeInterval(10day),
                                                      prefix = output_prefix,
                                                      field_slicer = nothing,
                                                      force = true)

run!(simulation)

#####
##### Animation!
#####

include("visualize_barotropic_gyre.jl")

# visualize_barotropic_gyre(simulation.output_writers[:fields])


text:
println("[TEST] Protobuf quick check on converting types and packing data")

module MyMod

  export MyType, MyTemplType, thing

  type MyType
    x::Array{Float64,1}
    val::Int64
    str::ASCIIString
    MyType() = new()
    MyType(x...) = new(x[1],x[2],x[3])
  end

  type MyTemplType{T}
    x::Array{Float64,1}
    val::Int64
    str::ASCIIString
    a::T
    MyTemplType() = new()
    MyTemplType(x...) = new(x[1],x[2],x[3],x[4])
  end

  function thing()
    println("Thing is running")
    nothing
  end

end

module ProtoTest
  using ProtoBuf
  export protostring, nothingspecialprotoread, registerType, runOutsideFunction, specialprotoread

  regTypes = Dict{ASCIIString, Type}()

  function protostring(datas)
    iob = PipeBuffer()
    @show typeof(datas)
    writeproto(iob, datas)
    return iob
  end

  function nothingspecialprotoread(iob, instance)
    return readproto(iob, instance)
  end

  typeConterters = Dict{AbstractString,Function}()
  function registerType(str::ASCIIString, spType::Type; converter::Union{Function,Union}=Union{})
    ProtoTest.regTypes[str] = spType
    if converter!=Union{}
      typeConverters[str] = converter
    end
    nothing
  end

  function specialprotoread(iob, instanceType::Type, Tstr::ASCIIString)
    it = instanceType{ ProtoTest.regTypes[Tstr] }()
    nothingspecialprotoread(iob, it)
  end

  function runOutsideFunction(f::Function)
    f()
    println("its done")
    nothing
  end

end

using MyMod, ProtoTest

mydatas = MyTemplType{MyType}(rand(2),1,"test01",MyType(rand(3),2, "testT"))

iob = protostring(mydatas)

registerType("MyType", MyType)
dd = specialprotoread(iob, MyTemplType, "MyType")

@show dd
@show mydatas

@test dd.x[1] == mydatas.x[1]

# dd = specialprotoread(iob, MyTemplType, MyType)

# dd = nothingspecialprotoread(iob, MyTemplType{MyType}())

# runOutsideFunction(MyMod.thing)

#


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_events_data
using AWS.Compat
using AWS.UUIDs

"""
    batch_acknowledge_alarm(acknowledge_action_requests)
    batch_acknowledge_alarm(acknowledge_action_requests, params::Dict{String,<:Any})

Acknowledges one or more alarms. The alarms change to the ACKNOWLEDGED state after you
acknowledge them.

# Arguments
- `acknowledge_action_requests`: The list of acknowledge action requests. You can specify
  up to 10 requests per operation.

"""
function batch_acknowledge_alarm(
    acknowledgeActionRequests; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "POST",
        "/alarms/acknowledge",
        Dict{String,Any}("acknowledgeActionRequests" => acknowledgeActionRequests);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_acknowledge_alarm(
    acknowledgeActionRequests,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/alarms/acknowledge",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("acknowledgeActionRequests" => acknowledgeActionRequests),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    batch_disable_alarm(disable_action_requests)
    batch_disable_alarm(disable_action_requests, params::Dict{String,<:Any})

Disables one or more alarms. The alarms change to the DISABLED state after you disable them.

# Arguments
- `disable_action_requests`: The list of disable action requests. You can specify up to 10
  requests per operation.

"""
function batch_disable_alarm(
    disableActionRequests; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "POST",
        "/alarms/disable",
        Dict{String,Any}("disableActionRequests" => disableActionRequests);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_disable_alarm(
    disableActionRequests,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/alarms/disable",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("disableActionRequests" => disableActionRequests),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    batch_enable_alarm(enable_action_requests)
    batch_enable_alarm(enable_action_requests, params::Dict{String,<:Any})

Enables one or more alarms. The alarms change to the NORMAL state after you enable them.

# Arguments
- `enable_action_requests`: The list of enable action requests. You can specify up to 10
  requests per operation.

"""
function batch_enable_alarm(
    enableActionRequests; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "POST",
        "/alarms/enable",
        Dict{String,Any}("enableActionRequests" => enableActionRequests);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_enable_alarm(
    enableActionRequests,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/alarms/enable",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("enableActionRequests" => enableActionRequests),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    batch_put_message(messages)
    batch_put_message(messages, params::Dict{String,<:Any})

Sends a set of messages to the AWS IoT Events system. Each message payload is transformed
into the input you specify (\"inputName\") and ingested into any detectors that monitor
that input. If multiple messages are sent, the order in which the messages are processed
isn't guaranteed. To guarantee ordering, you must send messages one at a time and wait for
a successful response.

# Arguments
- `messages`: The list of messages to send. Each message has the following format: '{
  \"messageId\": \"string\", \"inputName\": \"string\", \"payload\": \"string\"}'

"""
function batch_put_message(messages; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_events_data(
        "POST",
        "/inputs/messages",
        Dict{String,Any}("messages" => messages);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_put_message(
    messages,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/inputs/messages",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("messages" => messages), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    batch_reset_alarm(reset_action_requests)
    batch_reset_alarm(reset_action_requests, params::Dict{String,<:Any})

Resets one or more alarms. The alarms return to the NORMAL state after you reset them.

# Arguments
- `reset_action_requests`: The list of reset action requests. You can specify up to 10
  requests per operation.

"""
function batch_reset_alarm(
    resetActionRequests; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "POST",
        "/alarms/reset",
        Dict{String,Any}("resetActionRequests" => resetActionRequests);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_reset_alarm(
    resetActionRequests,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/alarms/reset",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("resetActionRequests" => resetActionRequests),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    batch_snooze_alarm(snooze_action_requests)
    batch_snooze_alarm(snooze_action_requests, params::Dict{String,<:Any})

Changes one or more alarms to the snooze mode. The alarms change to the SNOOZE_DISABLED
state after you set them to the snooze mode.

# Arguments
- `snooze_action_requests`: The list of snooze action requests. You can specify up to 10
  requests per operation.

"""
function batch_snooze_alarm(
    snoozeActionRequests; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "POST",
        "/alarms/snooze",
        Dict{String,Any}("snoozeActionRequests" => snoozeActionRequests);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_snooze_alarm(
    snoozeActionRequests,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/alarms/snooze",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("snoozeActionRequests" => snoozeActionRequests),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    batch_update_detector(detectors)
    batch_update_detector(detectors, params::Dict{String,<:Any})

Updates the state, variable values, and timer settings of one or more detectors (instances)
of a specified detector model.

# Arguments
- `detectors`: The list of detectors (instances) to update, along with the values to update.

"""
function batch_update_detector(detectors; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_events_data(
        "POST",
        "/detectors",
        Dict{String,Any}("detectors" => detectors);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_update_detector(
    detectors,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "POST",
        "/detectors",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("detectors" => detectors), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_alarm(alarm_model_name)
    describe_alarm(alarm_model_name, params::Dict{String,<:Any})

Retrieves information about an alarm.

# Arguments
- `alarm_model_name`: The name of the alarm model.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"keyValue"`: The value of the key used as a filter to select only the alarms associated
  with the key.
"""
function describe_alarm(alarmModelName; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_events_data(
        "GET",
        "/alarms/$(alarmModelName)/keyValues/";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_alarm(
    alarmModelName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "GET",
        "/alarms/$(alarmModelName)/keyValues/",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_detector(detector_model_name)
    describe_detector(detector_model_name, params::Dict{String,<:Any})

Returns information about the specified detector (instance).

# Arguments
- `detector_model_name`: The name of the detector model whose detectors (instances) you
  want information about.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"keyValue"`: A filter used to limit results to detectors (instances) created because of
  the given key ID.
"""
function describe_detector(
    detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "GET",
        "/detectors/$(detectorModelName)/keyValues/";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_detector(
    detectorModelName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "GET",
        "/detectors/$(detectorModelName)/keyValues/",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_alarms(alarm_model_name)
    list_alarms(alarm_model_name, params::Dict{String,<:Any})

Lists one or more alarms. The operation returns only the metadata associated with each
alarm.

# Arguments
- `alarm_model_name`: The name of the alarm model.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
"""
function list_alarms(alarmModelName; aws_config::AbstractAWSConfig=global_aws_config())
    return iot_events_data(
        "GET",
        "/alarms/$(alarmModelName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_alarms(
    alarmModelName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "GET",
        "/alarms/$(alarmModelName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_detectors(detector_model_name)
    list_detectors(detector_model_name, params::Dict{String,<:Any})

Lists detectors (the instances of a detector model).

# Arguments
- `detector_model_name`: The name of the detector model whose detectors (instances) are
  listed.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
- `"stateName"`: A filter that limits results to those detectors (instances) in the given
  state.
"""
function list_detectors(
    detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iot_events_data(
        "GET",
        "/detectors/$(detectorModelName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_detectors(
    detectorModelName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iot_events_data(
        "GET",
        "/detectors/$(detectorModelName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using BenchmarkTools
using SphereSurfaceHistogram
SSH = SphereSurfaceHistogram

# Needs to be defined, will be called on benchmark
const SUITE = BenchmarkGroup()

SUITE["Histogram creation"] = BenchmarkGroup()

suite = SUITE["Histogram creation"]["partition_sphere2"] = BenchmarkGroup()
for i in Int64[1e3, 1e4, 1e5, 1e6]
    suite["N = $i"] = @benchmarkable SSH.partition_sphere2($(4pi/i))
end

suite = SUITE["Histogram creation"]["partition_sphere_optim2"] = BenchmarkGroup()
for i in Int64[1e3, 1e4, 1e5, 1e6]
    suite["N = $i"] = @benchmarkable SSH.partition_sphere_optim2($(4pi/i))
end

suite = SUITE["Histogram creation"]["SSHBinner"] = BenchmarkGroup()
for i in Int64[1e3, 1e4, 1e5, 1e6]
    suite["N = $i"] = @benchmarkable SSHBinner($i)
end

suite = SUITE["Histogram creation"]["bin_many!"] = BenchmarkGroup()
for N_bins in Int64[1e2, 1e3, 1e4, 1e5]
    B = SSHBinner(N_bins)
    _suite = suite["N_bins = $N_bins"] = BenchmarkGroup()
    for N_spins in Int64[1e2, 1e3, 1e4, 1e5]
        spins = rand_spin(N_spins)

        _suite[
            "N_spins = $N_spins"
        ] = @benchmarkable bin_many!($B, $spins)
    end
end


text:
@testset "FWI" begin
	include("testGetData.jl")
	include("testSensitivity.jl")
	include("testInversion.jl")
	#include("testTimeDomain.jl")
end


text:
# G = (A*B) * (C*D)
@polly function kernel_3mm(E, A, B, F, C, D, G)
    ni,nk = size(A)
    nk,nj = size(B)
    nj,nm = size(C)
    nm,nl = size(D)

    # E := A*B
    for i = 1:ni, j = 1:nj
        E[i,j] = zero(eltype(E))
        for k = 1:nk
            E[i,j] += A[i,k] * B[k,j]
        end
    end

    # F := C*D
    for i = 1:nj, j = 1:nl
        F[i,j] = zero(eltype(F))
        for k = 1:nm
            F[i,j] += C[i,k] * D[k,j]
        end
    end

    # G := E*F
    for i = 1:ni, j = 1:nl
        G[i,j] = zero(eltype(G))
        for k = 1:nj
            G[i,j] += E[i,k] * F[k,j]
        end
    end
end

let
    ni = 800
    nj = 900
    nk = 1000
    nl = 1100
    nm = 1200

    E = zeros(Float32, ni, nj);
    A = zeros(Float32, ni, nk);
    B = zeros(Float32, nk, nj);
    F = zeros(Float32, nj, nl);
    C = zeros(Float32, nj, nm);
    D = zeros(Float32, nm, nl);
    G = zeros(Float32, ni, nl);

    for i = 1:ni, j = 1:nk
        A[i,j] = ((i*j+1) % ni) / (5*ni)
    end
    for i = 1:nk, j = 1:nj
        B[i,j] = ((i*(j+1)+2) % nj) / (5*nj)
    end
    for i = 1:nj, j = 1:nm
        C[i,j] = (i*(j+3) % nl) / (5*nl)
    end
    for i = 1:nm, j = 1:nl
        D[i,j] = ((i*(j+2)+2) % nk) / (5*nk)
    end

    SUITE["3mm"] = @benchmarkable kernel_3mm(E, A, B, F, C, D, G) setup = (E = copy($E); A = copy($A); B = copy($B); F = copy($F); C = copy($C); D = copy($D); G = copy($G))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 004a9900-576e-11eb-0a75-550da1a57e0d
begin
	using DrWatson

	@quickactivate "K-Ratio Project"
	using NeXLDatabase
	using Gadfly, DataFrames, PlutoUI, CSV
	db=openNeXLDatabase(datadir("exp_pro","kratio.db"))
	md"""
	# k-ratio Database Explorer
	Utilities for exploring the k-ratio database.
	"""
end

# ╔═╡ 6f8e06d0-576e-11eb-3dd5-9fc72c66f2c0

# ╔═╡ Cell order:
# ╠═004a9900-576e-11eb-0a75-550da1a57e0d
# ╠═6f8e06d0-576e-11eb-3dd5-9fc72c66f2c0


text:
function get_name_and_csv(path_to_filename)
    df = CSV.File(path_to_filename) |> DataFrames.DataFrame
    folder = splitdir(splitdir(path_to_filename)[1])[2]
    return folder, df
end

# Parser for Forecasts dat files
"""
Read all forecast CSV's in the path provided, the struct of the data should
follow this format
folder : PV
            file : DAY_AHEAD
            file : REAL_TIME
Folder name should be the device type
Files should only contain one real-time and day-ahead forecast
Args:
    files: A string
Returns:
    A dictionary with the CSV files as dataframes and folder names as keys
# TODO : Stochasti/Multiple scenarios
"""
 function read_data_files(rootpath::String; kwargs...)
   if :REGEX_FILE in keys(kwargs)
        REGEX_FILE = kwargs[:REGEX_FILE]
    else
        REGEX_FILE = r"(.*?)\.csv"
    end

    DATA = Dict{String, Any}()
    data = Dict{String, Any}()
    DATA["gen"] = data

    for (root, dirs, files) in walkdir(rootpath)

        for filename in files

            path_to_filename = joinpath(root, filename)
            if match(REGEX_FILE, path_to_filename) != nothing
                folder_name, csv_data = get_name_and_csv(path_to_filename)
                if folder_name == "load"
                    DATA["load"] = read_datetime(csv_data; kwargs...)
                else
                    data[folder_name] = read_datetime(csv_data; kwargs...)
                end
                @info "Successfully parsed $rootpath"
            else
                @warn "Unable to match regex with $path_to_filename"
            end
        end

    end

    return DATA
end

"""
Args:
    PowerSystems Dictionary
    Dictionary of all the data files
Returns:
    Returns an dictionary with Device name as key and PowerSystems Forecasts
    dictionary as values
"""
function assign_ts_data(ps_dict::Dict{String,Any},ts_dict::Dict{String,Any})
    if haskey(ts_dict,"load")
        ps_dict["load"] =  PowerSystems.add_time_series_load(ps_dict,ts_dict["load"])
    else
        @warn "Not assigning time series to loads"
    end

    if haskey(ts_dict,"gen")
        ts_map = _retrieve(ts_dict["gen"], Union{TimeSeries.TimeArray,DataFrames.DataFrame})
        for (key,val) in ts_map
            ts = _access(ts_dict["gen"],vcat(val,key))
            dd_map = _retrieve(ps_dict["gen"],key)
            if length(dd_map) > 0
                dd_map = string.(unique(values(dd_map))[1])
            else
                @warn("no $key entries in psdict")
                continue
            end
            dd = _access(ps_dict["gen"],vcat(dd_map,key))
            dd = PowerSystems.add_time_series(dd,ts)
        end
    else
        @warn "Not assigning time series to gens"
    end

    return ps_dict
end

function make_device_forecast(device::D, df::DataFrames.DataFrame, resolution::Dates.Period,horizon::Int) where {D<:Device}
    time_delta = Dates.Minute(df[2,:DateTime]-df[1,:DateTime])
    initialtime = df[1,:DateTime] # TODO :read the correct date/time when that was issued  forecast
    last_date = df[end,:DateTime]
    ts_dict = Dict{Any,Dict{Int,TimeSeries.TimeArray}}()
    ts_raw =  TimeSeries.TimeArray(df[1],df[2])
    for ts in initialtime:resolution:last_date
        ts_dict[ts] = Dict{Int,TimeSeries.TimeArray}(1 => ts_raw[ts:time_delta:(ts+resolution)])
    end
    forecast = Dict{String,Any}("horizon" =>horizon,
                                            "resolution" => resolution, #TODO : fix type conversion to JSON
                                            "interval" => time_delta,   #TODO : fix type conversion to JSON
                                            "initialtime" => initialtime,
                                            "device" => device,
                                            "data" => ts_dict
                                            )
    return forecast
end

 # -Parse csv file to dict
"""
Args:
    Dictionary of all the data files
    Length of the forecast - Week()/Dates.Day()/Dates.Hour()
    Forecast horizon in hours - Int64
    Array of PowerSystems devices in the systems - Renewable Generators and
    Loads
Returns:
    Returns an dictionary with Device name as key and PowerSystems Forecasts
    dictionary as values
"""
function make_forecast_dict(time_series::Dict{String,Any},
                            resolution::Dates.Period, horizon::Int,
                            Devices::Array{Generator,1})
    forecast = Dict{String,Any}()
    for device in Devices
        for (key_df,df) in time_series
            if device.name in map(String, names(df))
                for name in map(String, names(df))
                    if name == device.name
                        forecast[device.name] = make_device_forecast(device, df[[:DateTime, Symbol(device.name)]], resolution, horizon)
                    end
                end
            end
        end
        if !haskey(forecast,device.name)
            @info "No forecast found for $(device.name) "
        end
    end
    return forecast
end

"""
Args:
    Dictionary of all the data files
    Length of the forecast - Week()/Dates.Day()/Dates.Hour()
    Forecast horizon in hours - Int64
    Array of PowerSystems devices in the systems- Loads
Returns:
    Returns an dictionary with Device name as key and PowerSystems Forecasts
    dictionary as values
"""
function make_forecast_dict(time_series::Dict{String,Any},
                            resolution::Dates.Period, horizon::Int,
                            Devices::Array{ElectricLoad,1})
    forecast = Dict{String,Any}()
    for device in Devices
        if haskey(time_series,"load")
            if device.bus.name in  map(String, names(time_series["load"]))
                df = time_series["load"][[:DateTime,Symbol(device.bus.name)]]
                forecast[device.name] = make_device_forecast(device, df, resolution, horizon)
            end
        else
            @warn "No forecast found for Loads"
        end
    end
    return forecast
end

"""
Args:
    Dictionary of all the data files
    Length of the forecast - Week()/Dates.Day()/Dates.Hour()
    Forecast horizon in hours - Int64
    Array of PowerSystems devices in the systems- Loads
    Array of PowerSystems LoadZones

Returns:
    Returns an dictionary with Device name as key and PowerSystems Forecasts
    dictionary as values
"""
function make_forecast_dict(time_series::Dict{String,Any},
                            resolution::Dates.Period, horizon::Int,
                            Devices::Array{ElectricLoad,1},
                            LoadZones::Array{Device,1})
    forecast = Dict{String,Any}()
    for device in Devices
        if haskey(time_series,"load")
            for lz in LoadZones
                if device.bus in lz.buses
                    df = time_series["load"][[:DateTime,Symbol(lz.name)]]
                    forecast[device.name] = make_device_forecast(device, df, resolution, horizon)
                end
            end
        else
            @warn "No forecast found for Loads"
        end
    end
    return forecast
end

# - Parse Dict to Forecast Struct

"""
Args:
    A PowerSystems forecast dictionary
Returns:
    A PowerSystems forecast stuct array
"""
function make_forecast_array(dict)
    Forecasts = Array{Forecast}(undef, 0)
    for (device_key,device_dict) in dict
                push!(Forecasts,Deterministic(device_dict["device"],device_dict["horizon"],
                                device_dict["resolution"],device_dict["interval"],
                                device_dict["initialtime"],
                                device_dict["data"]
                                ))
    end
    return Forecasts
end

# Write dict to Json

function write_to_json(filename,Forecasts_dict)
    for (type_key,type_fc) in Forecasts_dict
        for (device_key,device_dicts) in type_fc
            stringdata =JSON.json(device_dicts, 3)
            open("$filename/$device_key.json", "w") do f
                write(f, stringdata)
             end
        end
    end
end

# Parse json to dict
#TODO : fix broken data formats
function parse_json(filename,device_names)
    Devices =Dict{String,Any}()
    for name in device_names
        if isfile("$filename/x$name.json")
            temp = Dict()
            open("$filename/x$name.json", "r") do f
            global temp
            dicttxt = readstring(f)  # file information to string
            temp=JSON.parse(dicttxt)  # parse and transform data
            Devices[name] = temp
            end
        end
    end
    return Devices
end


text:
@testset "GlobalCardinalityFixedOpen2Count: $(fct_type), $(array_size) items, $(sought_size) sought items, $(T)" for fct_type in ["vector of variables", "vector affine function"], array_size in [2, 3], sought_size in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(CountModel{T}())
    model = COIB.GlobalCardinalityFixedOpen2Count{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Count{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T},
    )

    x_counts, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:sought_size])
    sought_values = T(5) .+ collect(1:sought_size)

    if T == Int
        x_array, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:array_size])
    elseif T == Float64
        x_array = MOI.add_variables(model, array_size)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x_array..., x_counts...])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x_array..., x_counts...])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.GlobalCardinality{CP.OPEN_COUNTED_VALUES, T}(array_size, sought_values))

    for i in 1:array_size
        @test MOI.is_valid(model, x_array[i])
    end
    for i in 1:sought_size
        @test MOI.is_valid(model, x_counts[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}()) == sought_size

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}()) == bridge.cons_count
    end

    @testset "Sought array" begin
        @test length(bridge.cons_count) == sought_size
        for i in 1:sought_size
            @test MOI.is_valid(model, bridge.cons_count[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_count[i])
            @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_count[i]) == CP.Count(array_size, MOI.EqualTo(sought_values[i]))

            @test length(f.terms) == 1 + array_size
            for i in 1:(1 + array_size)
                @test f.terms[i].output_index == i
                @test f.terms[i].scalar_term.coefficient === one(T)
            end
            @test f.terms[1].scalar_term.variable === x_counts[i]
            for i in 1:array_size
                @test f.terms[1 + i].scalar_term.variable === x_array[i]
            end
        end
    end
end


text:
function dump_map(m::Map, filename::AbstractString)
    jsn = skeleton_dump(m)
    populate_routes!(jsn, m)

    json_dict = Dict{String,Any}("task_structure" => collect(values(jsn)))

    record_info!(json_dict, m)
    open(filename, "w") do f
        print(f, json(json_dict, 2))
    end
    return nothing
end

create_dict(d::Dict, k, v) = haskey(d, k) || (d[k] = v())

"""
    record_info!(json_dict, m::Map)

Record various info that may be helpful to the Project Manager.
"""
function record_info!(json_dict, m::Map)
    create_dict(json_dict, "info", Dict{String,Any})
    # Call the "check_routing" function and give it the "quiet" parameter to
    # keep get rid of redundant prints.
    json_dict["info"]["routing_success"] = Mapper2.check_routing(m; quiet = true)
end

################################################################################
struct RoutingTuple
    class           ::String
    source_task     ::_name_types
    source_index    ::_index_types
    dest_task       ::_name_types
    dest_index      ::_index_types
end

abstract type MapDump end
struct MapDumpNode <: MapDump
    name            ::String
    address         ::Tuple
    core_name       ::Union{String,Nothing}
    core_type       ::Union{String,Nothing}
    leaf_node_dict  ::Dict{String,Any}
end

# Constructor supplying empty endiding dictionary
function MapDumpNode(name, address::T, core_name, core_type) where T
    return MapDumpNode(
        name,
        T <: CartesianIndex ? address.I : address,
        core_name,
        core_type,
        Dict{String,Any}()
    )
end

struct MapDumpRoute <: MapDump
    network_id      ::Union{Int,Nothing}
    new_source_index::_index_types
    new_dest_index  ::_index_types
    source_task     ::_name_types
    source_index    ::_index_types
    dest_task       ::_name_types
    dest_index      ::_index_types
    offset_list     ::Vector
end

JSON.lower(m::MapDump) = Dict(string(f) => getfield(m,f) for f in fieldnames(typeof(m)))

"""
    skeleton_dump(m::Map)

Create entries in the final dictionary with just the task-name keys and the
address value populated.
"""
function skeleton_dump(m::Map)
    dict = Dict{String,MapDumpNode}()
    # Iterate through the dictionary in the mapping nodes
    for (name, address_path) in m.mapping.nodes
        node = getnode(m.taskgraph, name)

        addr = Mapper2.MapperCore.getaddress(m.toplevel,address_path)

        component = m.toplevel[address_path]
        core_name = component.metadata["pm_name"]
        core_type = component.metadata["pm_type"]

        dict[name] = MapDumpNode(name, addr, core_name, core_type)
    end
    return dict
end

function add_route!(n::MapDumpNode,
                    route,
                    class::String,
                    dir::String,
                    route_field::Symbol)

    leaf_dict = n.leaf_node_dict

    create_dict(leaf_dict, class, Dict{String,Any})
    class_dict = leaf_dict[class]

    create_dict(class_dict, dir, Dict{_index_types,Any})
    direction_dict = class_dict[dir]

    index = getfield(route, route_field)
    create_dict(direction_dict, index, Dict{String,Any})
    direction_dict[index] = route
end

function populate_routes!(jsn,m)
    routings = extract_routings(m)
    for (rtuple, route) in routings
        class = rtuple.class
        # Record the source index
        source_task = route.source_task
        source_node = jsn[source_task]

        add_route!(source_node, route, class, "output", :source_index)

        dest_task = route.dest_task
        dest_node = jsn[dest_task]
        add_route!(dest_node, route, class, "input", :dest_index)
    end
    return nothing
end

function extract_routings(m::Map)
    arch = m.toplevel
    routings = Dict{RoutingTuple,Any}()
    for (edge, graph) in zip(m.taskgraph.edges, m.mapping.edges)
        # Skip un-routed edges.
        needsrouting(rules(m), edge) || continue

        # Build the route tuple
        source_task = first(getsources(edge))
        dest_task   = first(getsinks(edge))
        offset_list = make_offset_list(arch, graph)

        source_index = edge.metadata["source_index"]
        dest_index   = edge.metadata["dest_index"]

        class = edge.metadata["pm_class"]

        # Get source and destination ports to extract;
        # - network_id
        # - new_source_index
        # - new_dest_index
        src_port_path = first(Mapper2.MapperGraphs.source_vertices(graph))
        dst_port_path = first(Mapper2.MapperGraphs.sink_vertices(graph))

        src_metadata = arch[src_port_path].metadata
        network_id = get(src_metadata,"network_id",nothing)
        new_source_index = src_metadata["index"]

        dst_metadata = arch[dst_port_path].metadata
        new_dest_index = dst_metadata["index"]

        index = edge.metadata["source_index"]

        key = RoutingTuple(class, source_task, source_index, dest_task, dest_index)
        routings[key] = MapDumpRoute(
            network_id,
            new_source_index,
            new_dest_index,
            source_task,
            source_index,
            dest_task,
            dest_index,
            offset_list
        )
    end
    return routings
end

function make_offset_list(arch, g)
    path = CartesianIndex{2}[]
    for p in Mapper2.MapperGraphs.linearize(g)
        # Skip global links since addresses they don't have an address
        isgloballink(p) && continue
        # Get the address from the path
        addr = Mapper2.MapperCore.getaddress(arch, p)
        if !in(addr,path)
            push!(path, addr)
        end
    end

    offset_list = [(path[i+1] - path[i]).I for i in 1:length(path)-1]
    return offset_list
end


text:
################################################################################
#
#   module LatPhysPlottingPyPlot
#   -> PyPlot
#   -> LatPhysBase
#   -> LinearAlgebra
#
#   --> PLOTTING OF LATTICES WITH PYPLOT
#           - 2D lattices
#           - 3D lattices
#
################################################################################

# start of module
module LatPhysReciprocalPlottingPyPlot

# used libraries
using LatPhysBase
using LatPhysReciprocal
using LinearAlgebra
using LatPhysPlottingPyPlot
using LaTeXStrings
using PyPlot

# explicitly import PyPlot.plot to overwrite
import PyPlot.plot
# explicitly import from basic plotting module
import LatPhysPlottingPyPlot.plotSite
import LatPhysPlottingPyPlot.plotSiteComplex

# Plotting of points
include("reciprocal_point_plotting.jl")

# Plotting of paths
include("reciprocal_path_plotting.jl")

# Plotting of Brillouin zones
include("brillouin_zone_plotting.jl")

# end of module
end


text:
# Types and constructors

struct QGradMonomialBasis{P,V,G,D} <: TensorPolynomialBasis{P,V,G}
  order::Int
  terms::CartesianIndices{D}
  perms::Matrix{Int}
end

function (::Type{QGradMonomialBasis{P,V}})(order::Int) where {P,V}
  D = _length(P)
  @assert _length(V) == D
  G = _gradient_type(V,P)
  n1d = order+1
  t = fill(n1d,D)
  t[1] = order
  terms = CartesianIndices(tuple(t...))
  perms = _prepare_perms(D)
  QGradMonomialBasis{P,V,G,D}(order,terms,perms)
end

# Implementation of the interface

length(b::QGradMonomialBasis{P,V,G,D}) where {P,V,G,D} = D*b.order*(b.order+1)^(D-1)

ndims(b::QGradMonomialBasis{P,V,G,D}) where {P,V,G,D} = D

function ScratchData(b::QGradMonomialBasis{P,V}) where {P,V}
  T = eltype(V)
  dim = _length(P)
  n1d = b.order+1
  c = zeros(dim,n1d)
  g = zeros(dim,n1d)
  MonomialBasisCache{T}(c,g)
end

function evaluate!(
  v::AbstractVector{V},
  b::QGradMonomialBasis{P,V,G,D},
  x::P,
  cache::MonomialBasisCache) where {P,V,G,D}
  _evaluate_nd_qgrad!(v,x,b.order,b.terms,b.perms,cache.c)
end

function gradient!(
  v::AbstractVector{G},
  b::QGradMonomialBasis{P,V,G,D},
  x::P,
  cache::MonomialBasisCache) where {P,V,G,D}
  _gradient_nd_qgrad!(v,x,b.order,b.terms,b.perms,cache.c,cache.g,V)
end

# Helpers

function _prepare_perms(D)
  perms = zeros(Int,D,D)
  for j in 1:D
    for d in j:D
      perms[d,j] =  d-j+1
    end
    for d in 1:(j-1)
      perms[d,j] =  d+(D-j)+1
    end
  end
  perms
end

function _evaluate_nd_qgrad!(
  v::AbstractVector{V},
  x,
  order,
  terms::CartesianIndices{D},
  perms::Matrix{Int},
  c::AbstractMatrix{T}) where {V,T,D}

  dim = D
  for d in 1:dim
    _evaluate_1d!(c,x,order,d)
  end

  o = one(T)
  k = 1
  m = zero(_mutable(V))
  js = eachindex(m)
  z = zero(T)

  for ci in terms

    for j in js

      @inbounds for i in js
        m[i] = z
      end

      s = o
      @inbounds for d in 1:dim
        s *= c[d,ci[perms[d,j]]]
      end

      m[j] = s
      v[k] = m
      k += 1

    end

  end

end

function _gradient_nd_qgrad!(
  v::AbstractVector{G},
  x,
  order,
  terms::CartesianIndices{D},
  perms::Matrix{Int},
  c::AbstractMatrix{T},
  g::AbstractMatrix{T},
  ::Type{V}) where {G,T,D,V}

  dim = D
  for d in 1:dim
    _evaluate_1d!(c,x,order,d)
    _gradient_1d!(g,x,order,d)
  end

  z = zero(_mutable(V))
  m = zero(_mutable(G))
  js = eachindex(z)
  mjs = eachindex(m)
  o = one(T)
  zi = zero(T)
  k = 1

  for ci in terms

    for j in js

      s = z
      for i in js
        s[i] = o
      end

      for q in 1:dim
        for d in 1:dim
          if d != q
            s[q] *= c[d,ci[perms[d,j]]]
          else
            s[q] *= g[d,ci[perms[d,j]]]
          end
        end
      end

      @inbounds for i in mjs
        m[i] = zi
      end

      for i in js
        m[i,j] = s[i]
      end
      v[k] = m
      k += 1

    end

  end

end