text:
using Test
using Quiqbox
using Quiqbox: isFull, BasisFuncMix, unpackBasisFuncs, inSymbols, varVal, ElementNames,
               sortBasisFuncs, ParamList, sumOf, expressionOfCore, mergeGaussFuncs,
               ijkOrbitalList, gaussProd, normOfGTOin
using Symbolics
using LinearAlgebra

@testset "Basis.jl" begin

# struct GaussFunc
exp1, con1 = (2, 1.0)
gf1 = GaussFunc(exp1, con1)
gf2 = GaussFunc(exp1, con1)
@test typeof(gf1) == typeof(gf2)
@test (gf1, gf2) isa Tuple{GaussFunc, GaussFunc}
@test (gf1, gf2) isa NTuple{2, GaussFunc}
gf1_ref = gf1
@test gf1.xpn == ParamBox(exp1, :xpn)
@test gf1.con == ParamBox(con1, :con)
@test gf1.param == (gf1.xpn, gf1.con)
@test hasEqual(gf1, gf2)
@test hasIdentical(gf1, gf1)
@test !hasIdentical(gf1, gf2)
@test gf1 !== gf2
@test (gf1 === gf1_ref)
@test hasIdentical(gf1, gf1_ref)
gf1_ref.xpn[] = 3.0
@test gf1.xpn[] == gf1.param[1][] == 3.0
gf1.param[2][] = 4
@test gf1_ref.con[] == gf1_ref.param[2][] == 4
@test gf1_ref === gf1
@test gf2 != gf1
x = [2.0, 1]
for i=1:2 gf1.param[i][] = x[i] end
@test gf1 === gf1_ref !== gf2
@test hasEqual(gf1, gf2)

# function genExponent genContraction
m1 = x->x^2
v1 = round(rand(), digits=6)
e1 = genExponent(v1, m1)
@test typeof(e1).parameters[2] == ParamList[:xpn]
@test e1[] == v1
@test e1() == m1(v1)

c1 = genContraction(v1, m1)
@test typeof(c1).parameters[2] == ParamList[:con]
@test c1[] == v1
@test e1() == m1(v1)

e2 = genExponent(e1)
e3 = genExponent(c1)
@test hasIdentical(e1, e2)
@test !hasEqual(e1, e3)
@test e1[] == e3[]
@test e1() == e3()
@test e1.map !== e3.map
v2 = rand()
@test e1.map(v2) == e3.map(v2)

c2 = genExponent(c1)
c3 = genExponent(e1)
@test hasIdentical(c1, c2)
@test !hasEqual(c1, c3)
@test c1[] == c3[]
@test c1() == c3()
@test c1.map !== c3.map
@test c1.map(v2) == c3.map(v2)

# function genSpatialPoint
k = fill(1.0)
v1 = genSpatialPoint([k,2,3  ])
v2 = genSpatialPoint([k,2,3.0])
@test v1 == v2
@test hasEqual(v1, v2)
@test !hasIdentical(v1, v2)
@test hasIdentical(v1[1], v2[1])
@test !hasIdentical(v1[2:3], v2[2:3])

# struct BasisFunc and function genBasisFunc
cen = [1,2,3]
cenParNames = [ParamList[:X], ParamList[:Y], ParamList[:Z]]
cenPar = ParamBox.(cen, cenParNames) |> Tuple
bf1 = genBasisFunc(cen, [gf1])
bf1_1 = genBasisFunc(cen, [gf1])
bf1_2 = genBasisFunc(bf1)
bf1_3 = genBasisFunc(cenPar, gf1)
bf1_4 = genBasisFunc(cenPar, gf1)
bf2 = genBasisFunc(cen, [gf2])
@test bf1 !== bf2
@test hasEqual(bf1, bf1_1, bf1_2, bf1_3, bf1_4, bf2)
@test hasIdentical(bf1, bf1_2)
@test !hasIdentical(bf1, bf1_1)
@test !hasIdentical(bf1, bf2)
@test !hasIdentical(bf1, bf1_3)
@test hasIdentical(bf1_3, bf1_4)
@test bf1.subshell == "S"
@test bf1 isa BasisFunc

bf11 = genBasisFunc(cen, [gf1, gf1])
bf1_3 = genBasisFunc(cen, (exp1, con1))
@test hasEqual(bf1, bf1_3)
@test !hasIdentical(bf1, bf1_3)
bf11_2 = genBasisFunc(cen, ([exp1, exp1], [con1, con1]))
@test hasEqual(bf11, bf11_2)
@test !hasIdentical(bf11, bf11_2)
bf2_P_norm2 = genBasisFunc(cen, [gf2], "P")
@test bf2_P_norm2.subshell == "P"
@test bf2_P_norm2 isa BasisFuncs

bf3_1 = genBasisFunc([0,0,0], (3,1))
bf3_2 = genBasisFunc([1,0,0], (2,1))
bf3_3 = genBasisFunc([0,0,0], (2,1), "P")
bf3_4 = genBasisFunc([0,0,0], (1,1), "D")
@test genBasisFunc([bf3_1, bf3_4, bf3_3, bf3_2]) == [bf3_1, bf3_3, bf3_4, bf3_2]

bf4_1 = genBasisFunc([0,0,0], "STO-3G")[]
bf4_2 = genBasisFunc([0,0,0], ("STO-3G", "He"))[]
bf4_3 = genBasisFunc([0,0,0], "STO-3G", nucleus = "He")[]
bf4s1 = genBasisFunc([0,0,0], ["STO-3G", "STO-3G"])
bf4s2 = genBasisFunc([0,0,0], ["STO-3G", "STO-3G"], nucleus = "He")
bf4s3 = genBasisFunc([0,0,0], [("STO-3G", "H"), ("STO-3G", "He")])
@test hasEqual.(Ref(bf4_1), bf4s1, Ref(bf4s3[1])) |> prod
@test hasEqual.(Ref(bf4_2), Ref(bf4_3), bf4s2, Ref(bf4s3[2])) |> prod
@test isapprox(1, overlap(bf4_1, bf4_1)[], atol=1e-8)
@test isapprox(1, overlap(bf4_2, bf4_2)[], atol=1e-8)

errorThreshold1 = 1e-11
bf3_2_2 = genBasisFunc([1,0,0], (2,1), normalizeGTO=true)
@test isapprox(1, overlap(bf3_2_2, bf3_2_2)[], atol=errorThreshold1)
@test isapprox(0.0553891828418, overlap(bf3_2, bf3_2)[], atol=errorThreshold1)
bf3_3_2 = genBasisFunc([0,0,0], (2,1), "P", normalizeGTO=true)
@test isapprox(LinearAlgebra.I, overlap(bf3_3_2, bf3_3_2), atol=errorThreshold1)
@test isapprox(0.0207709435653*LinearAlgebra.I, overlap(bf3_3, bf3_3), atol=errorThreshold1)
bf3_4_2 = genBasisFunc([0,0,0], (1,1), "D", normalizeGTO=true)
@test isapprox([1.0 0.0 0.0 1/3 0.0 1/3;
                0.0 1/3 0.0 0.0 0.0 0.0;
                0.0 0.0 1/3 0.0 0.0 0.0;
                1/3 0.0 0.0 1.0 0.0 1/3;
                0.0 0.0 0.0 0.0 1/3 0.0;
                1/3 0.0 0.0 1/3 0.0 1.0],
                overlap(bf3_4_2, bf3_4_2), atol=errorThreshold1)
@test isapprox([0.3691314831028692 0.0 0.0 0.1230438277009564 0.0 0.1230438277009564;
                0.0 0.1230438277009564 0.0 0.0 0.0 0.0;
                0.0 0.0 0.1230438277009564 0.0 0.0 0.0;
                0.1230438277009564 0.0 0.0 0.3691314831028692 0.0 0.1230438277009564;
                0.0 0.0 0.0 0.0 0.1230438277009564 0.0;
                0.1230438277009564 0.0 0.0 0.1230438277009564 0.0 0.3691314831028692],
                overlap(bf3_4, bf3_4), atol=errorThreshold1)

# function sortBasisFuncs
bfs1 = [genBasisFunc([1,1,1], (2,1), [1,0,0]), genBasisFunc([1,1,1], (3,1), [2,0,0]),
        genBasisFunc([1,1,2], (3,1), [0,0,0]), genBasisFunc([1,1,1], (3,1), "P")]
bfs2 = [genBasisFunc([1,1,1], (2,1), [1,0,0]), genBasisFunc([1,1,1], (3,1), "P"),
        genBasisFunc([1,1,1], (3,1), [2,0,0]), genBasisFunc([1,1,2], (3,1), [0,0,0])]
bfs3 = sortBasisFuncs(bfs1, groupCenters=true)
@test length.(bfs3) == [3,1]
bfs3 = bfs3 |> flatten
@test !hasEqual(bfs1, bfs2)
@test  hasEqual(bfs3, bfs2)

# function isFull
@test isFull(bf3_1) == true
@test isFull(bf3_3) == true
bf5 = genBasisFunc([0,0,0], (2,1), [1,0,0])
@test isFull(bf5) == false
@test isFull(1) == false

# function centerOf centerCoordOf
@test centerOf(bf5) == bf5.center
@test centerCoordOf(bf5) == [0,0,0]

# struct BasisFuncMix
bfm1 = BasisFuncMix(bf1)
@test bfm1 == BasisFuncMix([bf1])
@test BasisFuncMix(BasisFuncMix(bf3_1)) == BasisFuncMix(bf3_1)
bf5_2 = genBasisFunc([0,0,0], (2,1), [[1,0,0]])
bfm2 = BasisFuncMix(bf5)
@test hasEqual(bfm2, BasisFuncMix(bf5_2)[])

errorThreshold2 = 1e-15
bs1 = genBasisFunc.(gridCoords(GridBox(1,1.5)), Ref(GaussFunc(1.0, 0.5)))
nuc = ["H", "H"]
nucCoords = [rand(3), rand(3)]
bfm = BasisFuncMix(bs1)
S = overlaps([bfm])[]
@test S == overlap(bfm, bfm)[]
@test isapprox(S, overlaps(bs1) |> sum, atol=errorThreshold2)
T = elecKinetics([bfm])[]
@test T == elecKinetic(bfm, bfm)[]
@test isapprox(T, elecKinetics(bs1) |> sum, atol=errorThreshold2)
V = nucAttractions([bfm], nuc, nucCoords)[]
@test V == nucAttraction(bfm, bfm, nuc, nucCoords)[]
@test isapprox(V, nucAttractions(bs1, nuc, nucCoords) |> sum, atol=errorThreshold2)
eeI = eeInteractions([bfm])[]
@test eeI == eeInteraction(bfm, bfm, bfm, bfm)[]
@test isapprox(eeI, eeInteractions(bs1) |> sum, atol=errorThreshold2)

# function sumOf
bs2 = [genBasisFunc([1,1,1], (2,1), [1,0,0], normalizeGTO=true),
       genBasisFunc([1,1,1], (3,1), [2,0,0], normalizeGTO=true),
       genBasisFunc([1,1,2], (3,1), [0,0,0], normalizeGTO=true),
       genBasisFunc([1,1,1], (3,1), [0,1,0], normalizeGTO=true)]
bs2_2 = [genBasisFunc([1,1,1], (2,1), [1,0,0]),
         genBasisFunc([1,1,1], (3,1), [2,0,0]),
         genBasisFunc([1,1,2], (3,1), [0,0,0]),
         genBasisFunc([1,1,1], (3,1), [0,1,0])]
bs2_3 = [genBasisFunc([1,1,1], (2,1), [1,0,0]),
         genBasisFunc([1,1,1], (3,1), [2,0,0], normalizeGTO=true),
         genBasisFunc([1,1,2], (3,1), [0,0,0]),
         genBasisFunc([1,1,1], (3,1), [0,1,0], normalizeGTO=true)]
bs3 = [genBasisFunc([1,1,1], (2,1), [1,0,0], normalizeGTO=true),
       genBasisFunc([1,1,1], (3,1), [0,1,0], normalizeGTO=true),
       genBasisFunc([1,1,1], (3,1), [2,0,0], normalizeGTO=true),
       genBasisFunc([1,1,2], (3,1), [0,0,0], normalizeGTO=true)]
bs3_2 = [genBasisFunc([1,1,1], (2,1), [1,0,0]),
       genBasisFunc([1,1,1], (3,1), [0,1,0]),
       genBasisFunc([1,1,1], (3,1), [2,0,0]),
       genBasisFunc([1,1,2], (3,1), [0,0,0])]

bfm_1 = +(bs2...,)
bfm_2 = sumOf(bs2)
bfm_3 = BasisFuncMix(bs3)
bfm_4 = +(bs2_2...,)
bfm_5 = sumOf(bs2_2)
bfm_6 = BasisFuncMix(bs3_2)
bfm_7 = sumOf([bfm_6])
@test hasEqual(bfm_1, bfm_3)
@test hasEqual(bfm_2, bfm_3)
@test hasEqual(bfm_4, bfm_6)
@test hasEqual(bfm_5, bfm_6)
@test hasEqual(bfm_6, bfm_7)

# function mergeGaussFuncs
gf_merge1 = GaussFunc(2,1)
gf_merge2 = GaussFunc(2,1)
gf_merge3 = GaussFunc(2,1)

@test mergeGaussFuncs(gf_merge1) === gf_merge1

mgf1 = mergeGaussFuncs(gf_merge1, gf_merge1)[]
mgf2 = mergeGaussFuncs(gf_merge1, gf_merge2)[]
@test mgf1.xpn() == 2 == mgf1.con()
@test !hasIdentical(mgf1, mgf2)
gf_merge1_2 = GaussFunc(gf_merge1.xpn, gf_merge2.con)
gf_merge1_3 = GaussFunc(gf_merge2.xpn, gf_merge1.con)
mgf1_2 = mergeGaussFuncs(gf_merge1, gf_merge1_2)[]
mgf1_3 = mergeGaussFuncs(gf_merge1, gf_merge1_3)[]
@test mgf1_2.xpn() == 2 == mgf1_2.con()
@test !hasIdentical(mgf1, mgf1_2)
@test !hasIdentical(mgf1_2, mgf1_3)
@test hasEqual(mgf1, mgf1_2, mgf1_3)

gf_merge3 = GaussFunc(1.5,1)
@test hasIdentical(mergeGaussFuncs(gf_merge1, gf_merge3), [gf_merge1, gf_merge3])

# function add, mul, gaussProd
@test add(bs2[1]) === bs2[1]
bf1s = BasisFuncs(bf1.center, bf1.gauss, [ijkOrbitalList[bf1.ijk[1]]], bf1.normalizeGTO)
@test hasIdentical(add(bf1s), bf1)
bfm_add1 = BasisFuncMix(bs2)
@test add(bfm_add1) == sumOf(bs2)
bf_add1 = genBasisFunc([1,2,1], (2,1.1))
bf_add2 = genBasisFunc([1,1,1], (1,1.2))
bf_add1_2 = BasisFuncMix(bf_add1)
bf_add2_2 = BasisFuncMix(bf_add2)
@test hasEqual(add(bf_add1, bf_add2), add(bf_add2, bf_add1))
@test hasEqual(add(bf_add1_2, bf_add2_2), add(bf_add2_2, bf_add1_2))
@test hasEqual(add(bf_add1_2, bf_add2), add(bf_add1, bf_add2))
@test hasEqual(add(bf_add1, bf_add2_2), add(bf_add1, bf_add2))
@test hasEqual(add(bf_add1_2, bf_add2_2), add(bf_add1, bf_add2))

for bs in (bs2, bs2_2, bs2_3)
    X = overlaps(bs)^(-0.5)
    bsNew = [mul.(bs, @view X[:,i]) for i in 1:size(X, 2)] .|> sum
    SNew = overlaps(bsNew)
    @test isapprox(SNew, LinearAlgebra.I, atol=1e-14)
end
@test hasEqual(mul(1.0, bf1), bf1)
@test hasEqual(mul(1.1, bf1), mul(bf1, 1.1))

α₁, α₂ = rand(1:0.01:10, 2)
d₁, d₂ = rand(2)
R₁ = rand(-2:0.01:2, 3)
R₂ = rand(-2:0.01:2, 3)
xr = -10:0.1:10
yr = -10:0.1:10
zr = -10:0.1:10
gf_mul = mul(GaussFunc(α₁, d₁), GaussFunc(α₂, d₂))
@test isapprox(gf_mul.xpn[], α₁ + α₂, atol=1e-14)
@test isapprox(gf_mul.con[], d₁ * d₂, atol=1e-14)
bl = true
for x in xr, y in yr, z in zr
    xv = [x, y, z]
    α₃, d₃, R₃ = gaussProd((α₁, d₁, R₁), (α₂, d₂, R₂))
    bl *= isapprox(d₁*exp(-α₁*sum(abs2, xv-R₁)) * d₂*exp(-α₂*sum(abs2, xv-R₂)),
                   d₃*exp(-α₃*sum(abs2, xv-R₃)), atol=1e-16) # limit: atol=1e-20
end
@test bl

@test hasEqual(mul(bf1s, bf1s), [mul(bf1s[1], bf1s[1])])
@test hasEqual(mul(bf1s, bf1s, normalizeGTO=true)[],
               mul(bf1s[1], bf1s[1], normalizeGTO=true))
bf_mul1 = genBasisFunc([1,0,0], (2,3))
bf_mul2 = genBasisFunc([1,0,0], (1.5,1))
bf_mul2_2 = genBasisFunc([1,0,0], (1.5,1), normalizeGTO=true)
@test hasEqual(mul(bf_mul1, bf_mul2),  mul(BasisFuncMix([bf_mul1]), BasisFuncMix([bf_mul2])))
@test hasEqual(mul(bf_mul1, bf_mul2, normalizeGTO=true),
               mul(BasisFuncMix([bf_mul1]), BasisFuncMix([bf_mul2]), normalizeGTO=true))
bf_mul3 = genBasisFunc([1,0,0], ([1.0, 1.5], [0.3, 0.4]))
bf_mul4 = mul(bf_mul1, bf_mul2)
bf_mul4_0 = genBasisFunc([1,0,0], (3.5,3))
@test hasEqual(bf_mul4, bf_mul4_0)
bf_mul5 = mul(bf_mul1, bf_mul2, normalizeGTO=true)
bf_mul5_0 = genBasisFunc([1,0,0], (3.5,3), normalizeGTO=true)
@test hasEqual(bf_mul5, bf_mul5_0)
bf_mul6 = mul(bf_mul1, bf_mul2_2)
bf_mul6_0 = genBasisFunc([1,0,0], (3.5,3*normOfGTOin(bf_mul2)[]))
@test hasEqual(bf_mul6, bf_mul6_0)
bf_mul7 = mul(bf_mul1, bf_mul2_2, normalizeGTO=true)
bf_mul7_0 = genBasisFunc([1,0,0], (3.5,3*normOfGTOin(bf_mul2)[]), normalizeGTO=true)
@test hasEqual(bf_mul7, bf_mul7_0)
bf_mul8 = mul(bf_mul1, bf_mul3)
bf_mul8_0 = genBasisFunc([1,0,0], ([3, 3.5], [0.9, 1.2]))
@test hasEqual(bf_mul8, bf_mul8_0)
bf_mul9 = mul(bf_mul1, bf_mul3, normalizeGTO=true)
bf_mul9_0 = genBasisFunc([1,0,0], ([3, 3.5], [0.9, 1.2]), normalizeGTO=true)
@test hasEqual(bf_mul9, bf_mul9_0)
bf_mul10 = genBasisFunc([1,0,0], (1.2,0.3), [1,0,0])
bf_mul11 = genBasisFunc([0,1,0], (1.5,1))
bf_mul12 = genBasisFunc([0,1,0], (1.5,1), [1,0,0])
bf_mul13 = genBasisFunc([1,1,1], ([1.2, 1.0], [0.2, 0.5]))
bf_mul14 = genBasisFunc([1,0.2,1], ([1.2, 1.0], [0.2, 0.5]), [1,1,1])
bf_mul15 = mul(bf_mul10, BasisFuncMix([bf_mul1, bf_mul2]))
bf_mul15_0 = mul(bf_mul10, BasisFuncMix([bf_mul1, bf_mul2]), normalizeGTO=true)
bf_mul16 = mul(bf_mul15, bf_mul15_0)
bf_mul16_0 = mul(bf_mul15, bf_mul15_0, normalizeGTO=true)

testNorm = function (bf1, bf2)
    bf3 = mul(bf1, bf2)
    n1 = Quiqbox.getOverlap(bf3, bf3)
    n2 = Quiqbox.getOverlap(mul(bf1, bf1), mul(bf2, bf2))
    @test isapprox(n1, n2, atol=1e-12)
end

testNorm(bf_mul1, bf_mul10)
testNorm(bf_mul1, bf_mul11)
testNorm(bf_mul1, bf_mul12)
testNorm(bf_mul1, bf_mul13)
testNorm(bf_mul1, bf_mul14)
testNorm(bf_mul1, bf_mul15)
testNorm(bf_mul1, bf_mul15_0)
testNorm(bf_mul1, bf_mul16)
testNorm(bf_mul1, bf_mul16_0)
testNorm(bf_mul10, bf_mul11)
testNorm(bf_mul11, bf_mul13)
testNorm(bf_mul12, bf_mul13)
testNorm(bf_mul12, bf_mul14)
testNorm(bf_mul13, bf_mul14)
testNorm(bf_mul12, bf_mul15)
testNorm(bf_mul13, bf_mul15)
testNorm(bf_mul14, bf_mul15)
testNorm(bf_mul15, bf_mul15_0)
testNorm(bf_mul12, bf_mul16)
testNorm(bf_mul13, bf_mul16)
testNorm(bf_mul14, bf_mul16)
testNorm(bf_mul16, bf_mul16_0)

# (dev) function getOverlap
@test isapprox(Quiqbox.getOverlap(bf_mul2_2, bf_mul2_2), 1.0, atol=1e-10)
@test isapprox(Quiqbox.getOverlap(bf_mul5_0, bf_mul5_0), 9.0, atol=1e-10)
bf_mul18 = genBasisFunc(rand(3), "STO-3G")[]
@test isapprox(Quiqbox.getOverlap(bf_mul18, bf_mul18), 1.0, atol=1e-10)

# function shift
ijk = [1,0,0]
didjdk = [0,1,1]
bf_os1 = genBasisFunc([0,0,0], (2,1), ijk)
bf_os2 = genBasisFunc([0,0,0], (2,1), ijk, normalizeGTO=true)
bf_os1S = genBasisFunc([0,0,0], (2,1), ijk+didjdk)
bf_os2S = genBasisFunc([0,0,0], (2,1), ijk+didjdk, normalizeGTO=true)
@test hasEqual(shift(bf_os1, didjdk), bf_os1S)
@test hasEqual(shift(bf_os2, didjdk), bf_os2S)

# function unpackBasisFuncs
@test unpackBasisFuncs(bfm1)[1] == bf1
@test unpackBasisFuncs(bf1) == [bf1]
@test unpackBasisFuncs(0) == []

# function decompose
bf_d_1 = genBasisFunc([1,0,0], (1, 0.5))
bf_d_2 = genBasisFunc([1,0,0], ([2,1], [0.1, 0.5]))
bf_d_3 = genBasisFunc([1,1,0], ([2,1], [0.1, 0.2]), "P")
bm_d_1 = BasisFuncMix([bf_d_1, bf_d_2])
dm1 = reshape([bf_d_1], 1, 1)
@test hasEqual(decompose(bf_d_1), dm1)
@test hasEqual(decompose(bf_d_1, splitGaussFunc=true), dm1)
@test hasEqual(decompose(bf_d_2), reshape([bf_d_2], 1, 1))
dm2 = reshape([genBasisFunc([1,0,0], (2, 0.1)), bf_d_1], 2, 1)
@test hasEqual(decompose(bf_d_2, splitGaussFunc=true), dm2)
@test hasEqual(decompose(bf_d_3), reshape(bf_d_3[:], 1, 3))
@test hasEqual(decompose(bf_d_3, splitGaussFunc=true),
               hcat(decompose.(bf_d_3[:], splitGaussFunc=true)...))
@test decompose(bm_d_1) === bm_d_1
@test hasEqual(decompose(bm_d_1, splitGaussFunc=true),
              reshape(vcat(decompose.(bm_d_1.BasisFunc, splitGaussFunc=true)...), 3, 1))

# function basisSize
@test basisSize("P") == (3,)
@test basisSize(["S", "P", "D"]) == (1, 3, 6)
@test basisSize(bf1) == (1,)
@test basisSize(bfm1) == (1,) == basisSize(bfm2)
@test basisSize([bf1, bf2, bf3_3, bf5]) == (1,1,3,1)

# function genGaussFuncText
bfCoeff = [[6.163845031, 1.097161308], [0.4301284983, 0.6789135305],
           [0.245916322, 0.06237087296], [0.0494717692, 0.9637824081],
           [0.245916322, 0.06237087296], [0.5115407076, 0.6128198961]]
bfCoeff2 = vcat([[bfCoeff[2i+1]'; bfCoeff[2i+2]']' for i=0:2]...)
content = """
S    2   1.0
         6.163845031                   0.4301284983
         1.097161308                   0.6789135305
S    2   1.0
         0.245916322                   0.0494717692
         0.06237087296                 0.9637824081
P    2   1.0
         0.245916322                   0.5115407076
         0.06237087296                 0.6128198961
"""
lines = (content |> IOBuffer |> readlines)
@test map(i->Quiqbox.genGaussFuncText(bfCoeff2[i,:]...), 1:size(bfCoeff2)[1] |> collect) ==
      [lines[2], lines[3], lines[5], lines[6], lines[8], lines[9]].*"\n"

# function genBasisFuncText & genBFuncsFromText
randElement = ElementNames[rand(1:length(ElementNames))]
bs1 = genBasisFunc(missing, ("6-31G", "H"))
cens = [rand(3) for _=1:length(bs1)]
txt1 = genBasisFuncText(bs1, printCenter=false, groupCenters=false) |> join
txt2 = genBasisFuncText(bs1, printCenter=false) |> join
bs2_1 = genBFuncsFromText(txt1)
bs2_2 = genBFuncsFromText(txt2)
assignCenter!.(cens, bs1)
assignCenter!.(cens, bs2_1)
assignCenter!.(cens, bs2_2)
txt3 = genBasisFuncText(bs1) |> join
bs2_3 = genBFuncsFromText(txt3)
@test hasEqual(bs1, bs2_1, ignoreContainer=true)
@test hasEqual(bs1, bs2_2, ignoreContainer=true)
@test hasEqual.(sortBasisFuncs(bs1), bs2_3, ignoreFunction=true) |> prod
@test hasEqual.(bs1, bs2_1, ignoreFunction=true) |> prod
@test hasEqual.(bs1, bs2_2, ignoreFunction=true) |> prod
@test hasEqual.(sortBasisFuncs(bs1), bs2_3, ignoreFunction=true) |> prod

# function assignCenter!
bf6 = genBasisFunc(missing, "STO-3G")[]
coord = [1,0,0]
bf6_1 = genBasisFunc(coord, "STO-3G")[]
@test !hasEqual(bf6, bf6_1)
assignCenter!(coord, bf6)
@test hasEqual(bf6, bf6_1)

# function getParams
pb1 = ParamBox(2, :p, canDiff=false)
@test getParams(pb1) == pb1
@test getParams(pb1, :p) == pb1
@test getParams(pb1, :P) === nothing
@test getParams(pb1, onlyDifferentiable=true) === nothing

pb2 = ParamBox(2, :q)
@test getParams([pb1, pb2]) == [pb1, pb2]
@test getParams([pb1, pb2], :p) == [pb1]

gf_pbTest1 = GaussFunc(2,1)
@test getParams(gf_pbTest1) == gf_pbTest1.param |> collect

gf_pbTest2 = GaussFunc(1.5,0.5)
bf_pbTest1 = genBasisFunc([1,0,0], [gf_pbTest1, gf_pbTest2])
@test getParams(bf_pbTest1) == [bf_pbTest1.center...,
                                gf_pbTest1.param..., gf_pbTest2.param...]

alpha = ParamList[:xpn]
@test getParams(bf_pbTest1, alpha) ==
      vcat(getParams(gf_pbTest1, alpha), getParams(gf_pbTest2, alpha)) ==
      [gf_pbTest1.param[1], gf_pbTest2.param[1]]

cs = [pb1, bf_pbTest1]
ss = [:X, :Y, :Z]
@test [getParams(i, j) for i in cs, j in ss] |> flatten ==
      [nothing, bf_pbTest1.center[1], nothing, bf_pbTest1.center[2],
       nothing, bf_pbTest1.center[3]]
@test getParams(cs) == vcat(getParams(pb1), getParams(bf_pbTest1))
@test (getParams.(Ref(cs), ss) |> flatten) == (bf_pbTest1.center |> collect)

# function copyBasis
e = genExponent(3.0, x -> x^2 + 1)
c = genContraction(2.0)
gf_cb1 = GaussFunc(e, c)
gf_cb2 = copyBasis(gf_cb1)
@test !hasEqual(gf_cb1, gf_cb2)
@test hasEqual(gf_cb1.con, gf_cb2.con)
@test gf_cb1.xpn() == gf_cb2.xpn() == gf_cb2.xpn[] != gf_cb1.xpn[]
gf_cb3 = copyBasis(gf_cb1, false)
@test !hasEqual(gf_cb1, gf_cb3)
@test gf_cb1.xpn() != gf_cb3.xpn() == gf_cb3.xpn[] == gf_cb1.xpn[]
cen_cb1 = rand(3)

bf_cb1 = genBasisFunc(cen_cb1, [gf_cb1, gf_cb3])
bf_cb2 = copyBasis(bf_cb1)
bf_cb3 = copyBasis(bf_cb1, false)
testbf_cb = function (bf1, bf2)
    @test bf1.center == bf2.center
    @test bf1.ijk == bf2.ijk
    @test bf1.normalizeGTO == bf2.normalizeGTO
end
testbf_cb(bf_cb1, bf_cb2)
testbf_cb(bf_cb1, bf_cb3)
@test hasEqual(bf_cb2.gauss, (gf_cb2, gf_cb3))
@test hasEqual(bf_cb3.gauss, (gf_cb3, gf_cb3))

bfm_cb1 = BasisFuncMix([bf_cb1, bf_cb3])
bfm_cb2 = copyBasis(bfm_cb1)
bfm_cb3 = copyBasis(bfm_cb1, false)
testbf_cb.(bfm_cb1.BasisFunc, bfm_cb2.BasisFunc)
testbf_cb.(bfm_cb1.BasisFunc, bfm_cb3.BasisFunc)
@test hasEqual(bfm_cb2.BasisFunc |> collect, [bf_cb2, bf_cb3])
@test hasEqual(bfm_cb3.BasisFunc |> collect, [bf_cb3, bf_cb3])

# function uniqueParams! getVar getVarDict
e_gv1 = genExponent(2.0)
c_gv1 = genContraction(1.0)
gf_gv1 = GaussFunc(e_gv1, c_gv1)
e_gv2 = genExponent(2.5)
c_gv2 = genContraction(0.5)
gf_gv2 = GaussFunc(e_gv2, c_gv2)
e_gv3 = genExponent(1.05, x->x^2)
c_gv3 = genContraction(1.5)
gf_gv3 = GaussFunc(e_gv3, c_gv3)
gf_gv12 = GaussFunc(e_gv1, c_gv2)
gf_gv23 = GaussFunc(e_gv2, c_gv3)
gf_gv31 = GaussFunc(e_gv3, c_gv1)
x_gv1 = fill(1.0)
y_gv1 = fill(5.0)
cen_gv1 = genSpatialPoint([x_gv1, 2, 3])
cen_gv2 = genSpatialPoint([4, y_gv1, 6])
cen_gv3 = genSpatialPoint([x_gv1, y_gv1, 6])

bf_gv1 = genBasisFunc(cen_gv1, gf_gv1)
bf_gv2 = genBasisFunc(cen_gv2, gf_gv2)
bf_gv3 = genBasisFunc(cen_gv2, [gf_gv1, gf_gv2])
bf_gv4 = genBasisFunc(cen_gv1, gf_gv3)
bf_gv5 = genBasisFunc(cen_gv3, gf_gv1)
bf_gv6 = genBasisFunc(cen_gv3, [gf_gv12, gf_gv23, gf_gv31])
bfm_gv = BasisFuncMix([bf_gv2, bf_gv4, bf_gv6])

@test uniqueParams!(bf_gv4) == [cen_gv1..., e_gv3, c_gv3]

pbs_gv1 = uniqueParams!([bf_gv1, bf_gv2])
pbs_gv1_0 = [cen_gv1[1], cen_gv2[1],
             cen_gv1[2], cen_gv2[2],
             cen_gv1[3], cen_gv2[3],
                  e_gv1,      e_gv2,
                  c_gv1,      c_gv2]
@test pbs_gv1 == pbs_gv1_0
@test hasIdentical(pbs_gv1, pbs_gv1_0)

pbs_gv2 = uniqueParams!([bf_gv2, bf_gv1, bf_gv3])

pbs_gv2_0 = [cen_gv2[1], cen_gv1[1], cen_gv2[1],
             cen_gv2[2], cen_gv1[2], cen_gv2[2],
             cen_gv2[3], cen_gv1[3], cen_gv2[3],
             e_gv2, e_gv1, e_gv1, e_gv2,
             c_gv2, c_gv1, c_gv1, c_gv2]
@test pbs_gv2 == pbs_gv2_0
@test hasIdentical(pbs_gv2, pbs_gv2_0)

bfs_gv = sortBasisFuncs([bf_gv1, bf_gv2, bf_gv3, bf_gv4, bf_gv5, bf_gv6])
pbs_gv3 = uniqueParams!(bfs_gv, filterMapping=true)
pbs_gv4 = uniqueParams!(vcat(bfs_gv, bfm_gv), filterMapping=true)
@test pbs_gv3 == pbs_gv4
@test hasIdentical(pbs_gv3, pbs_gv4)

pbs_gv0 = [cen_gv1[1], cen_gv2[1],
           cen_gv1[2], cen_gv2[2],
           cen_gv1[3], cen_gv2[3], cen_gv3[3],
           e_gv1, e_gv2, e_gv3,
           c_gv1, c_gv2, c_gv3]
pbs_gv0_2 = sort(pbs_gv0, by=x->(typeof(x).parameters[2], x.index[]))
pbs_gv3_2 = sort(pbs_gv3, by=x->(typeof(x).parameters[2], x.index[]))
@test pbs_gv3_2 == pbs_gv0_2
@test hasEqual(pbs_gv3_2, pbs_gv0_2)
@test hasIdentical(pbs_gv3_2, pbs_gv0_2)

@test getVar(e_gv1).val.name == :α₁
@test getfield.(getfield.(getVar(bf_gv6), :val), :name) ==
      [:X₁, :Y₂, :Z₂, :α₁, :d₃, :α₃, :d₂, :x_α₂, :d₁]
@test getfield.(getfield.(getVar(bfm_gv), :val), :name) ==
      [:X₂, :Y₂, :Z₃, :α₃, :d₃, :X₁, :Y₁, :Z₁, :x_α₂,
       :d₂, :X₁, :Y₂, :Z₂, :α₁, :d₃, :α₃, :d₂, :x_α₂, :d₁]

@test getVarDict(e_gv1) == Dict(Symbolics.variable(:α, 1)=>2.0)

pairs_gv1 = sort(getVarDict(e_gv3, includeMapping=true) |> collect, by=x->string(x[1]))
strs_gv1 = getindex.(pairs_gv1, 1) .|> string
strs_gv2 = getindex.(pairs_gv1, 2) .|> string

@test strs_gv1[1][1:3] == "f_α" && strs_gv1[1][end-8:end] == "(x_α₂)"
@test strs_gv1[2:end] == ["x_α₂", "x_α₂^2"]
@test strs_gv2 == ["x_α₂^2", "1.05", "1.1025"]

@test getVarDict(gf_gv1) == Dict( Pair.(Symbolics.variable.([:α, :d], fill(1, 2)),
                                        [2.0, 1.0]) )

@test getVarDict(bf_gv1) == Dict( Pair.(Symbolics.variable.([:X, :Y, :Z, :α, :d],
                                                            fill(1, 5)),
                                        [1.0, 2.0, 3.0, 2.0, 1.0]) )

@test getVarDict(bf_gv1) == getVarDict([bf_gv1, gf_gv1])
@test getVarDict([bf_gv1, bf_gv2]) ==
      merge(getVarDict(bf_gv1), getVarDict(bf_gv2)) ==
      getVarDict([cen_gv1..., cen_gv2..., gf_gv1, gf_gv2]) ==
      getVarDict(vcat(bf_gv1.param|>collect, bf_gv2.param|>collect))
@test getVarDict(bfm_gv) ==
      getVarDict([bf_gv2, bf_gv4, bf_gv6]) ==
      getVarDict(pbs_gv0)

# function expressionOf, expressionOfCore
expr1 = expressionOf(gf1) |> string
@test expr1 == "exp(-2.0(r₁^2) - 2.0(r₂^2) - 2.0(r₃^2))" ||
      expr1 == "exp(-2.0(r₁^2) - (2.0(r₂^2)) - (2.0(r₃^2)))"

expr2 = expressionOf(bf1)[]|>string
@test expr2 == "exp(-2.0((r₁ - 1.0)^2) - 2.0((r₂ - 2.0)^2) - 2.0((r₃ - 3.0)^2))" ||
      expr2 == "exp(-2.0((r₁ - 1.0)^2) - (2.0((r₂ - 2.0)^2)) - (2.0((r₃ - 3.0)^2)))"
@test expressionOfCore(gf1) |> string == "d*exp(-α*(r₁^2 + r₂^2 + r₃^2))"
@test expressionOfCore(bf1)[]|>string == "d*exp(-α*((r₁ - X)^2 + (r₂ - Y)^2 + (r₃ - Z)^2))"

expStr = expressionOfCore(bf6)[] |> string
idx = findfirst('d', expStr)
@test isapprox(parse(Float64, expStr[1:idx-1]), 7.579425332952777, atol=1e-12)
@test expStr[idx:end] == "d*exp(-α*((r₁ - X)^2 + (r₂ - Y)^2 + (r₃ - Z)^2))*(α^0.75)" ||
      expStr[idx:end] == "d*(α^0.75)*exp(-α*((r₁ - X)^2 + (r₂ - Y)^2 + (r₃ - Z)^2))"

@test expressionOfCore(bfm1)[]|>string == "d*exp(-α*((r₁ - X)^2 + (r₂ - Y)^2 + (r₃ - Z)^2))"
expr_bfm2 =  expressionOfCore(bfm2)[] |> string
@test expr_bfm2 == "d*exp(-α*((r₁ - X)^2 + (r₂ - Y)^2 + (r₃ - Z)^2))*(r₁ - X)" ||
      expr_bfm2 == "d*(r₁ - X)*exp(-α*((r₁ - X)^2 + (r₂ - Y)^2 + (r₃ - Z)^2))"

# function inSymbols
sym1 = :a1
sym2 = :d
syms = [:a, :b, :c]
@test inSymbols(sym1, syms) == :a
@test !inSymbols(sym2, syms)

sym3 = Symbolics.variable(:a)
@test inSymbols(sym3, syms) == :a
@test inSymbols(sym3.val, syms) == :a
@test inSymbols(log(sym3).val, syms) == false

@test inSymbols(abs, syms) == false

# function varVal
vars = @variables X, Y, Z, F(X)
F = vars[end]
vars = vars[1:end-1]
f1(X,Y,Z) = log(X+0.5Y*Z)
expr = f1(X,Y,Z)
vals = [1, 3.334, -0.2]
d1 = Dict(vars .=> vals)
errT = 1e-10
@test varVal.(vars, Ref(d1)) == vals
@test isapprox(varVal(sum(vars), d1), sum(vals), atol=errT)
@test isapprox(varVal(prod(vars), d1), prod(vals), atol=errT)
@test isapprox(varVal(vars[1]^vars[2], d1), vals[1]^vals[2], atol=errT)
@test isapprox(varVal(expr, d1), f1(vals...), atol=errT)

f1 = Symbolics.variable(abs, T=Symbolics.FnType{Tuple{Any}, Real})(Z)
f2 = Symbolics.variable(:abs, T=Symbolics.FnType{Tuple{Any}, Real})(Z)
f1sym = f1.val
f2sym = f2.val
@test varVal(f1sym, d1) == abs(d1[Z])
@test varVal(f2sym, d1) == abs(d1[Z])

sp = 1.5
gb1 = GridBox(1,sp)
d2 = getVarDict(gb1.box |> flatten)
l = Symbolics.variable(:L,0)
vars2 = [i.val for i in keys(d2)]
diffs2 = Differential(l).(vars2)
exprs2 = map(x -> (x isa SymbolicUtils.Term) ? d2[x] : x, vars2)
diffExprs2 = Symbolics.derivative.(exprs2, Ref(l))
excs2 = Symbolics.build_function.(exprs2, l) .|> eval
vals2 = [f(sp) for f in excs2]
excdiffs2 = Symbolics.build_function.(diffExprs2, l) .|> eval
diffvals2 = [f(sp) for f in excdiffs2]

@test varVal.(vars2, Ref(d2)) == vals2
@test varVal.(diffs2, Ref(d2)) == diffvals2

end


text:
using Test, Dates, Julianda

foo = Julianda.Config.loadConfig("../config")
inst = "GBP_USD" # Instrument to test with

bar = Julianda.Instrument.getCandles(foo, inst, 1)
@test length(bar.candles) == 1
dt = DateTime(2019,7,1)
bar = Julianda.Instrument.getCandles(foo, inst, dt, dt)
@test length(bar.candles) == 500
bar = Julianda.Instrument.getCandles(foo, inst, dt, 1)
@test length(bar.candles) == 1
bar = Julianda.Instrument.getCandles(foo, inst, 1, dt)
@test length(bar.candles) == 1

bar = Julianda.Instrument.getOrderBook(foo, inst, dt)
@test length(bar.buckets) == 1745

bar = Julianda.Instrument.getPositionBook(foo, inst, dt)
@test length(bar.buckets) == 598

ext: jl

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
import YAML
using Twitter
using Printf

SECRET_KEYS_FILENAME = "secret_keys.yaml"

if !ispath(SECRET_KEYS_FILENAME)
    println(@sprintf "%s file not found! Please make sure to create the appropriate file containing your TWitter Developer keys." SECRET_KEYS_FILENAME)

else
    secrets = YAML.load_file(SECRET_KEYS_FILENAME)

    twitterauth(secrets["consumer_key"],
        secrets["consumer_token"],
        secrets["oauth_token"],
        secrets["oauth_secret"])

    user_timeline_of_tweets = get_user_timeline(screen_name = secrets["twitter_handle"], count=500)
    number_tweets_for_deletion = length(user_timeline_of_tweets)

    println(@sprintf "%d tweets found from account %s... Would you like to delete these tweets? (Y/N)" number_tweets_for_deletion secrets["twitter_handle"])
    response = readline()

    if (response == "Y")
        println(@sprintf """'%s' response detected. Tweets will now be deleted...""" response)
        for tweet in user_timeline_of_tweets
            tweet_id = tweet.id
            post_destroy_single_tweet_id(string(tweet_id))
        end
        println("The tweets have now been deleted. Exiting...")
    else
        println(@sprintf """ '%s' response detected. No tweets will be deleted. Exiting...""" response)
    end
end


text:
mutable struct canvas
    width::Int
    height::Int
    grid::Array{RayTracer.rt_color,2}
end

Base.convert(::Type{RayTracer.rt_color}, x::RayTracer.rt_color) = RayTracer.color(
    float(x.red), float(x.green), float(x.blue)
)

#Base.promote_rule(::Type{RayTracer.rt_color}, ::Type{RayTracer.rt_color{Integer}}) = RayTracer.rt_color{Float64}

function canvas(width::Int,height::Int)::canvas
    return canvas(width,height,fill(color(0.,0.,0.),(width,height)))
end

function write_pixel!(canv::canvas,x::Int,y::Int,c::RayTracer.rt_color)
    if x<=0 || x>canv.width
        println("Warning: not painting in x:$x, y:$y.")
        return nothing
    elseif y<=0 || y>canv.height
        println("Warning: not painting in x:$x, y:$y.")
        return nothing
    else
        canv.grid[x,y] = c
    end
end

function pixel_at(c::canvas, x::Int, y::Int)::RayTracer.rt_color
    return c.grid[x,y]
end

function canvas_to_ppm(c::canvas,filename::String="canvas.txt")::String
    payload = ""
    payload *= "P3\n"
    payload *= "$(c.width) $(c.height)\n"
    payload *= "255\n"
    for height in 1:c.height
        new_line = ""
        for width in 1:c.width
            px = pixel_at(c,width,height)
            for color in [px.red,px.green,px.blue]
                cl = min(max(convert(Int, round(color*255, digits=0)),0),255)
                if length(new_line*" "*string(cl))<71
                    new_line *= string(cl)*" "
                else
                    payload *= new_line[1:end-1]*"\n"
                    new_line = string(cl)*" "
                end
            end
        end
        payload *= new_line[1:end-1]*"\n"
    end
    open(filename, "w") do io
        write(io, payload)
    end
    return filename
end


text:
cd(@__DIR__)
import Pkg
Pkg.pkg"activate ."

ROOT_PATH = pwd()
push!(LOAD_PATH, ROOT_PATH, "src")

include(joinpath("src", "HelloWorldServer.jl"))
include("server.jl")


text:
# module Lupin

#=
function Base.getproperty(::Type{T}, sym::Symbol) where {T <: Enum}
    if sym in (:name, :parameters)
        getfield(T, sym)
    else
        for instance in instances(T)
            Symbol(instance) === sym && return instance
        end
        throw(ErrorException(string("type ", T, " has no field ", sym)))
    end
end
=#

# module Lupin


text:
@testset "Group algebras" begin
  G = small_group(8, 4)
  A = AlgGrp(FlintQQ, G)

  @testset "Regular matrix algebra" begin
    B, BtoA = Hecke.regular_matrix_algebra(A)

    @test dim(B) == dim(A)

    for i = 1:dim(A)
      for j = 1:dim(A)
        @test BtoA(B[i])*BtoA(B[j]) == BtoA(B[i]*B[j])
        @test (BtoA\A[i])*(BtoA\A[j]) == BtoA\(A[i]*A[j])
      end
    end
  end

  @testset "Generators" begin
    g, full_basis, v = gens(A, Val{true})

    @test length(full_basis) == dim(A)

    M = zero_matrix(FlintQQ, dim(A), dim(A))
    for i = 1:dim(A)
      Hecke.elem_to_mat_row!(M, i, full_basis[i])
    end
    @test rank(M) == dim(A)

    for i = 1:dim(A)
      b = full_basis[i]
      a = one(A)
      for (j, k) in v[i]
        a *= g[j]^k
      end
      @test b == a
    end
  end

end


text:
using HMatrices
using LoopVectorization
using StaticArrays
using LinearAlgebra

const EPS = 1e-8

# some simple implementation of laplace and helmholtz single layer matrices.
function laplace_matrix(X, Y)
    f = (x, y) -> begin
        EPS = 1e-8 # fudge factor to avoid division by zero
        d = norm(x - y) + EPS
        inv(4π * d)
    end
    KernelMatrix(f, X, Y)
    end

function helmholtz_matrix(X, Y, k)
    f = (x, y) -> begin
        EPS = 1e-8 # fudge factor to avoid division by zero
        d = norm(x - y) + EPS
        exp(im * k * d) * inv(4π * d)
    end
    KernelMatrix(f, X, Y)
end

function elastostatic_matrix(X, Y, μ, λ)
    f = (x, y) -> begin
        ν = λ / (2 * (μ + λ))
        r = x - y
        d = norm(r) + EPS
        RRT = r * transpose(r) # r ⊗ rᵗ
        return 1 / (16π * μ * (1 - ν) * d) * ((3 - 4 * ν) * I + RRT / d^2)
    end
    KernelMatrix(f, X, Y)
end

function elastosdynamic_matrix(X, Y, μ, λ, ω, ρ)
    f = (x, y) -> begin
        c1 = sqrt((λ + 2μ) / ρ)
        c2 = sqrt(μ / ρ)
        r = x - y
        d   = norm(r) + EPS
        RRT = r * transpose(r) # r ⊗ rᵗ
        s   = -im * ω
        z1  = s * d / c1
        z2  = s * d / c2
        α   = 4
        ψ   = exp(-z2) / d + (1 + z2) / (z2^2) * exp(-z2) / d - c2^2 / c1^2 * (1 + z1) / (z1^2) * exp(-z1) / d
        chi = 3 * ψ - 2 * exp(-z2) / d - c2^2 / c1^2 * exp(-z1) / d
        return 1 / (α * π * μ) * (ψ * I - chi * RRT / d^2) + x * transpose(y)
    end
    KernelMatrix(f, X, Y)
end

"""
    struct LaplaceMatrixVec{T,Td} <: AbstractKernelMatrix{T}

Vectorized version of [`laplace_matrix`](@ref).
"""
struct LaplaceMatrixVec{T,Td} <: AbstractKernelMatrix{T}
    X::Matrix{Td}
    Y::Matrix{Td}
    function LaplaceMatrixVec{T}(X::Matrix{Td}, Y::Matrix{Td}) where {T,Td}
        @assert size(X, 2) == size(Y, 2) == 3
        new{T,Td}(X, Y)
    end
end

Base.size(K::LaplaceMatrixVec) = size(K.X, 1), size(K.Y, 1)

# constructor based on Vector of StaticVector
function LaplaceMatrixVec{T}(_X::Vector{SVector{3,Td}}, _Y::Vector{SVector{3,Td}}) where {T,Td}
    X = reshape(reinterpret(Td, _X), 3, :) |> transpose |> collect
    Y = reshape(reinterpret(Td, _Y), 3, :) |> transpose |> collect
    LaplaceMatrixVec{T}(X, Y)
end
LaplaceMatrixVec(args...) = LaplaceMatrixVec{Float64}(args...) # default to Float64

function Base.getindex(K::LaplaceMatrixVec{T}, i::Int, j::Int)::T where {T}
    d2 = (K.X[i,1] - K.Y[j,1])^2 + (K.X[i,2] - K.Y[j,2])^2 + (K.X[i,3] - K.Y[j,3])^2
    d = sqrt(d2) + EPS
    return inv(4π * d)
end
function Base.getindex(K::LaplaceMatrixVec, I::UnitRange, J::UnitRange)
    T = eltype(K)
    m = length(I)
    n = length(J)
    Xv = view(K.X, I, :)
    Yv = view(K.Y, J, :)
    out = Matrix{T}(undef, m, n)
    @avx for j in 1:n
        for i in 1:m
            d2 = (Xv[i,1] - Yv[j,1])^2
            d2 += (Xv[i,2] - Yv[j,2])^2
            d2 += (Xv[i,3] - Yv[j,3])^2
            d = sqrt(d2) + EPS
            out[i,j] = inv(4 * π * d)
        end
    end
    return out
end
function Base.getindex(K::LaplaceMatrixVec, I::UnitRange, j::Int)
    return vec(K[I,j:j])
end
function Base.getindex(adjK::Adjoint{<:Any,LaplaceMatrixVec}, I::UnitRange, j::Int)
    K = parent(adjK)
    vec(K[j:j,I])
end

"""
    struct HelmholtzMatrix{T,Td,Tk} <: AbstractKernelMatrix{T}

Vectorized version of [`helmholtz_matrix`](@ref).
"""
struct HelmholtzMatrixVec{T,Td,Tk} <: AbstractKernelMatrix{T}
    X::Matrix{Td}
    Y::Matrix{Td}
    k::Tk
    function HelmholtzMatrixVec{T}(X::Matrix{Td}, Y::Matrix{Td}, k::Tk) where {T,Td,Tk}
        @assert size(X, 2) == size(Y, 2) == 3
        new{T,Td,Tk}(X, Y, k)
    end
end
HelmholtzMatrixVec(args...) = HelmholtzMatrixVec{ComplexF64}(args...)

Base.size(K::HelmholtzMatrixVec) = size(K.X, 1), size(K.Y, 1)

function HelmholtzMatrixVec{T}(_X::Vector{SVector{3,Td}}, _Y::Vector{SVector{3,Td}}, k) where {T,Td}
    X = reshape(reinterpret(Td, _X), 3, :) |> transpose |> collect
    Y = reshape(reinterpret(Td, _Y), 3, :) |> transpose |> collect
    HelmholtzMatrixVec{T}(X, Y, k)
end

function Base.getindex(K::HelmholtzMatrixVec{T}, i::Int, j::Int)::T where {T}
    d2 = (K.X[i,1] - K.Y[j,1])^2 + (K.X[i,2] - K.Y[j,2])^2 + (K.X[i,3] - K.Y[j,3])^2
    d  = sqrt(d2) + EPS
    return inv(4π * d) * exp(im * K.k * d)
end
function Base.getindex(K::HelmholtzMatrixVec{Complex{T}}, I::UnitRange, J::UnitRange) where {T}
    k = K.k
    m = length(I)
    n = length(J)
    Xv = view(K.X, I, :)
    Yv = view(K.Y, J, :)
    # since LoopVectorization does not (yet) support Complex{T} types, we will
    # reinterpret the output as a Matrix{T}, then use views. This can probably be
    # done better.
    out    = Matrix{Complex{T}}(undef, m, n)
    out_T  = reinterpret(T, out)
    out_r = @views out_T[1:2:end,:]
    out_i = @views out_T[2:2:end,:]
    @avx for j in 1:n
        for i in 1:m
            d2 = (Xv[i,1] - Yv[j,1])^2
            d2 += (Xv[i,2] - Yv[j,2])^2
            d2 += (Xv[i,3] - Yv[j,3])^2
            d  = sqrt(d2) + EPS
            s, c = sincos(k * d)
            zr = inv(4π * d) * c
            zi = inv(4π * d) * s
            out_r[i,j] = zr
            out_i[i,j] = zi
        end
    end
    return out
end
function Base.getindex(K::HelmholtzMatrixVec, I::UnitRange, j::Int)
    vec(K[I,j:j])
end
function Base.getindex(adjK::Adjoint{<:Any,<:HelmholtzMatrixVec}, I::UnitRange, j::Int)
    K = parent(adjK)
    conj!(K[j:j,I]) |> vec
end

function points_on_sphere(npts, R=1)
    theta = π * rand(npts)
    phi   = 2 * π * rand(npts)
    x     = @. sin(theta) * cos(phi)
    y     = @. R * sin(theta) * sin(phi)
    z     = @. R * cos(theta)
    data  = vcat(x', y', z')
    pts = reinterpret(SVector{3,Float64}, vec(data)) |> collect
    return pts
end

function points_on_cylinder(radius,n,shift=SVector(0,0,0))
    step     = 1.75*π*radius/sqrt(n)
    result          = Vector{SVector{3,Float64}}(undef,n)
    length          = 2*π*radius
    pointsPerCircle = length/step
    angleStep       = 2*π/pointsPerCircle
    for i=0:n-1
        x = radius * cos(angleStep*i)
        y = radius * sin(angleStep*i)
        z = step*i/pointsPerCircle
        result[i+1] = shift + SVector(x,y,z)
    end
    return result
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
module GalacticOptim

using Reexport
@reexport using DiffEqBase
@reexport using SciMLBase
using Requires
using DiffResults, ForwardDiff, Zygote, ReverseDiff, Tracker, FiniteDiff
@reexport using Optim, Flux
using Logging, ProgressLogging, Printf, ConsoleProgressMonitor, TerminalLoggers, LoggingExtras
using ArrayInterface, Base.Iterators

using ForwardDiff: DEFAULT_CHUNK_THRESHOLD
import DiffEqBase: OptimizationProblem, OptimizationFunction, AbstractADType
import SciMLBase: AbstractNoTimeSolution, build_solution, AbstractNonlinearProblem

import ModelingToolkit
import ModelingToolkit: AutoModelingToolkit
export AutoModelingToolkit

include("solve.jl")
include("function.jl")

export solve, EnsembleOptimizationProblem

export BBO, CMAEvolutionStrategyOpt

end # module


text:
import ..UncertainValues.TheoreticalDistributionScalarValue

import ..UncertainValues.AbstractUncertainOneParameterScalarValue
import ..UncertainValues.AbstractUncertainTwoParameterScalarValue
import ..UncertainValues.AbstractUncertainThreeParameterScalarValue

import ..UncertainValues.ConstrainedUncertainScalarValueOneParameter
import ..UncertainValues.ConstrainedUncertainScalarValueTwoParameter
import ..UncertainValues.ConstrainedUncertainScalarValueThreeParameter
import ..UncertainValues.TruncatedUncertainScalarKDE
import ..UncertainValues.AbstractUncertainScalarKDE

import ..UncertainValues.TheoreticalFittedUncertainScalar
import ..UncertainValues.ConstrainedUncertainScalarTheoreticalFit
import ..UncertainValues.FittedDistribution

import ..UncertainValues.getquantileindex

import Distributions.ValueSupport
import Distributions.Univariate
import Distributions.Distribution
import Distributions.support
import Distributions.Truncated
import Distributions: mean, std

import KernelDensity.UnivariateKDE

import StatsBase: quantile, Weights, mean, std

import Base.truncate

include("stats_on_truncated_distributions.jl")

include("truncate_kde.jl")
include("truncate_theoreticalscalar.jl")
include("truncate_theoretical_fitted.jl")
include("truncate_CertainValue.jl")
include("truncate_population.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ 8c51843c-9407-11ec-037a-0547f170ed7c
begin
	using PlutoUI
md"""
# Homework 6

Eric Nguyen
"""
end

# ╔═╡ 84fb9a5e-2ff5-44cf-b6d4-d483dd0ae5c1
md"""
### Problem 1

**(a).**
The solutions for ``(\alpha - 1)^2 > 4 \beta^2`` are ``y_1(x) = c_1 x^{r_1}`` and ``y_2(x) = c_2 x^{r_2}`` where

$r_1 = \frac{-(\alpha - 1) + \sqrt{(\alpha - 1)^2 - 4 \beta^2}}{2} \quad\text{ and }\quad r_2 = \frac{-(\alpha - 1) - \sqrt{(\alpha - 1)^2 - 4 \beta^2}}{2}$

The Wronskian of the two solutions gives us

$\begin{align*}
W[y_1(x), y_2(x)] &= \begin{vmatrix} y_1 & y_2 \\ y_1' & y_2' \end{vmatrix} \\
&= y_1(x) y_2'(x) - y_2(x) y_1'(x) \\
&= c_1 x^{r_1} ⋅ \frac{d}{dx} (c_2 x^{r_2}) - c_2 x^{r_2} ⋅ \frac{d}{dx} (c_1 x^{r_1}) \\
&= c_1 x^{r_1} ⋅ (r_2 c_2 x^{r_2 - 1}) - c_2 x^{r_2} ⋅ (c_1 r_1 x^{r_1 - 1}) \\

&= c_1 c_2 (r_2 x^{r_1} x^{r_2 - 1}) - c_1 c_2 (r_1 x^{r_2} x^{r_1 - 1}) \\
&= c_1 c_2 x^{r_1 + r_2 - 1} r_2 - c_1 c_2 x^{r_1 + r_2 - 1} r_1 \\
&= c_1 c_2 x^{r_1 + r_2 - 1} (r_2 - r_1)
\end{align*}$

Write ``r_1 + r_2`` and ``r_2 - r_1`` in terms of ``\alpha`` and ``\beta``.

$\begin{align*}
r_1 + r_2 &= \frac{\left[-(\alpha - 1) + \sqrt{(\alpha - 1)^2 - 4 \beta^2}\right] + \left[-(\alpha - 1) - \sqrt{(\alpha - 1)^2 - 4 \beta^2}\right]}{2} \\
&= \frac{-2(\alpha - 1)}{2} \\
&= 1 - \alpha
\end{align*}$

$\begin{align*}
r_2 - r_1 &= \frac{\left[-(\alpha - 1) - \sqrt{(\alpha - 1)^2 - 4 \beta^2}\right] - \left[-(\alpha - 1) + \sqrt{(\alpha - 1)^2 - 4 \beta^2}\right]}{2} \\
&= \frac{-2\sqrt{(\alpha - 1)^2 - 4 \beta^2}}{2} \\
&= -\sqrt{(\alpha - 1)^2 - 4 \beta^2}
\end{align*}$

$(\alpha - 1)^2 > 4 \beta^2 \implies r_2 - r_1 = -\sqrt{(\alpha - 1)^2 - 4 \beta^2} ≠ 0$

Substituting into the expression for the Wronskian, we get

$\begin{align*}
W[y_1(x), y_(x)] &= c_1 c_2 x^{(1 - \alpha) - 1} \left[-\sqrt{(\alpha - 1)^2 - 4 \beta^2}\right] \\
&= -c_1 c_2 x^{-\alpha} \sqrt{(\alpha - 1)^2 - 4 \beta^2} \\
&≠ 0 \implies \text{Linear independence}
\end{align*}$
"""

# ╔═╡ ea5236fd-9c69-4ffc-a73c-9f2b336f5444
md"""
**(b).**
For the case ``(\alpha - 1)^2 = 4 \beta^2`` then the quadratic equation only gives one solution.
Try

$y(x) = v(x) x^{r}$

$y' = v' x^r + rx^{r - 1} v$

$\begin{align*}
y'' &= v'' x^r + r x^{r - 1} v' + rx^{r - 1} v' + r(r - 1) x^{r - 2} v \\
&= x^r v'' + 2r x^{r - 1} v' + r(r - 1) x^{r - 2} v \\
&= v'' + 2r x^{-1} v' + r(r - 1) x^{-2} v \\
\end{align*}$

Find the homogeneous solution of ``v`` using the ansatz ``v = x^{\alpha}`` we get

$\begin{align*}
v'' + 2r x^{-1} v' + r(r - 1) x^{-2} v &= 0 \\
\left[\alpha(\alpha - 1) x^{\alpha - 2}\right] + 2r x^{-1} \left[\alpha x^{\alpha - 1}\right] + r(r - 1) x^{-2} [x^{\alpha}] &= 0 \\
\alpha^2 x^{\alpha - 2} - \alpha x^{\alpha - 2} + 2r \alpha x^{\alpha - 2} + r^2 x^{\alpha - 2} - r x^{\alpha - 2} &= 0 \\
\alpha^2 - \alpha + 2r \alpha + r^2 - r &= 0 \\
\alpha^2 + (2r - 1) \alpha + r^2 - r &= 0
\end{align*}$

$\begin{align*}
\alpha &= \frac{1 - 2r ± \sqrt{(2r - 1)^2 - 4(r^2 - r)}}{2} \\
&= \frac{1 - 2r ± \sqrt{4r^2 - 4r + 1 - 4r^2 + 4r}}{2} \\
&= \frac{1 - 2r ± 1}{2} \\
&= \frac{-2r}{2} = -r \quad\text{ or }\quad \frac{2 - 2r}{2} = 1 - r
\end{align*}$

$\alpha = -r \quad\text{ or }\quad \alpha = 1 - r$

So,

$v(x) = c_1 x^{-r} + c_2 x^{1 - r}$

which implies that the general solution to the differential equation is

$y(x) = (c_1 x^{-r} + c_2 x^{1 - r}) x^r = c_1 + c_2 x$

where

$y_1(x) = c_1 x^{-r} \qquad\text{and}\qquad y_2(x) = c_2 x^{1 - r}$

The Wronskian of the two solutions is

$\begin{align*}
W[y_1(x),y_2(x)] &= \begin{vmatrix} y_1 & y_2 \\ y_1' & y_2' \end{vmatrix} \\
&= y_1(x) y_2'(x) - y_2(x) y_1'(x) \\
&= c_1 x^{-r} (1 - r) c_2 x^{-r} - c_2 x^{1 - r} (-r) c_1 x^{-r - 1} \\
&= c_1 c_2 x^{-2r} (1 - r) + c_1 c_2 x^{-2r} r \\
&= c_1 c_2 x^{-2r} \left[(1 - r) + r\right] \\
&= c_1 c_2 x^{-2r} \\
&= c_1 c_2 x^{-2 \left[\frac{-(\alpha - 1)}{2}\right]} \\
&= c_1 c_2 x^{\alpha - 1} \\
&≠ 0 \implies \text{Linear independence for } x > 0
\end{align*}$
"""

# ╔═╡ e2708eda-5707-4c8f-a9b8-74929998d37c
md"""
### Problem 2

The Bessel function of the integer order ``m`` is defined as

$J_m(x) = \sum_{k = 0}^∞ \frac{(-1)^k}{k! (m + k)!} \left(\frac{x}{2}\right)^{m + 2k}$

We can then find the Bessel function of the integer order ``-m``, using the ``\Gamma``-Function to interpret ``n!`` for ``n < 0``.

$\begin{align*}
J_{-m}(x) &= \sum_{k = 0}^∞ \frac{(-1)^k}{k! (k - m)!} \left(\frac{x}{2}\right)^{2k - m} \\
&= \sum_{k = m}^∞ \frac{(-1)^{k}}{\Gamma(k + 1) \Gamma(-m + k + 1)} \left(\frac{x}{2}\right)^{2k - m} \\
&= \sum_{k = 0}^∞ \frac{(-1)^{k+m}}{\Gamma(k + m + 1) \Gamma(-m + k + m + 1)} \left(\frac{x}{2}\right)^{2(k + m) - m} \\
&= \left(\frac{x}{2}\right)^{-m} \sum_{k = 0}^∞ \frac{(-1)^{k+m}}{\Gamma(k + m + 1) \Gamma(k + 1)} \left(\frac{x}{2}\right)^{2k + 2m} \\
&= \left(\frac{x}{2}\right)^{-m} \left(\frac{x}{2}\right)^{m} \sum_{k = 0}^∞ \frac{(-1)^{k}(-1)^m}{\Gamma(k + m + 1) \Gamma(k + 1)} \left(\frac{x}{2}\right)^{2k+m} \\
&= (-1)^m \sum_{k = 0}^∞ \frac{(-1)^{k}(-1)^m}{\Gamma(k + m + 1) \Gamma(k + 1)} \left(\frac{x}{2}\right)^{2k+m} \\
&= (-1)^m \sum_{k = 0}^∞ \frac{(-1)^{k}}{(k + m)! k!} \left(\frac{x}{2}\right)^{2k+m} \\
&= (-1)^m J_m(x)
\end{align*}$

This means that

$y(x) = c_1 J_m(x) + c_2 J_{-m}(x)$

cannot be a general solution to Bessel's equation for ``m ∈ ℤ`` because the Gamma function is undefined for ``-m``.
"""

# ╔═╡ 6b011578-0d89-49dc-8130-886a6c9518b2
md"""
### Problem 3

$\begin{align*}
j_0(x) &= \frac{\sin{x}}{x} \\
j_0'(x) &= \frac{x \cos{x} - \sin{x}}{x^2} \\
j_0''(x) &= \frac{x^2 (-x \sin{x}) - (x \cos{x} - \sin{x}) ⋅ 2x}{x^4} \\
&= \frac{-x^3 \sin{x} - 2x^2 \cos{x} + 2x \sin{x}}{x^4} \\
&= \frac{-x^2 \sin{x} - 2x \cos{x} + 2 \sin{x}}{x^3} \\
\end{align*}$

$\begin{align*}
R(r) &= j_0(kr) = \frac{\sin{kr}}{kr} \\
R'(r) &= j_0'(kr) = \frac{kr \cos{kr} - \sin{kr}}{k^2 r^2} \\
R''(r) &= j_0''(kr) = \frac{-k^2 r^2 \sin{kr} - 2kr \cos{kr} + 2 \sin{kr}}{k^3 r^3} \\
\end{align*}$

$\begin{align*}
k^2 r^2 R(r) &= kr \sin{kr} \\
2r R'(r) &= \frac{2kr \cos{kr} - \sin{kr}}{k^2 r} \\
r^2 R''(r) &= \frac{-k^2 r^2 \sin{kr} - 2kr \cos{kr} + 2 \sin{kr}}{k^3 r} \\
\end{align*}$

It doesn't work unless I multiply the terms by ``k^2, k, 1`` respectively...

$\begin{align*}
&r^2 k^2 R''(r) + 2rk R'(r) + k^2 r^2 R(r) \\
&= \frac{-k^2 r^2 \sin{kr} - 2kr \cos{kr} + 2 \sin{kr}}{k r} + \frac{2kr \cos{kr} - 2 \sin{kr}}{k r} + kr \sin{kr} \\
&= \frac{-k^2 r^2 \sin{kr} - 2kr \cos{kr} + 2 \sin{kr}}{kr} + \frac{2kr \cos{kr} - 2 \sin{kr}}{kr} + \frac{k^2 r^2 \sin{kr}}{kr} \\
&= -k^2 r^2 \sin{kr} - 2kr \cos{kr} + 2 \sin{kr} + 2kr \cos{kr} - 2 \sin{kr} + k^2 r^2 \sin{kr} \\
&= 0
\end{align*}$

Indeed,

$j_0(x) = \frac{\sin{x}}{x}$

solves the Helmholtz Equation where ``ℓ = 0``.
"""

# ╔═╡ abcc322f-8dac-4529-bd50-c1917c47d29c
md"""
### Problem 4

The Rodrigues Formula is

$P_ℓ(x) = \frac{1}{2^ℓ ℓ!} \frac{d^ℓ}{dx^ℓ} (x^2 - 1)^ℓ = 1$

We can use this formula to derive the first three Legendre Polynomials ``P_0(x)``, ``P_1(x)``, and ``P_2(x)``.

$\begin{align*}
P_0(x) &= \frac{1}{2^0 0!} \frac{d^0}{dx^0} (x^2 - 1)^0 = 1 \\
P_1(x) &= \frac{1}{2^1 1!} \frac{d^1}{dx^1} (x^2 - 1)^1 = \frac{1}{2} (2x) = x \\
P_2(x) &= \frac{1}{2^2 2!} \frac{d^2}{dx^2} (x^2 - 1)^2 \\
&= \frac{1}{8} \frac{d}{dx}[2(x^2 - 1) ⋅ 2x] \\
&= \frac{1}{2} \frac{d}{dx}(x^3 - x) \\
&= \frac{1}{2} (3x^2 - 1) \\
\end{align*}$

Indeed, these results match the ones given in the course notes.
"""

# ╔═╡ 39f17932-dd54-49d0-baa1-786bb78a218d
md"""
### Problem 5

**(a).**

$\begin{align*}
F_1 = 3m x_1''(t) &= -kx_1 + k(x_2 - x_1) \\
0 &= 3mx_1''(t) + kx_1 - k(x_2 - x_1) \\
0 &= 3mx_1''(t) - kx_2 \\
0 &= 3x_1''(t) + 2 {\omega_0}^2 x_1 - {\omega_0}^2 x_2
\end{align*}$

$\begin{align*}
F_2 = 2m x_2''(t) &= -k(x_2 - x_1) \\
0 &= 2m x_2''(t) + k(x_2 - x_1) \\
0 &= 2m x_2''(t) + kx_2 - kx_1 \\
0 &= 2x_2''(t) + {\omega_0}^2 x_2 - {\omega_0}^2 x_1
\end{align*}$

**(b).**
Using the ansatz

$x_1(t) = a_1 e^{i \omega t} \qquad\text{and}\qquad x_2(t) = a_2 e^{i \omega t}$

we get

$\begin{align*}
-3 \omega^2 a_1 + 2 {\omega_0}^2 a_1 - {\omega_0}^2 a_2 &= 0 \\
-2 \omega^2 a_2 + {\omega_0}^2 a_2 - {\omega_0}^2 a_1 &= 0
\end{align*}$

Combining terms we get

$\begin{align*}
(-3 \omega^2 + 2 {\omega_0}^2) a_1 - {\omega_0}^2 a_2 &= 0 \\
-{\omega_0}^2 a_1 + (- 2 \omega^2 + {\omega_0}^2) a_2 &= 0
\end{align*}$

The setting the ratios of ``a_1`` and ``a_2`` equal to each other gives us

$\frac{-3 \omega^2 + 2 {\omega_0}^2}{-{\omega_0}^2} = \frac{-{\omega_0}^2}{-2\omega^2 + {\omega_0}^2}$

Solving for ``{\omega}^2`` we get

$\begin{align*}
(-3 \omega^2 + 2 {\omega_0}^2)(-2\omega^2 + {\omega_0}^2) &= {\omega_0}^4 \\
6 \omega^4 - 7 \omega^2 {\omega_0}^2 + 2 {\omega_0}^4 &= {\omega_0}^4 \\
6 \omega^4 - 7 \omega^2 {\omega_0}^2 + {\omega_0}^4 &= 0 \\
\end{align*}$

$\omega^2 = \frac{7 {\omega_0}^2 ± \sqrt{49 {\omega_0}^4 - 4(6) {\omega_0}^4}}{2(6)} = \frac{7 {\omega_0}^2 ± \sqrt{25 {\omega_0}^4}}{12} = \frac{7 {\omega_0}^2 ± 5 {\omega_0}^2}{12}$

$\omega^2 = \frac{1}{6} {\omega_0}^2 \qquad\text{or}\qquad \omega^2 = {\omega_0}^2$

Thus the eigenfrequencies are

$\begin{align*}
\omega_1 &= \frac{1}{\sqrt{6}} \omega_0 \\
\omega_2 &= -\frac{1}{\sqrt{6}} \omega_0 \\
\omega_3 &= \omega_0 \\
\omega_4 &= -\omega_0
\end{align*}$

The eigenmode ``x_{-}(t)`` oscillates at a frequency of ``\omega_0`` and the eigenmode ``x_{+}(t)`` oscillates at a frequency of ``\omega_0 / \sqrt{6}``.

**(c).**

$x_1(t) = ae^{i \omega_0 t / \sqrt{6}} + be^{-i \omega_0 t / \sqrt{6}} + ce^{i \omega_0 t} + de^{-i \omega_0 t}$

$x_2(t) = ae^{i \omega_0 t / \sqrt{6}} + be^{-i \omega_0 t / \sqrt{6}} - ce^{i \omega_0 t} - de^{-i \omega_0 t}$

**(d).**
"""

# ╔═╡ da7fe8c2-be38-4247-ab84-29ddfe0efd1c
PlutoUI.LocalResource("img/hw_6_p5_d.png", :width => 450)

# ╔═╡ 8b5f957a-d209-4dcb-ab27-4088b95e4ad4
md"""
**(e).**
"""

# ╔═╡ 6859c1ba-f29f-466c-9ee5-ecf99f033de1
PlutoUI.LocalResource("img/hw_6_p5_e.png", :width=>450)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.35"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.2"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "85bf3e4bd279e405f91489ce518dedb1e32119cb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.35"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─8c51843c-9407-11ec-037a-0547f170ed7c
# ╟─84fb9a5e-2ff5-44cf-b6d4-d483dd0ae5c1
# ╟─ea5236fd-9c69-4ffc-a73c-9f2b336f5444
# ╟─e2708eda-5707-4c8f-a9b8-74929998d37c
# ╟─6b011578-0d89-49dc-8130-886a6c9518b2
# ╟─abcc322f-8dac-4529-bd50-c1917c47d29c
# ╟─39f17932-dd54-49d0-baa1-786bb78a218d
# ╟─da7fe8c2-be38-4247-ab84-29ddfe0efd1c
# ╟─8b5f957a-d209-4dcb-ab27-4088b95e4ad4
# ╟─6859c1ba-f29f-466c-9ee5-ecf99f033de1
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
mutable struct VisitedList
    list::Vector{UInt8}
    visited_value::UInt8
end

"""
    VisitedList(num_elements)
List with `num_elements` entries to keep track of wether they have been visited or not.
Check status with
    isvisited(vl::VisitedList, idx) → Bool

and visit with
    visit!(vl::VisitedList, idx)

To reset the list, call 'reset!(vl)'.
"""
VisitedList(num_elements) = VisitedList(Vector{UInt8}(fill(zero(UInt8),num_elements)),1)

function reset!(vl::VisitedList)
    vl.visited_value += one(UInt8)
    if vl.visited_value == zero(UInt8)
        vl.list .= zero(UInt8)
        vl.visited_value += one(UInt8)
    end
end
isvisited(vl::VisitedList, idx::Integer) = vl.list[idx] == vl.visited_value
isvisited(vl::VisitedList, q::Neighbor) = isvisited(vl, q.idx)
visit!(vl::VisitedList, idx::Integer) = vl.list[idx] = vl.visited_value
visit!(vl::VisitedList, q::Neighbor) = visit!(vl, q.idx)

## Multithreaded pool-management of VisitedList's
struct VisitedListPool
    pool::Vector{VisitedList}
    num_elements::Int
end

"""
    VisitedListPool(num_lists, num_elements)
A thread-stable container for multiple `VisitedList`s initialized with
`num_lists` lists with each `num_elements` entries.

To retrieve a list, call `get_list(vlp::VisitedListPool)`,
and to release ist, call `release_list(vlp, vl::VisitedList)`.
"""
function VisitedListPool(num_lists::Real, num_elements::Real)
    pool = [VisitedList(num_elements) for n=1:num_lists]
    VisitedListPool(pool, num_elements)
end

function get_list(vlp::VisitedListPool)
    if length(vlp.pool) > 0
        vl = pop!(vlp.pool)
        reset!(vl)
    else
        vl = VisitedList(vlp.num_elements)
    end
    return vl
end

function release_list(vlp::VisitedListPool, vl::VisitedList)
    push!(vlp.pool,vl)
end


text:
Random.seed!(0)

#sig = [1 2; 2 2]
#sig = [1 1; 2 3]
#sig = [2 3]
#sig = [2 2]
sig = [3 2; 2 3]

S = random_S0Graph(sig)
T = complement(S)

n = S.n
D = S.D
J = block_expander(S.sig)

w = random_bounded(S.n)

@time opt0, x1, Z1 = dsw(S, w, eps=solver_eps)
if true
    @time opt1, _, x2, Z2 = dsw_via_complement(T, w, eps=solver_eps)
    @test opt1 ≈ opt0  atol=tol
else
    @time opt1, _, y = dsw_via_complement(T, w, eps=solver_eps)
    @test opt1 ≈ opt0  atol=tol
    @time opt2, x2, Z2 = dsw(T, y, eps=solver_eps)
    @test opt2 ≈ 1  atol=tol
end

x1 /= opt0
Z1 /= opt0
Z1 = J * Z1 * J'
Z2 = J * Z2 * J'
@test Z1 ≈ Z1'
@test Z2 ≈ Z2'
Z1 = Hermitian(Z1)
Z2 = Hermitian(Z2)

@test tr(√D * x1 * √D * x2) ≈ 1  atol=tol
@test partialtrace(Z1, 1, [n,n]) ≈ transpose(x1)  atol=tol
@test partialtrace(Z2, 1, [n,n]) ≈ transpose(x2)  atol=tol

v1 = reshape(conj(x1), n^2)
v2 = reshape(conj(x2), n^2)
Q1 = [1 v1'; v1 Z1]
Q2 = [1 v2'; v2 Z2]
@test Q1 ≈ Q1'
@test Q2 ≈ Q2'
Q1 = Hermitian(Q1)
Q2 = Hermitian(Q2)
D2 = cat([ 1, -kron(√D, √D) ]..., dims=(1,2))

@test minimum(eigvals(Q1)) > -tol
@test minimum(eigvals(Q2)) > -tol
@test tr(Q1 * D2 * Q2 * D2) ≈ 0  atol=tol
@test Z1 * kron(√D, √D) * v2 ≈ v1  atol=tol
@test Z2 * kron(√D, √D) * v1 ≈ v2  atol=tol
@test Z1 * kron(√D, √D) * Z2 ≈ v1 * v2'  atol=tol
@test Z1 * kron(eye(n), D) * Z2 ≈ v1 * v2'  atol=tol

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
module RobotDynamics

using Rotations
using StaticArrays
using LinearAlgebra
using ForwardDiff
using FiniteDiff
using RecipesBase
using SparseArrays

using Rotations: skew
using StaticArrays: SUnitRange

include("utils.jl")
include("knotpoint.jl")
include("functionbase.jl")
include("statevectortype.jl")
include("scalar_function.jl")
include("dynamics.jl")
include("discrete_dynamics.jl")

include("jacobian_gen.jl")
include("discretized_dynamics.jl")
include("integration.jl")

include("liestate.jl")
include("rbstate.jl")
include("rigidbody.jl")

include("jacobian.jl")
include("trajectories.jl")
include("plot_recipes.jl")

include("deprecate.jl")

export KnotPoint

const DataVector{T} = Union{Vector{T},StaticVector{<:Any,T},SubArray{T,1}}

using FiniteDiff: compute_epsilon
function FiniteDiff.finite_difference_gradient!(
    df::StridedVector{<:Number},
    f,
    x::StaticVector,
    cache::FiniteDiff.GradientCache{T1,T2,T3,T4,fdtype,returntype,inplace};
    relstep=FiniteDiff.default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    dir=true) where {T1,T2,T3,T4,fdtype,returntype,inplace}

    # c1 is x1 if we need a complex copy of x, otherwise Nothing
    # c2 is Nothing
    fx, c1, c2, c3 = cache.fx, cache.c1, cache.c2, cache.c3
    copyto!(c3,x)
    if fdtype == Val(:forward)
        for i ∈ eachindex(x)
            epsilon = compute_epsilon(fdtype, x[i], relstep, absstep, dir)
            x_old = x[i]
            fx0 = f(x)
            c3[i] += epsilon
            dfi = (f(c3) - fx0) / epsilon
            c3[i] = x_old

            df[i] = real(dfi)
        end
    else
        FiniteDiff.fdtype_error(returntype)
    end
    df
end

end # module


text:
# Variables

"""
    struct AddVariableNotAllowed <: NotAllowedError
        message::String # Human-friendly explanation why the attribute cannot be set
    end

An error indicating that variables cannot be added to the model.
"""
struct AddVariableNotAllowed <: NotAllowedError
    message::String # Human-friendly explanation why the attribute cannot be set
end
AddVariableNotAllowed() = AddVariableNotAllowed("")

operation_name(::AddVariableNotAllowed) = "Adding variables"

"""
    add_variables(model::ModelLike, n::Int)::Vector{VariableIndex}

Add `n` scalar variables to the model, returning a vector of variable indices.

A [`AddVariableNotAllowed`](@ref) error is thrown if adding variables cannot be
done in the current state of the model `model`.
"""
add_variables(model::ModelLike, n) = throw(AddVariableNotAllowed())

"""
    add_variable(model::ModelLike)::VariableIndex

Add a scalar variable to the model, returning a variable index.

A [`AddVariableNotAllowed`](@ref) error is thrown if adding variables cannot be
done in the current state of the model `model`.
"""
add_variable(model::ModelLike) = throw(AddVariableNotAllowed())


text:
using MathOptInterface
const MOI = MathOptInterface

# Test, whether the cut callback is actually called and whether
# `callback_value` works as intended.
@testset "obtaining the LP-solution" begin
    atol, rtol = 1e-6, 1e-6

    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    calls = 0
    x_val::Float64 = -1
    y_val::Float64 = -1
    function cutcallback(cb_data)
        x_val = MOI.get(optimizer,
                        MOI.CallbackVariablePrimal{SCIP.CutCbData}(cb_data),
                        x)
        y_val = MOI.get(optimizer,
                        MOI.CallbackVariablePrimal{SCIP.CutCbData}(cb_data),
                        y)
        calls += 1
    end
    MOI.set(optimizer, MOI.UserCutCallback(), cutcallback)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # The cut callback was called and obtaining the LP-solution worked.
    @test calls >= 1
    @test x_val + y_val >= 1.0 - min(atol, 1.0 * rtol)

    # SCIP found an optimal solution
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol

    # free the problem
    finalize(mscip)
end

# Test, whether adding cuts within cut callbacks via `submit` works [1/2].
@testset "cutting one optimal solution" begin
    atol, rtol = 1e-6, 1e-6

    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    calls = 0
    function cutcallback(cb_data)
        MOI.submit(optimizer, MOI.UserCut{SCIP.CutCbData}(cb_data),
                   MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0),
                   MOI.LessThan(0.0))
        calls += 1
    end
    MOI.set(optimizer, MOI.UserCutCallback(), cutcallback)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # The cut callback was called.
    @test calls >= 1

    # SCIP found the single remaining optimal solution
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol

    # free the problem
    finalize(mscip)
end

# Test, whether adding cuts within cut callbacks via `submit` works [2/2].
@testset "cutting another optimal solution" begin
    atol, rtol = 1e-6, 1e-6

    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    calls = 0
    function cutcallback(cb_data)
        MOI.submit(optimizer, MOI.UserCut{SCIP.CutCbData}(cb_data),
                   MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0),
                   MOI.LessThan(0.0))
        calls += 1
    end
    MOI.set(optimizer, MOI.UserCutCallback(), cutcallback)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # The cut callback was called.
    @test calls >= 1

    # SCIP found the single remaining optimal solution
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol

    # free the problem
    finalize(mscip)
end


text:
module CSetAutomorphisms

export common, canonical_hash, apply_automorphism, autos, CDict, to_vizstate,
       color_saturate, canonical_iso_nauty, init_graphs, test_iso, to_lg,
       graph_to_cset, graph_to_lg, Labeled, TheoryDecGraph

include("./Perms.jl")
include("./ColorRefine.jl")
include("./NautyInterface.jl")
include("./Canonical.jl")
include("./Visualization.jl")
include("./TestHelp.jl")

end # module


text:
using BenchmarkTools, RangeEnclosures
using RangeEnclosures: relative_precision

SUITE = BenchmarkGroup()

include("Daisy/daisy.jl")


text:
"""
    coef(MLM)

Extracts coefficients from Mlm object

# Arguments

- MLM = Mlm object

# Value

2d array of floats

"""
function coef(MLM::Mlm)

    return MLM.B
end

"""
    predict(MLM, newPredictors)

Calculates new predictions based on Mlm object

# Arguments

- MLM = Mlm object
- newPredictors = Predictors object. Defaults to the `data.predictors` field
  in the Mlm object used to fit the model.

# Value

Response object

"""
function predict(MLM::Mlm, newPredictors::Predictors=MLM.data.predictors)

  	# Include X and Z intercepts in new data if necessary
  	if MLM.data.predictors.hasXIntercept==true &&
       newPredictors.hasXIntercept==false
    	newPredictors.X = add_intercept(newPredictors.X)
    	newPredictors.hasXIntercept = true
    	println("Adding X intercept to newPredictors.")
  	end
  	if MLM.data.predictors.hasZIntercept==true &&
       newPredictors.hasZIntercept==false
    	newPredictors.Z = add_intercept(newPredictors.Z)
    	newPredictors.hasZIntercept = true
    	println("Adding Z intercept to newPredictors.")
  	end

  	# Remove X and Z intercepts in new data if necessary
  	if MLM.data.predictors.hasXIntercept==false &&
       newPredictors.hasXIntercept==true
    	newPredictors.X = remove_intercept(newPredictors.X)
    	newPredictors.hasXIntercept = false
    	println("Removing X intercept from newPredictors.")
  	end
  	if MLM.data.predictors.hasZIntercept==false &&
       newPredictors.hasZIntercept==true
    	newPredictors.Z = remove_intercept(newPredictors.Z)
    	newPredictors.hasZIntercept = false
    	println("Removing Z intercept from newPredictors.")
  	end

    # Calculate predictions
    return Response(calc_preds(newPredictors.X, newPredictors.Z, coef(MLM)))
end

"""
    fitted(MLM)

Calculate fitted values of an Mlm object

# Arguments

- MLM = Mlm object

# Value

Response object

"""
function fitted(MLM::Mlm)

    # Call the predict function with default newPredictors
    return predict(MLM)
end

"""
    resid(MLM, newData)

Calculates residuals of an Mlm object

# Arguments

- MLM = Mlm object
- newData = RawData object. Defaults to the `data` field in the Mlm object
  used to fit the model.

# Value

2d array of floats

"""
function resid(MLM::Mlm, newData::RawData=MLM.data)

    # Include X and Z intercepts in new data if necessary
    if MLM.data.predictors.hasXIntercept==true &&
       newData.predictors.hasXIntercept==false
        newData.predictors.X = add_intercept(newData.predictors.X)
        newData.predictors.hasXIntercept = true
        newData.p = newData.p + 1
        println("Adding X intercept to newData.")
    end
    if MLM.data.predictors.hasZIntercept==true &&
       newData.predictors.hasZIntercept==false
        newData.predictors.Z = add_intercept(newData.predictors.Z)
        newData.predictors.hasZIntercept = true
        newData.q = newData.q + 1
        println("Adding Z intercept to newData.")
    end

    # Remove X and Z intercepts in new data if necessary
    if MLM.data.predictors.hasXIntercept==false &&
       newData.predictors.hasXIntercept==true
        newData.predictors.X = remove_intercept(newData.predictors.X)
        newData.predictors.hasXIntercept = false
        newData.p = newData.p - 1
        println("Removing X intercept from newData.")
    end
    if MLM.data.predictors.hasZIntercept==false &&
       newData.predictors.hasZIntercept==true
        newData.predictors.Z = remove_intercept(newData.predictors.Z)
        newData.predictors.hasZIntercept = false
        newData.q = newData.q - 1
        println("Removing Z intercept from newData.")
    end

	# Calculate residuals
    return calc_resid(get_X(newData), get_Y(newData), get_Z(newData),
                      coef(MLM))
end


text:
"""
`framedict[method]` returns the `FrameCode` for `method`. For `@generated` methods,
see [`genframedict`](@ref).
"""
const framedict = Dict{Method,FrameCode}()                # essentially a method table for lowered code

"""
`genframedict[(method,argtypes)]` returns the `FrameCode` for a `@generated` method `method`,
for the particular argument types `argtypes`.

The framecodes stored in `genframedict` are for the code returned by the generator
(i.e, what will run when you call the method on particular argument types);
for the generator itself, its framecode would be stored in [`framedict`](@ref).
"""
const genframedict = Dict{Tuple{Method,Type},FrameCode}() # the same for @generated functions

"""
`meth ∈ compiled_methods` indicates that `meth` should be run using [`Compiled`](@ref)
rather than recursed into via the interpreter.
"""
const compiled_methods = Set{Method}()

"""
`mod ∈ compiled_modules` indicates that any method in `mod` should be run using [`Compiled`](@ref)
rather than recursed into via the interpreter.
"""
const compiled_modules = Set{Module}()

const junk_framedata = FrameData[] # to allow re-use of allocated memory (this is otherwise a bottleneck)
const junk_frames = Frame[]
debug_recycle() = false
@noinline function _check_frame_not_in_junk(frame)
    @assert frame.framedata ∉ junk_framedata
    @assert frame ∉ junk_frames
end

@inline function recycle(frame)
    debug_recycle() && _check_frame_not_in_junk(frame)
    push!(junk_framedata, frame.framedata)
    push!(junk_frames, frame)
end

function return_from(frame::Frame)
    recycle(frame)
    frame = caller(frame)
    frame === nothing || (frame.callee = nothing)
    return frame
end

function clear_caches()
    empty!(junk_framedata)
    empty!(framedict)
    empty!(genframedict)
    empty!(junk_frames)
    for bp in breakpoints()
        empty!(bp.instances)
    end
end

const empty_svec = Core.svec()

function namedtuple(kwargs)
    names, types, vals = Symbol[], [], []
    for pr in kwargs
        if isa(pr, Expr)
            push!(names, pr.args[1])
            val = pr.args[2]
            push!(types, typeof(val))
            push!(vals, val)
        elseif isa(pr, Pair)
            push!(names, pr.first)
            val = pr.second
            push!(types, typeof(val))
            push!(vals, val)
        else
            error("unhandled entry type ", typeof(pr))
        end
    end
    return NamedTuple{(names...,), Tuple{types...}}(vals)
end

get_source(meth::Method) = Base.uncompressed_ast(meth)

function get_source(g::GeneratedFunctionStub, env)
    b = g(env..., g.argnames...)
    b isa CodeInfo && return b
    return eval(b)
end

function copy_codeinfo(code::CodeInfo)
    @static if VERSION < v"1.1.0-DEV.762"
        newcode = ccall(:jl_new_struct_uninit, Any, (Any,), CodeInfo)::CodeInfo
        for (i, name) in enumerate(fieldnames(CodeInfo))
            if isdefined(code, name)
                val = getfield(code, name)
                ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), newcode, i-1, val===nothing || isa(val, Union{Type, Method}) ? val : copy(val))
            end
        end
        return newcode
    else
        # Inline this when support for VERSION above is dropped
        return copy(code)
    end
end

"""
    frun, allargs = prepare_args(fcall, fargs, kwargs)

Prepare the complete argument sequence for a call to `fcall`. `fargs = [fcall, args...]` is a list
containing both `fcall` (the `#self#` slot in lowered code) and the positional
arguments supplied to `fcall`. `kwargs` is a list of keyword arguments, supplied either as
list of expressions `:(kwname=kwval)` or pairs `:kwname=>kwval`.

For non-keyword methods, `frun === fcall`, but for methods with keywords `frun` will be the
keyword-sorter function for `fcall`.

# Example

```jldoctest
julia> mymethod(x) = 1
mymethod (generic function with 1 method)

julia> mymethod(x, y; verbose=false) = nothing
mymethod (generic function with 2 methods)

julia> JuliaInterpreter.prepare_args(mymethod, [mymethod, 15], ())
(mymethod, Any[mymethod, 15])

julia> JuliaInterpreter.prepare_args(mymethod, [mymethod, 1, 2], [:verbose=>true])
(var"#mymethod##kw"(), Any[var"#mymethod##kw"(), (verbose = true,), mymethod, 1, 2])
```
"""
function prepare_args(@nospecialize(f), allargs, kwargs)
    if !isempty(kwargs)
        f = Core.kwfunc(f)
        allargs = Any[f, namedtuple(kwargs), allargs...]
    elseif f === Core._apply
        f = to_function(allargs[2])
        allargs = append_any((allargs[2],), allargs[3:end]...)
    end
    return f, allargs
end

if VERSION < v"1.2-" || !isdefined(Core.Compiler, :specialize_method)
    specialize_method(method::Method, @nospecialize(atypes), sparams::SimpleVector) =
        Core.Compiler.code_for_method(method, atypes, sparams, typemax(UInt))
else
    const specialize_method = Core.Compiler.specialize_method
end

function prepare_framecode(method::Method, @nospecialize(argtypes); enter_generated=false)
    sig = method.sig
    if method.module ∈ compiled_modules || method ∈ compiled_methods
        return Compiled()
    end
    # Get static parameters
    (ti, lenv::SimpleVector) = ccall(:jl_type_intersection_with_env, Any, (Any, Any),
                        argtypes, sig)::SimpleVector
    enter_generated &= is_generated(method)
    if is_generated(method) && !enter_generated
        framecode = get(genframedict, (method, argtypes::Type), nothing)
    else
        framecode = get(framedict, method, nothing)
    end
    if framecode === nothing
        if is_generated(method) && !enter_generated
            # If we're stepping into a staged function, we need to use
            # the specialization, rather than stepping through the
            # unspecialized method.
            code = Core.Compiler.get_staged(specialize_method(method, argtypes, lenv))
            code === nothing && return nothing
            generator = false
        else
            if is_generated(method)
                code = get_source(method.generator, lenv)
                generator = true
            else
                code = get_source(method)
                generator = false
            end
        end
        code = code::CodeInfo
        # Currenly, our strategy to deal with llvmcall can't handle parametric functions
        # (the "mini interpreter" runs in module scope, not method scope)
        if (!isempty(lenv) && (hasarg(isidentical(:llvmcall), code.code) ||
                              hasarg(a->is_global_ref(a, Base, :llvmcall), code.code))) ||
                hasarg(isidentical(:iolock_begin), code.code)
            return Compiled()
        end
        framecode = FrameCode(method, code; generator=generator)
        if is_generated(method) && !enter_generated
            genframedict[(method, argtypes)] = framecode
        else
            framedict[method] = framecode
        end
    end
    return framecode, lenv
end

function get_framecode(method)
    framecode = get(framedict, method, nothing)
    if framecode === nothing
        code = get_source(method)
        framecode = FrameCode(method, code; generator=false)
        framedict[method] = framecode
    end
    return framecode
end

"""
    framecode, frameargs, lenv, argtypes = prepare_call(f, allargs; enter_generated=false)

Prepare all the information needed to execute lowered code for `f` given arguments `allargs`.
`f` and `allargs` are the outputs of [`prepare_args`](@ref).
For `@generated` methods, set `enter_generated=true` if you want to extract the lowered code
of the generator itself.

On return `framecode` is the [`FrameCode`](@ref) of the method.
`frameargs` contains the actual arguments needed for executing this frame (for generators,
this will be the types of `allargs`);
`lenv` is the "environment", i.e., the static parameters for `f` given `allargs`.
`argtypes` is the `Tuple`-type for this specific call (equivalent to the signature of the `MethodInstance`).

# Example

```jldoctest
julia> mymethod(x::Vector{T}) where T = 1
mymethod (generic function with 1 method)

julia> framecode, frameargs, lenv, argtypes = JuliaInterpreter.prepare_call(mymethod, [mymethod, [1.0,2.0]]);

julia> framecode
  1  1  1 ─     return 1

julia> frameargs
2-element Array{Any,1}:
 mymethod (generic function with 1 method)
 [1.0, 2.0]

julia> lenv
svec(Float64)

julia> argtypes
Tuple{typeof(mymethod),Array{Float64,1}}
```
"""
function prepare_call(@nospecialize(f), allargs; enter_generated = false)
    # Can happen for thunks created by generated functions
    if isa(f, Core.Builtin) || isa(f, Core.IntrinsicFunction)
        return nothing
    elseif any(is_vararg_type, allargs)
        return nothing  # https://github.com/JuliaLang/julia/issues/30995
    end
    argtypesv = Any[_Typeof(a) for a in allargs]
    argtypes = Tuple{argtypesv...}
    method = whichtt(argtypes)
    if method === nothing
        # Call it to generate the exact error
        f(allargs[2:end]...)
    end
    ret = prepare_framecode(method, argtypes; enter_generated=enter_generated)
    # Exceptional returns
    if ret === nothing
        # The generator threw an error. Let's generate the same error by calling it.
        f(allargs[2:end]...)
    end
    isa(ret, Compiled) && return ret, argtypes
    # Typical return
    framecode, lenv = ret
    if is_generated(method) && enter_generated
        allargs = Any[_Typeof(a) for a in allargs]
    end
    return framecode, allargs, lenv, argtypes
end

function prepare_framedata(framecode, argvals::Vector{Any}, lenv::SimpleVector=empty_svec, caller_will_catch_err::Bool=false)
    src = framecode.src
    slotnames = src.slotnames::SlotNamesType
    ssavt = src.ssavaluetypes
    ng, ns = isa(ssavt, Int) ? ssavt : length(ssavt::Vector{Any}), length(src.slotflags)
    if length(junk_framedata) > 0
        olddata = pop!(junk_framedata)
        locals, ssavalues, sparams = olddata.locals, olddata.ssavalues, olddata.sparams
        exception_frames, last_reference = olddata.exception_frames, olddata.last_reference
        last_exception = olddata.last_exception
        callargs = olddata.callargs
        resize!(locals, ns)
        fill!(locals, nothing)
        resize!(ssavalues, 0)
        resize!(ssavalues, ng)
        # for check_isdefined to work properly, we need sparams to start out unassigned
        resize!(sparams, 0)
        empty!(exception_frames)
        resize!(last_reference, ns)
        last_exception[] = nothing
    else
        locals = Vector{Union{Nothing,Some{Any}}}(nothing, ns)
        ssavalues = Vector{Any}(undef, ng)
        sparams = Vector{Any}(undef, 0)
        exception_frames = Int[]
        last_reference = Vector{Int}(undef, ns)
        callargs = Any[]
        last_exception = Ref{Any}(nothing)
    end
    fill!(last_reference, 0)
    if isa(framecode.scope, Method)
        meth = framecode.scope::Method
        nargs, meth_nargs = length(argvals), Int(meth.nargs)
        islastva = meth.isva && nargs >= meth_nargs
        for i = 1:meth_nargs-islastva
            if nargs >= i
                locals[i], last_reference[i] = Some{Any}(argvals[i]), 1
            else
                locals[i] = Some{Any}(())
            end
        end
        if islastva
            locals[meth_nargs] =  (let i=meth_nargs; Some{Any}(ntupleany(k->argvals[i+k-1], nargs-i+1)); end)
            last_reference[meth_nargs] = 1
        end
    end
    resize!(sparams, length(lenv))
    # Add static parameters to environment
    for i = 1:length(lenv)
        T = lenv[i]
        isa(T, TypeVar) && continue  # only fill concrete types
        sparams[i] = T
    end
    FrameData(locals, ssavalues, sparams, exception_frames, last_exception, caller_will_catch_err, last_reference, callargs)
end

"""
    frame = prepare_frame(framecode::FrameCode, frameargs, lenv)

Construct a new `Frame` for `framecode`, given lowered-code arguments `frameargs` and
static parameters `lenv`. See [`JuliaInterpreter.prepare_call`](@ref) for information about how to prepare the inputs.
"""
function prepare_frame(framecode::FrameCode, args::Vector{Any}, lenv::SimpleVector, caller_will_catch_err::Bool=false)
    framedata = prepare_framedata(framecode, args, lenv, caller_will_catch_err)
    return Frame(framecode, framedata)
end

function prepare_frame_caller(caller::Frame, framecode::FrameCode, args::Vector{Any}, lenv::SimpleVector)
    caller_will_catch_err = !isempty(caller.framedata.exception_frames) || caller.framedata.caller_will_catch_err
    caller.callee = frame = prepare_frame(framecode, args, lenv, caller_will_catch_err)
    frame.caller = caller
    return frame
end

"""
    ExprSplitter(mod::Module, ex::Expr; lnn=nothing)

Create an iterable that returns individual expressions together with their module of evaluation.
Optionally supply an initial `LineNumberNode` `lnn`.

# Example

```
julia> expr = quote
           public(x::Integer) = true
           module Private
           private(y::String) = false
           end
           const threshold = 0.1
       end;

julia> for (mod, ex) in ExprSplitter(Main, expr)
           @show mod ex
       end
mod = Main
ex = quote
    #= REPL[7]:2 =#
    public(x::Integer) = begin
            #= REPL[7]:2 =#
            true
        end
end
mod = Main.Private
ex = quote
    #= REPL[7]:4 =#
    private(y::String) = begin
            #= REPL[7]:4 =#
            false
        end
end
mod = Main
ex = :($(Expr(:toplevel, :(#= REPL[7]:6 =#), :(const threshold = 0.1))))
```

Note that `Main.Private` was created for you so that its internal expressions could be evaluated.
`ExprSplitter` will check to see whether the module already exists and if so return it rather than
try to create a new module with the same name.

In general each returned expression is a block with two parts: a `LineNumberNode` followed by a single expression.
In some cases the returned expression may be `:toplevel`, as shown in the `const` declaration,
but otherwise it will be a `:block`.

# World age, frame creation, and evaluation

The primary purpose of `ExprSplitter` is to allow sequential return to top-level (e.g., the REPL)
after evaluation of each expression. Returning to top-level allows the world age to update, and hence allows one to call
methods and use types defined in earlier expressions in a block.

For evaluation by JuliaInterpreter, the returned module/expression pairs can be passed directly to
the `Frame` constructor. However, some expressions cannot be converted into `Frame`s and may need
special handling:

```julia
julia> for (mod, ex) in ExprSplitter(Main, expr)
           if ex.head === :global
               # global declarations can't be lowered to a CodeInfo.
               # In this demo we choose to evaluate them, but you can do something else.
               Core.eval(mod, ex)
               continue
           end
           frame = Frame(mod, ex)
           debug_command(frame, :c, true)
       end

julia> threshold
0.1

julia> public(3)
true
```

If you're parsing package code, `ex` might be a docstring-expression; you may wish
to check for such expressions and take distinct actions.

See [`Frame(mod::Module, ex::Expr)`](@ref) for more information about frame creation.
"""
mutable struct ExprSplitter
    # Non-mutating fields
    stack::Vector{Tuple{Module,Expr}}   # mod[i] is module of evaluation for
    index::Vector{Int}    # next-to-handle argument index for :block or :toplevel exprs
    # Mutating fields
    lnn::Union{LineNumberNode,Nothing}
end
function ExprSplitter(mod::Module, ex::Expr; lnn=nothing)
    iter = ExprSplitter(Tuple{Module,Expr}[], Int[], lnn)
    push_modex!(iter, mod, ex)
    queuenext!(iter)
    return iter
end

Base.IteratorSize(::Type{ExprSplitter}) = Base.SizeUnknown()
Base.eltype(::Type{ExprSplitter}) = Tuple{Module,Expr}

function push_modex!(iter::ExprSplitter, mod::Module, ex::Expr)
    push!(iter.stack, (mod, ex))
    if ex.head === :toplevel || ex.head === :block
        # Issue #427
        modifies_scope = false
        for a in ex.args
            if isa(a, Expr) && a.head ∈ (:local, :global)
                modifies_scope = true
                break
            end
        end
        push!(iter.index, modifies_scope ? 0 : 1)
    end
    return iter
end

function pop_modex!(iter)
    mod, ex = pop!(iter.stack)
    if ex.head === :toplevel || ex.head === :block
        pop!(iter.index)
    end
    return mod, ex
end

# Load the next-to-evaluate expression into `iter.stack[end]`.
function queuenext!(iter::ExprSplitter)
    isempty(iter.stack) && return nothing
    mod, ex = iter.stack[end]
    head = ex.head
    if head === :module
        # Find or create the module
        newname = ex.args[2]::Symbol
        if isdefined(mod, newname)
            newmod = getfield(mod, newname)
            newmod isa Module || throw(ErrorException("invalid redefinition of constant $(newname)"))
            mod = newmod
        else
            newnamestr = String(newname)
            id = Base.identify_package(mod, newnamestr)
            # If we're in a test environment and Julia's internal stdlibs are not a declared dependency of the package,
            # we might fail to find it. Try really hard to find it.
            if id === nothing && mod === Base.__toplevel__
                for loaded_id in keys(Base.loaded_modules)
                    if loaded_id.name == newnamestr
                        id = loaded_id
                        break
                    end
                end
            end
            if id !== nothing && haskey(Base.loaded_modules, id)
                mod = Base.root_module(id)::Module
            else
                loc = firstline(ex)
                mod = Core.eval(mod, Expr(:module, ex.args[1], ex.args[2], Expr(:block, loc, loc)))::Module
            end
        end
        # We've handled the module declaration, remove it and queue the body
        pop!(iter.stack)
        ex = ex.args[3]::Expr
        push_modex!(iter, mod, ex)
        return queuenext!(iter)
    elseif VERSION < v"1.2" && is_function_def(ex)
        # Newer Julia versions insert a LineNumberNode between statements, but at least Julia 1.0 does not.
        # Scan the function body for a LNN
        lnn = firstline(ex)
        if isa(lnn, LineNumberNode)
            if iter.lnn === nothing || iter.lnn.line < lnn.line
                iter.lnn = LineNumberNode(lnn.line-1, lnn.file)
            end
        end
    elseif head === :macrocall
        iter.lnn = ex.args[2]::LineNumberNode
    elseif head === :block || head === :toplevel
        # Container expression
        idx = iter.index[end]
        if idx == 0
            # return the whole block (issue #427)
            return nothing
        end
        while idx <= length(ex.args)
            a = ex.args[idx]
            if isa(a, LineNumberNode)
                iter.lnn = a
            elseif isa(a, Expr)
                iter.index[end] = idx + 1
                push_modex!(iter, mod, a)
                return queuenext!(iter)
            end
            idx += 1
        end
        # We exhausted the expression without returning anything to evaluate
        pop!(iter.stack)
        pop!(iter.index)
        return queuenext!(iter)
    end
    return nothing  # mod, ex will be returned by iterate
end

function Base.iterate(iter::ExprSplitter, state=nothing)
    isempty(iter.stack) && return nothing
    mod, ex = pop_modex!(iter)
    lnn = iter.lnn
    if is_doc_expr(ex)
        body = ex.args[4]
        if isa(body, Expr) && body.head === :module
            # Just document the module itself and push the module def onto the stack
            excopy = Expr(ex.head, ex.args[1], ex.args[2], ex.args[3])
            push!(excopy.args, body.args[2])
            append!(excopy.args, ex.args[5:end])   # there should only be at most a 5th, but just for robustness
            ex = excopy
            push_modex!(iter, mod, body)
        end
    end
    if ex.head === :block || ex.head === :toplevel
        # This was a block that we couldn't safely descend into (issue #427)
        if !isempty(iter.index) && iter.index[end] > length(iter.stack[end][2].args)
            pop!(iter.stack)
            pop!(iter.index)
            queuenext!(iter)
        end
        return (mod, ex), nothing
    end
    queuenext!(iter)
    # :global expressions can't be lowered. For debugging it might be nice
    # to still return the lnn, but then we have to work harder on detecting them.
    ex.head === :global && return (mod, ex), nothing
    return (mod, Expr(:block, lnn, ex)), nothing
end

"""
    framecode, frameargs, lenv, argtypes = determine_method_for_expr(expr; enter_generated = false)

Prepare all the information needed to execute a particular `:call` expression `expr`.
For example, try `JuliaInterpreter.determine_method_for_expr(:(\$sum([1,2])))`.
See [`JuliaInterpreter.prepare_call`](@ref) for information about the outputs.
"""
function determine_method_for_expr(expr; enter_generated = false)
    f = to_function(expr.args[1])
    allargs = expr.args
    # Extract keyword args
    kwargs = Expr(:parameters)
    if length(allargs) > 1 && isexpr(allargs[2], :parameters)
        kwargs = splice!(allargs, 2)::Expr
    end
    f, allargs = prepare_args(f, allargs, kwargs.args)
    return prepare_call(f, allargs; enter_generated=enter_generated)
end

"""
    frame = enter_call_expr(expr; enter_generated=false)

Build a `Frame` ready to execute the expression `expr`. Set `enter_generated=true`
if you want to execute the generator of a `@generated` function, rather than the code that
would be created by the generator.

# Example

```jldoctest
julia> mymethod(x) = x+1
mymethod (generic function with 1 method)

julia> JuliaInterpreter.enter_call_expr(:(\$mymethod(1)))
Frame for mymethod(x) in Main at none:1
  1* 1  1 ─ %1 = (+)(x, 1)
  2  1  └──      return %1
x = 1

julia> mymethod(x::Vector{T}) where T = 1
mymethod (generic function with 2 methods)

julia> a = [1.0, 2.0]
2-element Array{Float64,1}:
 1.0
 2.0

julia> JuliaInterpreter.enter_call_expr(:(\$mymethod(\$a)))
Frame for mymethod(x::Array{T,1}) where T in Main at none:1
  1* 1  1 ─     return 1
x = [1.0, 2.0]
T = Float64
```

See [`enter_call`](@ref) for a similar approach not based on expressions.
"""
function enter_call_expr(expr; enter_generated = false)
    clear_caches()
    r = determine_method_for_expr(expr; enter_generated = enter_generated)
    if r !== nothing && !isa(r[1], Compiled)
        return prepare_frame(Base.front(r)...)
    end
    nothing
end

"""
    frame = enter_call(f, args...; kwargs...)

Build a `Frame` ready to execute `f` with the specified positional and keyword arguments.

# Example

```jldoctest
julia> mymethod(x) = x+1
mymethod (generic function with 1 method)

julia> JuliaInterpreter.enter_call(mymethod, 1)
Frame for mymethod(x) in Main at none:1
  1* 1  1 ─ %1 = (+)(x, 1)
  2  1  └──      return %1
x = 1

julia> mymethod(x::Vector{T}) where T = 1
mymethod (generic function with 2 methods)

julia> JuliaInterpreter.enter_call(mymethod, [1.0, 2.0])
Frame for mymethod(x::Array{T,1}) where T in Main at none:1
  1* 1  1 ─     return 1
x = [1.0, 2.0]
T = Float64
```

For a `@generated` function you can use `enter_call((f, true), args...; kwargs...)`
to execute the generator of a `@generated` function, rather than the code that
would be created by the generator.

See [`enter_call_expr`](@ref) for a similar approach based on expressions.
"""
function enter_call(@nospecialize(finfo), @nospecialize(args...); kwargs...)
    clear_caches()
    if isa(finfo, Tuple)
        f = finfo[1]
        enter_generated = finfo[2]::Bool
    else
        f = finfo
        enter_generated = false
    end
    f, allargs = prepare_args(f, Any[f, args...], kwargs)
    # Can happen for thunks created by generated functions
    if isa(f, Core.Builtin) || isa(f, Core.IntrinsicFunction)
        error(f, " is a builtin or intrinsic")
    end
    r = prepare_call(f, allargs; enter_generated=enter_generated)
    if r !== nothing && !isa(r[1], Compiled)
        return prepare_frame(Base.front(r)...)
    end
    return nothing
end

# This is a version of InteractiveUtils.gen_call_with_extracted_types, except that is passes back the
# call expression for further processing.
function extract_args(__module__, ex0)
    if isa(ex0, Expr)
        if any(a->(isexpr(a, :kw) || isexpr(a, :parameters)), ex0.args)
            arg1, args, kwargs = gensym("arg1"), gensym("args"), gensym("kwargs")
            return quote
                $arg1 = $(ex0.args[1])
                $args, $kwargs = $separate_kwargs($(ex0.args[2:end]...))
                tuple(Core.kwfunc($arg1), $kwargs, $arg1, $args...)
            end
        elseif ex0.head === :.
            return Expr(:tuple, :getproperty, ex0.args...)
        elseif ex0.head === :(<:)
            return Expr(:tuple, :(<:), ex0.args...)
        else
            return Expr(:tuple,
                mapany(x->isexpr(x,:parameters) ? QuoteNode(x) : x, ex0.args)...)
        end
    end
    if isexpr(ex0, :macrocall) # Make @edit @time 1+2 edit the macro by using the types of the *expressions*
        return error("Macros are not supported in @enter")
    end
    ex = Meta.lower(__module__, ex0)
    if !isa(ex, Expr)
        return error("expression is not a function call or symbol")
    elseif ex.head === :call
        return Expr(:tuple,
            mapany(x->isexpr(x, :parameters) ? QuoteNode(x) : x, ex.args)...)
    elseif ex.head === :body
        a1 = ex.args[1]
        if isexpr(a1, :call)
            a11 = a1.args[1]
            if a11 === :setindex!
                return Expr(:tuple,
                    mapany(x->isexpr(x, :parameters) ? QuoteNode(x) : x, arg.args)...)
            end
        end
    end
    return error("expression is not a function call, "
               * "or is too complex for @enter to analyze; "
               * "break it down to simpler parts if possible")
end

"""
    @interpret f(args; kwargs...)

Evaluate `f` on the specified arguments using the interpreter.

# Example

```jldoctest
julia> a = [1, 7]
2-element Array{Int64,1}:
 1
 7

julia> sum(a)
8

julia> @interpret sum(a)
8
```
"""
macro interpret(arg)
    args = try
        extract_args(__module__, arg)
    catch e
        return :(throw($e))
    end
    quote
        local theargs = $(esc(args))
        local frame = JuliaInterpreter.enter_call_expr(Expr(:call, theargs...))
        if frame === nothing
            eval(Expr(:call, map(QuoteNode, theargs)...))
        elseif shouldbreak(frame, 1)
            frame, BreakpointRef(frame.framecode, 1)
        else
            local ret = finish_and_return!(frame)
            # We deliberately return the top frame here; future debugging commands
            # via debug_command may alter the leaves, we want the top frame so we can
            # ultimately do `get_return`.
            isa(ret, BreakpointRef) ? (frame, ret) : ret
        end
    end
end