text:
module m_stokes_driven_tht6_veclap
using Test
using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities
using MeshSteward: T6block, T6toT3
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs
using Elfel.FEIterators: jacjac
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
# using UnicodePlots

mu = 0.25 # dynamic viscosity
A = 1.0 # length of the side of the square
N = 100;# number of subdivisions along the sides of the square domain

function genmesh()
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(A, A, N, N), "velocity")
    # This mesh will be used for the pressures. Notice that it needs to be
        # "compatible" with the velocity mesh in the sense that they need to share
        # the nodes at the corners of the triangles.
        pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    function integrateK!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(pel, pqp)
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * (2 * gradNux[i][1] * gradNux[j][1] + gradNux[i][2] * gradNux[j][2])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * (gradNuy[i][1] * gradNuy[j][1] + 2 * gradNuy[i][2] * gradNuy[j][2])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(uxfesp), FEIterator(uyfesp), FEIterator(pfesp))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(uxfesp, qargs), QPIterator(uyfesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrateK!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function test()
    vmesh, pmesh = genmesh()
    # Velocity spaces
    uxfesp = FESpace(Float64, vmesh, FEH1_T6(), 1)
    uyfesp = FESpace(Float64, vmesh, FEH1_T6(), 1)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # Part of the boundary that is immovable
    boxes = [[0.0 A 0.0 0.0], [0.0 0.0 0.0 A], [A A 0.0 A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(uxfesp, 0, i, 1, 0.0)
            setebc!(uyfesp, 0, i, 1, 0.0)
        end
    end
    # The lid: driven in the X direction
    uxbar = 1.0
    box = [0.0 A A A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(uxfesp, 0, i, 1, uxbar)
        setebc!(uyfesp, 0, i, 1, 0.0)
    end
    # Pressure space
    pfesp = FESpace(Float64, pmesh, FEH1_T3(), 1)
    setebc!(pfesp, 0, 1, 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([uxfesp, uyfesp, pfesp])
    tndof = ndofs(uxfesp) + ndofs(uyfesp) + ndofs(pfesp)
    tnunk = nunknowns(uxfesp) + nunknowns(uyfesp) + nunknowns(pfesp)
    @test (tndof, tnunk) == (91003, 89402)
    # Assemble the coefficient matrix
    K = assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [uxfesp, uyfesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([uxfesp, uyfesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(uxfesp, "ux", 1)
    makeattribute(uyfesp, "uy", 1)
    vtkwrite("m_stokes_driven_tht6_veclap-p", baseincrel(pmesh), [(name = "p",), ])
    vtkwrite("m_stokes_driven_tht6_veclap-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    try rm("m_stokes_driven_tht6_veclap-p.vtu"); catch end
    try rm("m_stokes_driven_tht6_veclap-v.vtu"); catch end
    true
end

end
using .m_stokes_driven_tht6_veclap
m_stokes_driven_tht6_veclap.test()

module m_stokes_driven_t3b

using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities, ir_identity
using MeshSteward: T3block
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T3_BUBBLE, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs
using Elfel.FEIterators: jacjac
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
# using UnicodePlots

mu = 0.25 # dynamic viscosity
A = 1.0 # length of the side of the square
N = 100;# number of subdivisions along the sides of the square domain

function genmesh()
    # This mesh will be both for the velocities and for the pressure
    mesh = Mesh()
    attach!(mesh, T3block(A, A, N, N), "velocity+pressure")
    ir = baseincrel(mesh)
    eidir = ir_identity(ir)
    attach!(mesh, eidir)
    return mesh
end

function assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    function integrateK!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxuy = LocalMatrixAssembler(uxnedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxuy, eldofs(uxel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(pel, pqp)
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * (2 * gradNux[i][1] * gradNux[j][1] + gradNux[i][2] * gradNux[j][2])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * (gradNuy[i][1] * gradNuy[j][1] + 2 * gradNuy[i][2] * gradNuy[j][2])
                end
                for j in 1:uynedof, i in 1:uxnedof
                    kuxuy[i, j] += (mu * JxW) * (gradNux[i][1] * gradNuy[j][2])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuxuy)
            assemble!(ass, transpose(kuxuy))
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(uxfesp), FEIterator(uyfesp), FEIterator(pfesp))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(uxfesp, qargs), QPIterator(uyfesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrateK!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function test()
    mesh = genmesh()
    # Velocity spaces
    uxfesp = FESpace(Float64, mesh, FEH1_T3_BUBBLE(), 1)
    uyfesp = FESpace(Float64, mesh, FEH1_T3_BUBBLE(), 1)
    locs = geometry(mesh)
    inflate = A / N / 100
    # Part of the boundary that is immovable
    boxes = [[0.0 A 0.0 0.0], [0.0 0.0 0.0 A], [A A 0.0 A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(uxfesp, 0, i, 1, 0.0)
            setebc!(uyfesp, 0, i, 1, 0.0)
        end
    end
    # The lid
    uxbar = 1.0
    box = [0.0 A A A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(uxfesp, 0, i, 1, uxbar)
        setebc!(uyfesp, 0, i, 1, 0.0)
    end
    # Pressure space
    pfesp = FESpace(Float64, mesh, FEH1_T3(), 1)
    setebc!(pfesp, 0, 1, 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([uxfesp, uyfesp, pfesp])
    tndof = ndofs(uxfesp) + ndofs(uyfesp) + ndofs(pfesp)
    tnunk = nunknowns(uxfesp) + nunknowns(uyfesp) + nunknowns(pfesp)
    # Assemble the coefficient matrix
    K = assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [uxfesp, uyfesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([uxfesp, uyfesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(uxfesp, "ux", 1)
    makeattribute(uyfesp, "uy", 1)
    vtkwrite("m_stokes_driven_t3b-p", baseincrel(mesh), [(name = "p",), ])
    vtkwrite("m_stokes_driven_t3b-v", baseincrel(mesh), [(name = "ux",), (name = "uy",)])
    try rm("m_stokes_driven_t3b-p.vtu"); catch end
    try rm("m_stokes_driven_t3b-v.vtu"); catch end
    true
end

end
using .m_stokes_driven_t3b
m_stokes_driven_t3b.test()

module mstok1

using Test

"""
    th_p2_p1

The manufactured-solution colliding flow example from Elman et al 2014. The
Hood-Taylor formulation with quadratic triangles for the velocity and continuous
pressure on linear triangles.

The formulation is the one derived in Reddy, Introduction to the finite element
method, 1993. Page 486 ff.
"""
module th_p2_p1

using Test
using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities, ir_identity, attribute, VecAttrib
using MeshSteward: T6block, T6toT3
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary, transform
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs, eldofvals
using Elfel.FEIterators: jacjac, location
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(2 * A, 2 * A, N, N), "velocity")
    ir = baseincrel(vmesh)
    transform(ir, x -> x .- A)
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    function integrateK!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxuy = LocalMatrixAssembler(uxnedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxuy, eldofs(uxel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(pel, pqp)
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * (2 * gradNux[i][1] * gradNux[j][1] + gradNux[i][2] * gradNux[j][2])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * (gradNuy[i][1] * gradNuy[j][1] + 2 * gradNuy[i][2] * gradNuy[j][2])
                end
                for j in 1:uynedof, i in 1:uxnedof
                    kuxuy[i, j] += (mu * JxW) * (gradNux[i][1] * gradNuy[j][2])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuxuy)
            assemble!(ass, transpose(kuxuy))
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(uxfesp), FEIterator(uyfesp), FEIterator(pfesp))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(uxfesp, qargs), QPIterator(uyfesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrateK!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(pfesp)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(pfesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(pfesp, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(uxfesp, uyfesp)
    function integrate!(elits, qpits, trueux, trueuy)
        uxnedof, uynedof = ndofsperel.(elits)
        E = 0.0
        for el in zip(elits...)
            uxel, uyel = el
            uxdofvals = eldofvals(uxel)
            uydofvals = eldofvals(uyel)
            for qp in zip(qpits...)
                uxqp, uyqp = qp
                Jac, J = jacjac(uxel, uxqp)
                JxW = J * weight(uxqp)
                Np = bfun(uxqp)
                uxt = trueux(location(uxel, uxqp)...)
                uyt = trueuy(location(uyel, uyqp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:uxnedof
                    uxa += (uxdofvals[j] * Np[j])
                    uya += (uydofvals[j] * Np[j])
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elits = (FEIterator(uxfesp), FEIterator(uyfesp),)
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(uxfesp, qargs), QPIterator(uyfesp, qargs),)
    return integrate!(elits, qpits, trueux, trueuy)
end

function run(N)
    vmesh, pmesh = genmesh(N)
    # Velocity spaces
    uxfesp = FESpace(Float64, vmesh, FEH1_T6(), 1)
    uyfesp = FESpace(Float64, vmesh, FEH1_T6(), 1)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(uxfesp, 0, i, 1, trueux(locs[i]...))
            setebc!(uyfesp, 0, i, 1, trueuy(locs[i]...))
        end
    end
    # Pressure space
    pfesp = FESpace(Float64, pmesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(pmesh); nearestto = [0.0, 0.0])
    setebc!(pfesp, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([uxfesp, uyfesp, pfesp])
    tndof = ndofs(uxfesp) + ndofs(uyfesp) + ndofs(pfesp)
    tnunk = nunknowns(uxfesp) + nunknowns(uyfesp) + nunknowns(pfesp)
    # Assemble the coefficient matrix
    K = assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [uxfesp, uyfesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([uxfesp, uyfesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(uxfesp, "ux", 1)
    makeattribute(uyfesp, "uy", 1)
    ep = evaluate_pressure_error(pfesp)
    ev = evaluate_velocity_error(uxfesp, uyfesp)
    # vtkwrite("th_p2_p1-p", baseincrel(pmesh), [(name = "p",), ])
    # vtkwrite("th_p2_p1-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    # geom = geometry(pfesp.mesh)
    # ir = baseincrel(pfesp.mesh)
    # pt = VecAttrib([truep(geom[i]...) for i in 1:length(geom)])
    # ir.right.attributes["pt"] = pt
    # vtkwrite("th_p2_p1-pt", baseincrel(pmesh), [(name = "pt",), ])
    (ep, ev)
end

end

function test()
    ref = [(3.5171450671095306, 0.2968271617227661),
    (0.5999467323539439, 0.03781189670123018),
    (0.12350320261417459, 0.004741849976722882)]
    N = 4
    for loop in 1:3
        ep, ev = th_p2_p1.run(N)
        @test isapprox(vec([ep, ev]), vec([ref[loop]...]))
        N = N * 2
    end
end

end

using .mstok1
mstok1.test()

"""
    m_th_p2_p1_veclap_alt

The manufactured-solution colliding flow example from Elman et al 2014. The
Hood-Taylor formulation with quadratic triangles for the velocity and continuous
pressure on linear triangles.

This implementation is an alternative: for the velocity, a single finite element
space with multiple components (2) is used instead of multiple finite element
spaces.

The formulation is the one derived in Donea, Huerta, Introduction to the finite element
method, 1993. Page 486 ff, and Elman, et al., Finite elements and fast
iterative solvers, p. 132. In other words, it is the vector Laplacian version.
"""
module m_th_p2_p1_veclap_alt

using Test
using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities, ir_identity, attribute, VecAttrib
using MeshSteward: T6block, T6toT3
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary, transform
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs, eldofvals
using Elfel.FEIterators: jacjac, location
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(2 * A, 2 * A, N, N), "velocity")
    ir = baseincrel(vmesh)
    transform(ir, x -> x .- A)
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(ufesp, pfesp, tndof, mu)
    function integrate!(ass, elits, qpits, mu)
        unedof, pnedof = ndofsperel.(elits)
        uedofcomp = edofcompnt(ufesp)
        kuu = LocalMatrixAssembler(unedof, unedof, 0.0)
        kup = LocalMatrixAssembler(unedof, pnedof, 0.0)
        for el in zip(elits...)
            uel, pel = el
            init!(kuu, eldofs(uel), eldofs(uel))
            init!(kup, eldofs(uel), eldofs(pel))
            for qp in zip(qpits...)
                uqp, pqp = qp
                Jac, J = jacjac(uel, uqp)
                JxW = J * weight(uqp)
                gradNu = bfungrad(uqp, Jac)
                Np = bfun(pqp)
                for j in 1:unedof, i in 1:unedof
                    if uedofcomp[i] == uedofcomp[j]
                        kuu[i, j] += (mu * JxW) * (dot(gradNu[i], gradNu[j]))
                    end
                end
                for j in 1:pnedof, i in 1:unedof
                    kup[i, j] += (-JxW * Np[j]) * gradNu[i][uedofcomp[i]]
                end
            end
            assemble!(ass, kuu)
            assemble!(ass, kup)
            assemble!(ass, transpose(kup))
        end
        return ass
    end

    elits = (FEIterator(ufesp), FEIterator(pfesp))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(ufesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrate!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(pfesp)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(pfesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(pfesp, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(ufesp)
    function integrate!(elit, qpit, trueux, trueuy)
        unedof = ndofsperel(elit)
        uedofcomp = edofcompnt(ufesp)
        E = 0.0
        for el in elit
            udofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Nu = bfun(qp)
                uxt = trueux(location(el, qp)...)
                uyt = trueuy(location(el, qp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:unedof
                    (uedofcomp[j] == 1) && (uxa += (udofvals[j] * Nu[j]))
                    (uedofcomp[j] == 2) && (uya += (udofvals[j] * Nu[j]))
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(ufesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(ufesp, qargs)
    return integrate!(elit, qpit, trueux, trueuy)
end

function run(N)
    vmesh, pmesh = genmesh(N)
    # Velocity spaces
    ufesp = FESpace(Float64, vmesh, FEH1_T6(), 2)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(ufesp, 0, i, 1, trueux(locs[i]...))
            setebc!(ufesp, 0, i, 2, trueuy(locs[i]...))
        end
    end
    # Pressure space
    pfesp = FESpace(Float64, pmesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(pmesh); nearestto = [0.0, 0.0])
    setebc!(pfesp, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([ufesp, pfesp])
    tndof = ndofs(ufesp) + ndofs(pfesp)
    tnunk = nunknowns(ufesp) + nunknowns(pfesp)
    # Assemble the coefficient matrix
    K = assembleK(ufesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [ufesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([ufesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(ufesp, "ux", 1)
    makeattribute(ufesp, "uy", 2)
    ep = evaluate_pressure_error(pfesp)
    ev = evaluate_velocity_error(ufesp)
    # vtkwrite("m_th_p2_p1_veclap_alt-p", baseincrel(pmesh), [(name = "p",), ])
    # vtkwrite("m_th_p2_p1_veclap_alt-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    return (ep, ev)
end

test() = begin
    ep, ev = m_th_p2_p1_veclap_alt.run(4)
    @test isapprox([ep, ev], [2.596076907594511, 0.3001331486426876])
end

end

using .m_th_p2_p1_veclap_alt
m_th_p2_p1_veclap_alt.test()

module m_th_p2_p1_veclap_alt_e

using Test
using LinearAlgebra
using StaticArrays
using MeshCore.Exports
using MeshSteward.Exports
using Elfel.Exports
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(2 * A, 2 * A, N, N), "velocity")
    ir = baseincrel(vmesh)
    transform(ir, x -> x .- A)
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(ufesp, pfesp, tndof, mu)
    function integrate!(ass, elits, qpits, mu)
        unedof, pnedof = ndofsperel.(elits)
        uedofcomp = edofcompnt(ufesp)
        kuu = LocalMatrixAssembler(unedof, unedof, 0.0)
        kup = LocalMatrixAssembler(unedof, pnedof, 0.0)
        for el in zip(elits...)
            uel, pel = el
            init!(kuu, eldofs(uel), eldofs(uel))
            init!(kup, eldofs(uel), eldofs(pel))
            for qp in zip(qpits...)
                uqp, pqp = qp
                Jac, J = jacjac(uel, uqp)
                JxW = J * weight(uqp)
                gradNu = bfungrad(uqp, Jac)
                Np = bfun(pqp)
                for j in 1:unedof, i in 1:unedof
                    if uedofcomp[i] == uedofcomp[j]
                        kuu[i, j] += (mu * JxW) * (dot(gradNu[i], gradNu[j]))
                    end
                end
                for j in 1:pnedof, i in 1:unedof
                    kup[i, j] += (-JxW * Np[j]) * gradNu[i][uedofcomp[i]]
                end
            end
            assemble!(ass, kuu)
            assemble!(ass, kup)
            assemble!(ass, transpose(kup))
        end
        return ass
    end

    elits = (FEIterator(ufesp), FEIterator(pfesp))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(ufesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrate!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(pfesp)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(pfesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(pfesp, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(ufesp)
    function integrate!(elit, qpit, trueux, trueuy)
        unedof = ndofsperel(elit)
        uedofcomp = edofcompnt(ufesp)
        E = 0.0
        for el in elit
            udofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Nu = bfun(qp)
                uxt = trueux(location(el, qp)...)
                uyt = trueuy(location(el, qp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:unedof
                    (uedofcomp[j] == 1) && (uxa += (udofvals[j] * Nu[j]))
                    (uedofcomp[j] == 2) && (uya += (udofvals[j] * Nu[j]))
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(ufesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(ufesp, qargs)
    return integrate!(elit, qpit, trueux, trueuy)
end

function run(N)
    vmesh, pmesh = genmesh(N)
    # Velocity spaces
    ufesp = FESpace(Float64, vmesh, FEH1_T6(), 2)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(ufesp, 0, i, 1, trueux(locs[i]...))
            setebc!(ufesp, 0, i, 2, trueuy(locs[i]...))
        end
    end
    # Pressure space
    pfesp = FESpace(Float64, pmesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(pmesh); nearestto = [0.0, 0.0])
    setebc!(pfesp, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([ufesp, pfesp])
    tndof = ndofs(ufesp) + ndofs(pfesp)
    tnunk = nunknowns(ufesp) + nunknowns(pfesp)
    # Assemble the coefficient matrix
    K = assembleK(ufesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [ufesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([ufesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(ufesp, "ux", 1)
    makeattribute(ufesp, "uy", 2)
    ep = evaluate_pressure_error(pfesp)
    ev = evaluate_velocity_error(ufesp)
    # vtkwrite("m_th_p2_p1_veclap_alt_e-p", baseincrel(pmesh), [(name = "p",), ])
    # vtkwrite("m_th_p2_p1_veclap_alt_e-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    return (ep, ev)
end

test() = begin
    ep, ev = m_th_p2_p1_veclap_alt_e.run(4)
    @test isapprox([ep, ev], [2.596076907594511, 0.3001331486426876])
end

end

using .m_th_p2_p1_veclap_alt_e
m_th_p2_p1_veclap_alt_e.test()


text:
@doc raw"""
```
t = find_threshold(Yen(), histogram, edges)
```

Computes the threshold value using Yen's maximum correlation criterion for
bilevel thresholding.

# Output

Returns a real number `t` in `edges`. The `edges` parameter represents an
`AbstractRange` which specifies the intervals associated with the histogram bins.

# Details

This algorithm uses the concept of *entropic correlation* of a gray level histogram to produce a threshold
value.

Let ``f_1, f_2, \ldots, f_I`` be the frequencies in the various bins of the
histogram and ``I`` the number of bins. With ``N = \sum_{i=1}^{I}f_i``, let
``p_i = \frac{f_i}{N}`` (``i = 1, \ldots, I``) denote the probability
distribution of gray levels. From this distribution one derives two additional
distributions. The first defined for discrete values ``1`` to ``s`` and the
other, from ``s+1`` to ``I``. These distributions are

```math
A: \frac{p_1}{P_s}, \frac{p_2}{P_s}, \ldots, \frac{p_s}{P_s}
\quad \text{and} \quad
B: \frac{p_{s+1}}{1-P_s}, \ldots, \frac{p_n}{1-P_s}
\quad \text{where} \quad
P_s = \sum_{i=1}^{s}p_i.
```
The entropic correlations associated with each distribution are

```math
C(A) = -\ln \sum_{i=1}^{s} \left( \frac{p_i}{P_s} \right)^2 \quad \text{and} \quad C(B) = -\ln \sum_{i=s+1}^{I} \left( \frac{p_i}{1 - P_s} \right)^2.
```

Combining these two entropic correlation functions we have

```math
\psi(s) = -\ln \sum_{i=1}^{s} \left( \frac{p_i}{P_s} \right)^2 -\ln \sum_{i=s+1}^{I} \left( \frac{p_i}{1 - P_s} \right)^2.
```
Finding the discrete value ``s`` which maximises the function ``\psi(s)`` produces
the sought-after threshold value (i.e. the bin which determines the threshold).

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia

using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img, 256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Yen(), counts[1:end], edges)
```

# Reference

1. Yen JC, Chang FJ, Chang S (1995), “A New Criterion for Automatic Multilevel Thresholding”, IEEE Trans. on Image Processing 4 (3): 370-378, [doi:10.1109/83.366472](https://doi.org/10.1109/83.366472)
"""
function find_threshold(algorithm::Yen, histogram::AbstractArray, edges::AbstractRange)
    total = sum(histogram)
    m = length(histogram)
    p = zeros(m)

    # Calulate probability.
    for i in eachindex(p)
        p[i] = histogram[i] / total
    end

    # Setup sums.
    Pₛ = 0
    Gₛ = 0
    G′ₛ = sum(p.^2)
    maxval = typemin(Float64)
    t = firstindex(edges)
    for s in eachindex(p[1:m-1])
        # Update sums.
        Pₛ += p[s]
        Gₛ += p[s]^2
        G′ₛ -= p[s]^2

        # Calculate total correlation.
        TC = -log(Gₛ * G′ₛ) + 2 * log(Pₛ * (1 - Pₛ))
        if TC > maxval
            maxval = TC
            t = s
        end
    end

    return edges[t]
end

ext: jl


text:
#! format: off

using SIIPExamples
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder
using DataFrames

using HiGHS # mip solver
solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.05)

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

print_tree(PowerSimulations.PM.AbstractPowerModel)

PTDF_matrix = PTDF(sys)

template = template_unit_commitment(
    network = NetworkModel(
        StandardPTDFModel,
        PTDF = PTDF_matrix,
        duals = [CopperPlateBalanceConstraint],
        use_slacks = false,
    ),
    use_slacks = true,
)
for (k, v) in template.branches
    v.duals = [NetworkFlowConstraint]
end

problem = DecisionModel(template, sys, horizon = 24, optimizer = solver)
build!(problem, output_dir = mktempdir())

solve!(problem)

res = ProblemResults(problem)
duals = Dict([
    k => read_dual(res, k) for
    k in list_dual_keys(res) if PSI.get_entry_type(k) == NetworkFlowConstraint
])
λ = read_dual(res, "CopperPlateBalanceConstraint__System")[:, 2]
flow_duals = outerjoin(values(duals)..., on = :DateTime)
μ = Matrix(flow_duals[:, PTDF_matrix.axes[1]])

LMP = flow_duals[:, [:DateTime]]
for bus in get_components(Bus, sys)
    LMP[:, get_name(bus)] = λ .+ μ * PTDF_matrix[:, get_number(bus)]
end

LMP

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
@testset "ExponentialLikelihood" begin
    lik = ExponentialLikelihood()
    test_interface(lik, SqExponentialKernel(), exp.(randn(10)))
end


text:
# Mixing rules: (ρin, ρout) => ρnext, where ρout is produced by diagonalizing the Hamiltonian at ρin
# These define the basic fix-point iteration, that are then combined with acceleration methods (eg anderson)
# All these methods attempt to approximate the inverse jacobian of the SCF step, J^-1 = (1 - χ0 vc)^-1
# Note that "mixing" is sometimes used to refer to the combined process of formulating the fixed-point
# and solving it; we call "mixing" only the first part

# The interface is `mix(m, basis, ρin, ρout) -> ρnext`

using LinearMaps
using IterativeSolvers

@doc raw"""
Kerker mixing: ``J^{-1} ≈ \frac{α G^2}{(k_F^2 + G^2}``
where ``k_F`` is the Thomas-Fermi screening constant.

Notes:
  - Abinit calls ``1/k_F`` the dielectric screening length (parameter *dielng*)
"""
struct KerkerMixing
    α::Real
    kF::Real
end
# Default parameters suggested by Kresse, Furthmüller 1996 (α=0.8, kF=1.5 Ǎ^{-1})
# DOI 10.1103/PhysRevB.54.11169
KerkerMixing(;α=0.8, kF=0.8) = KerkerMixing(α, kF)
function mix(mixing::KerkerMixing, basis, ρin::RealFourierArray, ρout::RealFourierArray; kwargs...)
    T = eltype(basis)
    Gsq = [sum(abs2, basis.model.recip_lattice * G) for G in G_vectors(basis)]
    ρin = ρin.fourier
    ρout = ρout.fourier
    ρnext = @. ρin + T(mixing.α) * (ρout - ρin) * Gsq / (T(mixing.kF)^2 + Gsq)
    # take the correct DC component from ρout; otherwise the DC component never gets updated
    ρnext[1] = ρout[1]
    from_fourier(basis, ρnext; assume_real=true)
end

"""
Simple mixing: J^-1 ≈ α
"""
struct SimpleMixing
    α::Real
end
SimpleMixing(;α=1) = SimpleMixing(α)
function mix(mixing::SimpleMixing, basis, ρin::RealFourierArray, ρout::RealFourierArray; kwargs...)
    if mixing.α == 1
        return ρout
    else
        T = eltype(basis)
        ρin + T(mixing.α) * (ρout - ρin)
    end
end


text:
using StochasticDiffEq, DiffEqDevTools, Test, LinearAlgebra, Random

Random.seed!(100)
const σ_const = 0.87
const μ = 1.01

u0 = rand(2)
A = [-2.0 1.0;1.0 -2.0]
B = Diagonal([σ_const for i in 1:2])

function f_nondiag(u,p,t)
  return A*u + 1.01*u
end

function f_nondiag_iip(du,u,p,t)
  mul!(du,A,u)
  du .+= 1.01u
end

function f_analytic_nondiag(u0,p,t,W)
 tmp = (A+1.01I-(B^2))*t + B*sum(W)
 exp(tmp)*u0
end

function g_nondiag(u,p,t)
  du = zeros(2,2)
  du[1,1] = σ_const*u[1]
  du[1,2] = σ_const*u[1]
  du[2,1] = σ_const*u[2]
  du[2,2] = σ_const*u[2]
  return du
end

function g_nondiag_iip(du,u,p,t)
  du[1,1] = σ_const*u[1]
  du[1,2] = σ_const*u[1]
  du[2,1] = σ_const*u[2]
  du[2,2] = σ_const*u[2]
end

coeff = 2*σ_const^2 #To not compute the same coefficient in the sde.
function ggprime(u, p, t)
  return coeff*u
end

function ggprime(du, u, p, t)
  du .= coeff*u
end

prob = SDEProblem(SDEFunction(f_nondiag,g_nondiag,analytic=f_analytic_nondiag),g_nondiag,u0,(0.0,1.0),noise_rate_prototype=zeros(2,2))
probiip = SDEProblem(SDEFunction(f_nondiag_iip,g_nondiag_iip,analytic=f_analytic_nondiag),g_nondiag_iip,u0,(0.0,1.0),noise_rate_prototype=zeros(2,2))

## Just solve to test compatibility
IEM = solve(probiip,ImplicitEM())
IEM = solve(probiip,ImplicitEulerHeun())
IEM = solve(probiip,ISSEM())
IEM = solve(probiip,ISSEulerHeun())

## use BenchmarkTools to benchmark speed here if desired.
dt = 0.002
seed = 1
# solEMiip = solve(probiip, EM(), dt=dt, seed=seed)
# solRKMil = solve(probiip, RKMilCommute(), dt=dt, seed=seed)
solPCEuler = solve(prob, PCEuler(ggprime), dt=dt, seed=seed)
solPCEuleriip = solve(probiip, PCEuler(ggprime), dt=dt, seed=seed)

@test solPCEuler.u ≈ solPCEuleriip.u atol=1e-10

## plot to see performance of PCEuler
# first_elem(x) = x[1]
# using PyPlot
# figure()
# plot(solEMiip.t, first_elem.(solEMiip.u_analytic), label="True", "--")
# plot(solEMiip.t, first_elem.(solEMiip.u), label="EM")
# plot(solRKMil.t, first_elem.(solRKMil.u), label="RKMil")
# plot(solPCEuleriip.t, first_elem.(solPCEuleriip.u), label="PC")
# legend()

##
dts = (1/2) .^ (10:-1:5) #14->7 good plot
trajectories = 50
simEM = test_convergence(dts,probiip,EM(),trajectories=trajectories)
simPCEuler = test_convergence(dts,probiip,PCEuler(ggprime),trajectories=trajectories)
#simRKMil = test_convergence(dts,probiip,RKMilCommute(),trajectories=trajectories)
@test all(simPCEuler.errors[:l2] .< simEM.errors[:l2])

## Plotting script to see the order 1 scaling of PCEuler
# figure()
# loglog(dts, simEM.errors[:l2], label="EM")
# loglog(dts, simRKMil.errors[:l2], label="RKMil")
# loglog(dts, simPCEuler.errors[:l2], label="PC")
# ylabel("l2 error")
# xlabel("dt")
# legend()

f_morenoise = (du,u,p,t) -> du.=1.01u
g_morenoise = function (du,u,p,t)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f_morenoise,g_morenoise,ones(2),(0.0,1.0),
                  noise_rate_prototype=zeros(2,4))

sol =solve(prob,dt=1/2^(3),EM())
sol =solve(prob,dt=1/2^(3),ISSEM())
sol =solve(prob,dt=1/2^(3),ImplicitEM())
sol =solve(prob,dt=1/2^(3),EulerHeun())
sol =solve(prob,dt=1/2^(3),ImplicitEulerHeun())
sol =solve(prob,dt=1/2^(3),ISSEulerHeun())


text:
Pkg.activate("./docs");

using Documenter, CollegeStratCollege, FilesLH

makedocs(
    modules = [CollegeStratCollege],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "hendri54",
    sitename = "CollegeStratCollege.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

pkgDir = rstrip(normpath(@__DIR__, ".."), '/');
@assert endswith(pkgDir, "CollegeStratCollege")
deploy_docs(pkgDir; trialRun = false);

Pkg.activate(".");

# ------------


text:
module ccp

#include("WishboneDevice.jl")
#include("TransportTapcp.jl")

end # module

max_stars_repo_name: jkocz/ccp

max_issues_repo_name: jkocz/ccp

max_forks_repo_name: jkocz/ccp


text:
"""
`SimplifiedString(∇V, integrate!, reparameterize!, dist, Δt)` -  Set up the
simplified string method data structure

### Fields

* `∇V`   - Gradient of the potential
* `integrate!` - Integrator for ẋ = - ∇V(x)
* `reparameterize!` - Choice of reparametrization
* `dist` - Distance function used in reparametrization and convergence testing
* `pin` - Boolean for pinning the endpoints of hte string
* `Δt`    - Time step for the integrator
"""
struct SimplifiedString{TGV, TI, TR, TD, TP<:Bool, TF<:AbstractFloat} <: SimplifiedStringMethod
    ∇V::TGV
    integrate!::TI
    reparameterize!::TR
    dist::TD
    pin::TP
    Δt::TF
end

"""
`simplified_string` - Run the simplified string method for an energy landscape
set up with the `SimplifiedString` data structure.  This can return the entire
string evolution.

### Fields
* `U` - Initial string
* `S` - Simplified string data structure
  ### Optional Fields
* `options` - String method options
"""
function simplified_string(U₀, S::TS; options=StringOptions()) where {TS <: SimplifiedString}
    U = deepcopy(U₀)
    U_new = deepcopy(U);

    U_trajectory = typeof(U)[deepcopy(U)];

    n_images = length(U);
    Δs = 1/(n_images-1);
    err_est = 0.;

    for n in 1:options.nmax
        # time stepping routine

        if S.pin
            Threads.@threads for j in 2:n_images-1
                S.integrate!(U_new[j], S.∇V, S.Δt);
            end
        else
            Threads.@threads for j in 1:n_images
                S.integrate!(U_new[j], S.∇V, S.Δt);
            end
        end
        # reparametrization step
        S.reparameterize!(U_new, S.dist, S.pin)

        err_est = maximum([S.dist(U_new[i], U[i]) for i in 1:n_images])/S.Δt;
        if(options.verbose && mod(n, options.print_iters)==0)
            @printf("[%d]: error = %g\n", n, err_est);
        end

        copy!.(U, U_new);
        if(options.save_trajectory)
            push!(U_trajectory, deepcopy(U));
        end
        if (err_est< options.tol)
            if(!options.save_trajectory)
                push!(U_trajectory, deepcopy(U));
            end
            break
        end
    end

    if(err_est >= options.tol)
        @printf("ERROR: Did not converge after %d iterartions", options.nmax)
    end
    return U_trajectory
end

"""
`simplified_string!` - Run the simplified string method for an energy landscape
set up with the `SimplifiedString` data structure.  This is an in place transform.

### Fields
* `U` - Initial string
* `S` - Simplified string data structure
  ### Optional Fields
* `options` - String method options
"""
function simplified_string!(U, S::TS; options=StringOptions()) where {TS <: SimplifiedString}
    U_new = deepcopy(U);

    n_images = length(U);
    Δs = 1/(n_images-1);
    err_est = 0.;

    for n in 1:options.nmax
        # time stepping routine
        if S.pin
            Threads.@threads for j in 2:n_images-1
                S.integrate!(U_new[j], S.∇V, S.Δt);
            end
        else
            Threads.@threads for j in 1:n_images
                S.integrate!(U_new[j], S.∇V, S.Δt);
            end
        end
        # reparametrization step
        S.reparameterize!(U_new, S.dist, S.pin)

        err_est = maximum([S.dist(U_new[i], U[i]) for i in 1:n_images])/S.Δt;
        if(options.verbose && mod(n, options.print_iters)==0)
            @printf("[%d]: error = %g\n", n, err_est);
        end

        copy!.(U, U_new);
        if (err_est< options.tol)
            break
        end
    end

    if(err_est >= options.tol)
        @printf("ERROR: Did not converge after %d iterartions", options.nmax)
    end
    U
end


text:
using ScoreDrivenERGM
using Test


text:
function grammar_vladislavleva_c(; erc=false)
    Base.eval(Main, :(protectedDiv(x, y) = Benchmarks.protectedDiv(x, y)))
    Base.eval(Main, :(protectedSin(x) = Benchmarks.protectedSin(x)))
    Base.eval(Main, :(protectedCos(x) = Benchmarks.protectedCos(x)))
    grammar = erc ?
        (@grammar begin
            R = x
            R = y
            R = R + R
            R = R - R
            R = R * R
            R = protectedDiv(R,R)
            R = R^2
            R = exp(R)
            R = exp(-R)
            R = protectedSin(R)
            R = protectedCos(R)
            R = abs(R)^E
            R = R + E
            R = R * E
            E = _(5*rand())
        end) :
        (@grammar begin
            R = x
            R = y
            R = R + R
            R = R - R
            R = R * R
            R = protectedDiv(R,R)
            R = R^2
            R = exp(R)
            R = exp(-R)
            R = protectedSin(R)
            R = protectedCos(R)
            R = abs(R)^E
            R = R + E
            R = R * E
            E = |(-5.0:1.0:5.0)
        end)
    grammar
end

gt_vladislavleva_3(x,y) = exp(-x)*x^3*(protectedCos(x)*protectedSin(x))*(protectedCos(x)*protectedSin(x)^2-1)*(y-5)
function loss_vladislavleva_3(S::SymbolTable, tree::RuleNode, grammar::Grammar)
    ex = get_executable(tree, grammar)
    los = 0.0
    n = 0
    for x in 0.05:0.1:10.0, y in 0.05:2.0:10.05
        S[:x] = x
        S[:y] = y
        los += abs2(Core.eval(S,ex) - gt_vladislavleva_3(x,y))
        n += 1
    end
    los / n
end

function run_vladislavleva_3(seed::Int=0, n_pop::Int=300)
    Random.seed!(seed)

    grammar = grammar_vladislavleva_c(erc=true)
    S = SymbolTable(grammar, ExprOptimization.Benchmarks)
    loss(tree::RuleNode, grammar::Grammar) = loss_vladislavleva_3(S, tree, grammar)

    init_method = RandomInit()
    select_method = TournamentSelection(3)
    p = GeneticProgram(n_pop, #pop_size
                       40,  #iterations
                       10,  #max_depth
                       0.4, #p_reproduction
                       0.5, #p_crossover
                       0.1, #p_mutation
                       init_method=init_method,
                       select_method=select_method)
    tstart = CPUtime_us()
    result = optimize(p, grammar, :R, loss)
    dt = CPUtime_us() - tstart
    result, dt
end


text:
using ITensors
using Test

@testset "Test debug checks on IndexSet construction" begin
  i = Index(2, "i")

  @test !ITensors.using_debug_checks()
  # Test that no error is thrown in constructor
  @test IndexSet(i, i) isa IndexSet{2}
  @test IndexSet(i, i') isa IndexSet{2}

  # Turn on debug checks
  ITensors.enable_debug_checks()
  @test ITensors.using_debug_checks()
  @test_throws ErrorException IndexSet(i, i)
  # Test that no error is thrown in constructor
  @test IndexSet(i, i') isa IndexSet{2}

  # Turn off debug checks
  ITensors.disable_debug_checks()
  @test !ITensors.using_debug_checks()
  # Test that no error is thrown in constructor
  @test IndexSet(i, i) isa IndexSet{2}
  @test IndexSet(i, i') isa IndexSet{2}
end


text:
# helper function to compare NaN
nearlysame(x, y) = x ≈ y || (isnan(x) & isnan(y))
nearlysame(A::AbstractArray, B::AbstractArray) = all(map(nearlysame, A, B))

# AbstractGray and Color3 tests should be generated seperately
function generate_test_types(number_types::AbstractArray{<:DataType}, color_types::AbstractArray{<:UnionAll})
    test_types = map(Iterators.product(number_types, color_types)) do T
        try
            T[2]{T[1]}
        catch err
            !isa(err, TypeError) && rethrow(err)
        end
    end
    test_types = filter(x->x != false, test_types)
    if isempty(filter(x->x<:Color3, test_types))
        test_types = [number_types..., test_types...]
    end
    test_types
end

"""
    test_colwise(dist, n, sz, T)

test if `colwise` works as expected in the following two inputs:

* `(n,)Vector` - `(n,)Vector`
* `(1,n)Matrix` - `(1,n)Matrix`

and throw error for ambigious inputs:

* `(m,n)Matrix` - `(m,b)Matrix`  where `m!=1`

"""
function test_colwise(dist, n, sz, T)
    @testset "colwise" begin
        vec_imgsA = [rand(T, sz) for _ in 1:n]
        vec_imgsB = [rand(T, sz) for _ in 1:n]
        mat_imgsA = reshape(vec_imgsA, (1, n))
        mat_imgsB = reshape(vec_imgsB, (1, n))
        mat_imgsA_wrong = reshape(vec_imgsA, (n, 1))
        mat_imgsB_wrong = reshape(vec_imgsB, (n, 1))

        RT = result_type(dist, vec_imgsA, vec_imgsB)
        r1 = zeros(RT, n)
        r2 = zeros(RT, n)
        for j = 1:n
            r1[j] = dist(vec_imgsA[j], vec_imgsB[j])
            r2[j] = dist(mat_imgsA[1,j], mat_imgsB[1,j])
        end
        @test nearlysame(colwise(dist, vec_imgsA, vec_imgsB), r1)
        @test nearlysame(colwise(dist, mat_imgsA, mat_imgsB), r2)
        @test_throws DimensionMismatch colwise(dist, mat_imgsA_wrong, mat_imgsB_wrong)
    end
end

"""
    test_pairwise(dist, nA, nB, sz, T)

test if `pairwise` works as expected in the following two inputs:

* `(n,)Vector`-`(n,)Vector`
"""
function test_pairwise(dist, nA, nB, sz, T)
    @testset "pairwise" begin
        vec_imgsA = [rand(T, sz) for _ in 1:nA]
        vec_imgsB = [rand(T, sz) for _ in 1:nB]

        RT = result_type(dist, vec_imgsA, vec_imgsB)
        rAB = zeros(RT, nA, nB)
        rAA = zeros(RT, nA, nA)
        for j = 1:nB, i = 1:nA
            rAB[i, j] = dist(vec_imgsA[i], vec_imgsB[j])
        end
        for j = 1:nA, i = 1:nA
            rAA[i, j] = dist(vec_imgsA[i], vec_imgsA[j])
        end
        @test nearlysame(pairwise(dist, vec_imgsA, vec_imgsB), rAB)
        @test nearlysame(pairwise(dist, vec_imgsA, vec_imgsA), rAA)
        @test nearlysame(pairwise(dist, vec_imgsA), rAA)
    end
end

_base_colorant_type(::Type{<:Number}) = Gray
_base_colorant_type(::Type{T}) where T<:Colorant = base_colorant_type(T)
"""
    test_numeric(dist, a, b, T; filename=nothing)

simply test that `dist` works for 2d image as expected, more tests go to `Distances.jl`
"""
function test_numeric(dist, a, b, T; filename=nothing)
    size(a) == size(b) || error("a and b should be the same size")
    if filename == nothing
        filename = "references/$(typeof(dist))_$(ndims(a))d"

        if eltype(a) <: Color3
            filename = filename * "_$(_base_colorant_type(T))"
        elseif eltype(a) <: Union{Number, AbstractGray}
            filename = filename * "_$(_base_colorant_type(T))"
        end
    end
    @testset "numeric" begin
        @testset "$T" begin
            @test_reference "$(filename).txt" Float64(dist(a, b)) by=(x,y)->isapprox(x,y; atol=1e-8)
        end
    end
end

"""
    test_cross_type(dist, a, b, type_list)

simply test if operations between `N0f8`, `Bool` and `Float32` types works as expected.
`a` and `b` should be simple enough to get rid of InexactError.
"""
function test_cross_type(dist, a, b, type_list)
    size(a) == size(b) || error("a and b should be the same size")
    rsts = [[dist(Ta.(a), Tb.(b)),
             dist(Tb.(a), Ta.(b))] for (Ta, Tb) in subsets(type_list, 2)]
    rsts = hcat(rsts...)
    @test all(isapprox.(rsts, rsts[1]; rtol=1e-5))
end

function test_Metric(d, sz, T)
    x = rand(T, sz)
    y = rand(T, sz)
    z = rand(T, sz)
    @test d(x, y) >= 0
    @test d(x, x) ≈ 0
    @test d(x, y) ≈ d(y, x)
    @test d(x, y) + d(y, z) >= d(x, z)
end

function test_SemiMetric(d, sz, T)
    x = rand(T, sz)
    y = rand(T, sz)
    z = rand(T, sz)
    @test d(x, y) >= 0
    @test d(x, x) ≈ 0
    @test d(x, y) ≈ d(y, x)
end

function test_ndarray(d, sz, T)
    x = rand(T, sz)
    y = rand(T, sz)
    @test_nowarn d(x, y)

    T <: AbstractGray || return nothing
    @test d(x, y) ≈ d(channelview(x), channelview(y))
end


text:
module Antennas

# Top level imports
using DimensionalData
using Unitful

# Top level types
"""
    RadiationPattern
Holds the pattern data, isa AbstractDimArray
"""
struct RadiationPattern{T<:Number,
    N,
    D<:Tuple,
    R<:Tuple,
    A<:AbstractArray{T,N},
    Na,
    Me} <: AbstractDimArray{T,N,D,A}
data::A
dims::D
refdims::R
name::Na
metadata::Me
end

"""
    AntennaArray
Holds the antenna array data
"""
struct AntennaArray
    excitations
    locations
end

include("Patterns.jl")
include("Arrays.jl")
include("Plotting.jl")

end


text:
#
#    extractor(formula_def::String)
#
# Extract the formula name and the variables from `formula_def`.
#
function extractor(formula_def)
    try
        ast = Meta.parse(string(formula_def))
        #dump(ast)
        formula::Formula = Formula(formula_def, ast)
        find_symbol(ast, formula, 0)
        return formula
    catch ex
        throw(WrongFormula(formula_def))
    end
end

# Look up symbol and return value, or throw.
function find_symbol(s::Symbol, formula::Formula, level)
    push!(formula.inputs, String(s))
end

function find_symbol(x::Number, formula::Formula, level)
    # do nothing
end

# To parse an expression, convert the head to a singleton
# tag, so that Julia can dispatch on that tag.
function find_symbol(e::Expr, formula::Formula, exprlevel)

    find_symbol(Val(e.head), e.args, formula, exprlevel+1)
end

# Call the function named in args[1]
function find_symbol(::Val{:call}, args, formula::Formula, level)
    op = args[1]
    if hasproperty(Statistics, op) || op in [:sum]
        # A progressive has a statistic operator at level 1
        formula.progressive = true
    end

    find_symbol(Val(args[1]), args[2:end], formula, level)
end

function find_symbol(::Val{:.}, args, formula::Formula, level)
    qn = args[2].value
    push!(formula.inputs, "$(args[1]).$qn")
end

# formula definition
function find_symbol(::Val{:(=)}, args, formula::Formula, level)
    # formula name
    formula.output = String(args[1])
    formula.expr = args[2]
    find_symbol(args[2], formula, level)
end

function find_symbol(op::Val, args, formula::Formula, level)
    for arg ∈ args
        find_symbol(arg, formula, level)
    end
end


text:
abstract type DistanceTransform end

struct Chamfer{T} <: DistanceTransform
    dt
end

"""
    Chamfer(img, dt=zeros(Float32, size(img)))

Prepares an array to be `transform`ed using the 3-4 chamfer algorithm
laid out in 'Distance transformations in digital images,
Computer Vision, Graphics, and Image Processing'
[Gunilla Borgefors](https://studentportalen.uu.se/uusp-filearea-tool/download.action?nodeId=214320&toolAttachmentId=64777)

# Arguments
- img: 2D or 3D array to be transformed based on location
    to the nearest background (0) pixel
- tfm: `zeros(Float32, size(img))`
"""
Chamfer(img, dt=zeros(Float32, size(img))) = Chamfer{typeof(dt)}(dt)

"""
    transform(img::AbstractMatrix{T}, tfm::Chamfer)
    transform(img::AbstractArray{T,3}, tfm::Chamfer)

Applies a 3-4 chamfer distance transform to an input image.
Returns an array with spatial information embedded in the
array elements.

# Arguments
- img: 2D or 3D array to be transformed based on location
    to the nearest background (0) pixel

# Citation
'Distance transformations in digital images,
Computer Vision, Graphics, and Image Processing'
[Gunilla Borgefors](DOI: https://doi.org/10.1016/S0734-189X(86)80047-0.)
"""
function transform(img::AbstractMatrix{T}, tfm::Chamfer) where {T}
    dt = tfm.dt
    w, h = size(img)
    # Forward pass
    x = 1
    y = 1
    if img[x, y] == 0
        dt[x, y] = 65535 # some large value
    end
    for x in 1:(w - 1)
        if img[x + 1, y] == 0
            dt[x + 1, y] = 3 + dt[x, y]
        end
    end
    for y in 1:(h - 1)
        x = 1
        if img[x, y + 1] == 0
            dt[x, y + 1] = min(3 + dt[x, y], 4 + dt[x + 1, y])
        end
        for x in 1:(w - 2)
            if img[x + 1, y + 1] == 0
                dt[x + 1, y + 1] = min(
                    4 + dt[x, y], 3 + dt[x + 1, y], 4 + dt[x + 2, y], 3 + dt[x, y + 1]
                )
            end
        end
        x = w

        if img[x, y + 1] == 0
            dt[x, y + 1] = min(4 + dt[x - 1, y], 3 + dt[x, y], 3 + dt[x - 1, y + 1])
        end
    end

    # Backward pass
    for x in (w - 1):-1:1
        y = h
        if img[x, y] == 0
            dt[x, y] = min(dt[x, y], 3 + dt[x + 1, y])
        end
    end
    for y in (h - 1):-1:1
        x = w
        if img[x, y] == 0
            dt[x, y] = min(dt[x, y], 3 + dt[x, y + 1], 4 + dt[x - 1, y + 1])
        end
        for x in 1:(w - 2)
            if img[x + 1, y] == 0
                dt[x + 1, y] = min(
                    dt[x + 1, y],
                    4 + dt[x + 2, y + 1],
                    3 + dt[x + 1, y + 1],
                    4 + dt[x, y + 1],
                    3 + dt[x + 2, y],
                )
            end
        end
        x = 1
        if img[x, y] == 0
            dt[x, y] = min(
                dt[x, y], 4 + dt[x + 1, y + 1], 3 + dt[x, y + 1], 3 + dt[x + 1, y]
            )
        end
    end
    return dt
end

# Helper function for nd arrays
function _transform(img::AbstractMatrix{T}, tfm::Chamfer, dt) where {T}
    w, h = size(img)
    # Forward pass
    x = 1
    y = 1
    if img[x, y] == 0
        dt[x, y] = 65535 # some large value
    end
    for x in 1:(w - 1)
        if img[x + 1, y] == 0
            dt[x + 1, y] = 3 + dt[x, y]
        end
    end
    for y in 1:(h - 1)
        x = 1
        if img[x, y + 1] == 0
            dt[x, y + 1] = min(3 + dt[x, y], 4 + dt[x + 1, y])
        end
        for x in 1:(w - 2)
            if img[x + 1, y + 1] == 0
                dt[x + 1, y + 1] = min(
                    4 + dt[x, y], 3 + dt[x + 1, y], 4 + dt[x + 2, y], 3 + dt[x, y + 1]
                )
            end
        end
        x = w

        if img[x, y + 1] == 0
            dt[x, y + 1] = min(4 + dt[x - 1, y], 3 + dt[x, y], 3 + dt[x - 1, y + 1])
        end
    end

    # Backward pass
    for x in (w - 1):-1:1
        y = h
        if img[x, y] == 0
            dt[x, y] = min(dt[x, y], 3 + dt[x + 1, y])
        end
    end
    for y in (h - 1):-1:1
        x = w
        if img[x, y] == 0
            dt[x, y] = min(dt[x, y], 3 + dt[x, y + 1], 4 + dt[x - 1, y + 1])
        end
        for x in 1:(w - 2)
            if img[x + 1, y] == 0
                dt[x + 1, y] = min(
                    dt[x + 1, y],
                    4 + dt[x + 2, y + 1],
                    3 + dt[x + 1, y + 1],
                    4 + dt[x, y + 1],
                    3 + dt[x + 2, y],
                )
            end
        end
        x = 1
        if img[x, y] == 0
            dt[x, y] = min(
                dt[x, y], 4 + dt[x + 1, y + 1], 3 + dt[x, y + 1], 3 + dt[x + 1, y]
            )
        end
    end
    return dt
end

function transform(img::AbstractArray{T,3}, tfm::Chamfer) where {T}
    dt = tfm.dt
    for z in 1:size(img)[3]
        dt[:, :, z] = _transform(img[:, :, z], tfm, dt[:, :, z])
    end
    return dt
end


text:
# code for loading .continuous files
### Constants ###
const CONT_REC_TIME_BITTYPE = Int64
const CONT_REC_N_SAMP = 1024
const CONT_REC_N_SAMP_BITTYPE = UInt16
const CONT_REC_REC_NO_BITTYPE = UInt16
const CONT_REC_SAMP_BITTYPE = Int16
const CONT_REC_BYTES_PER_SAMP = sizeof(CONT_REC_SAMP_BITTYPE)
const CONT_REC_END_MARKER = UInt8[0x00, 0x01, 0x02, 0x03, 0x04,
                                  0x05, 0x06, 0x07, 0x08, 0xff]
const CONT_REC_HEAD_SIZE = mapreduce(sizeof, +, [CONT_REC_TIME_BITTYPE,
                             CONT_REC_N_SAMP_BITTYPE, CONT_REC_REC_NO_BITTYPE])
const CONT_REC_BODY_SIZE = CONT_REC_N_SAMP * sizeof(CONT_REC_SAMP_BITTYPE)
const CONT_REC_TAIL_SIZE = sizeof(CONT_REC_END_MARKER)
const CONT_REC_BLOCK_SIZE = CONT_REC_HEAD_SIZE + CONT_REC_BODY_SIZE +
    CONT_REC_TAIL_SIZE

### Types ###
"Represents the header of each data block"
mutable struct BlockHeader
    timestamp::CONT_REC_TIME_BITTYPE
    nsample::CONT_REC_N_SAMP_BITTYPE
    recordingnumber::CONT_REC_REC_NO_BITTYPE
end
BlockHeader() = BlockHeader(0, 0, 0)

"Represents the entirety of a data block"
mutable struct DataBlock
    head::BlockHeader
    body::Vector{UInt8}
    data::Vector{CONT_REC_SAMP_BITTYPE}
    tail::Vector{UInt8}
    function DataBlock(
        head::BlockHeader,
        body::Vector{UInt8},
        data::Vector{CONT_REC_SAMP_BITTYPE},
        tail::Vector{UInt8}
    )
        length(body) == CONT_REC_BODY_SIZE ||
            error("body length must be ", CONT_REC_BODY_SIZE)
        length(data) == CONT_REC_N_SAMP ||
            error("data length must be ", CONT_REC_N_SAMP)
        length(tail) == CONT_REC_TAIL_SIZE ||
            error("data length must be ", CONT_REC_TAIL_SIZE)
        return new(head, body, data, tail)
    end
end
function DataBlock()
    head = BlockHeader()
    body = @compat Vector{UInt8}(undef, CONT_REC_BODY_SIZE)
    data = @compat Vector{CONT_REC_SAMP_BITTYPE}(undef, CONT_REC_N_SAMP)
    tail =  @compat Vector{UInt8}(undef, CONT_REC_TAIL_SIZE)
    DataBlock(head, body, data, tail)
end

"""
    ContinuousFile(io::IOStream)
Type for an open continuous file.

# Fields

**`io`** `IOStream` object.

**`filepath`** Path to underlying file, possibly empty

**`nsample`** number of samples in a file.

**`nblock`** number of data blocks in a file.

**`header`** [`OriginalHeader`](@ref) of the current file.
"""
struct ContinuousFile
    "IOStream for open continuous file"
    io::IOStream
    "Path to file, possibly empty"
    filepath::String
    "Number of samples in file"
    nsample::Int
    "Number of data blocks in file"
    nblock::Int
    "File header"
    header::OriginalHeader
end
function ContinuousFile(io::IOStream, filepath::AbstractString = "")
    header = OriginalHeader(io) # Read header
    nblock = count_blocks(io)
    nsample = count_data(nblock)
    return ContinuousFile(io, string(filepath), nsample, nblock, header)
end

"""
Abstract array for file-backed OpenEphys data.

All subtypes support a ready-only array interface and should
be constructable with a single IOStream argument.
"""
abstract type OEArray{T} <: AbstractArray{T, 1} end

setindex!(::OEArray, ::Int) = throw(ReadOnlyMemoryError())

"""
Abstract array for file-backed continuous OpenEphys data.

Will throw [`CorruptedException`](@ref) if the data file has
a corrupt [`OriginalHeader`](@ref), is not the correct size
for an `.continuous` file, or contains corrupt data blocks.

Subtype of abstract type [`OEArray`](@ref) are read only,
and have with the following fields:

# Fields

**`contfile`** [`ContinuousFile`](@ref) for the current file.

**`block`** buffer object for the data blocks in the file.

**`blockno`** the current block being access in the file.

**`check`** `Bool` to check each data block's validity.
"""
abstract type OEContArray{T} <: OEArray{T} end

### Stuff for code generation ###
sampletype = Real
timetype = Real
rectype = Integer
jointtype = Tuple{sampletype, timetype, rectype}
arraytypes = ((:SampleArray, sampletype, DataBlock, Float64),
              (:TimeArray, timetype, BlockHeader, Float64),
              (:RecNoArray, rectype, BlockHeader, Int),
              (:JointArray, jointtype, DataBlock, Tuple{Float64, Float64, Int}))
### Generate array datatypes ###
for (typename, typeparam, buffertype, defaulttype) = arraytypes
    @eval begin
        mutable struct $(typename){T<:$(typeparam)} <: OEContArray{T}
            contfile::ContinuousFile
            block::$(buffertype)
            blockno::UInt
            check::Bool
        end
        function $(typename)(
            ::Type{T}, contfile::ContinuousFile, check::Bool = true
        ) where {T}
            if check
                if ! check_filesize(contfile.io)
                    throw(CorruptedException(string(
                        "\nThe size of this file indicates that it cannot be well-formed.\n",
                        "This likely means that the last data block is missing samples.\n",
                        "HINT: To attempt access to the remaining contents of this file, use:\n",
                        $(typename),
                        "(T, io, false) to turn off this check."
                    )))
                end
            end
            block = $(buffertype)()
            return $(typename){T}(contfile, block, 0, check)
        end
        function $(typename)(
            ::Type{T},
            io::IOStream,
            check::Bool = true,
            filepath::AbstractString = ""
        ) where {T}
            return $(typename)(T, ContinuousFile(io, filepath), check)
        end
        function $(typename)(
            ::Type{T},
            filepath::AbstractString,
            check::Bool = true
        ) where {T}
            ior = open(filepath, "r")
            atexit(() -> close(ior))
            return $(typename)(T, ior, check, filepath)
        end
        function $(typename)(firstarg::Union{IO, AbstractString}, args...)
            return $(typename)($(defaulttype), firstarg, args...)
        end
    end
end

const arrayargs = "([type::Type{T},] file::Union{IO, String}, [check::Bool, filepath::String])"
const arraypreamble =
    "Subtype of [`OEContArray`](@ref) to provide file backed access to OpenEphys"
@doc """
    SampleArray$arrayargs
$arraypreamble sample values. If `type` is a floating
point type, then the sample value will be converted to voltage (in uV). Otherwise,
the sample values will remain the raw ADC integer readings.
""" SampleArray
@doc """
    TimeArray$arrayargs
$arraypreamble time stamps. If `type` is a floating
point type, then the time stamps will be converted to seconds. Otherwise,
the time stamp will be the sample number.
""" TimeArray
@doc """
    RecNoArray$arrayargs
$arraypreamble numbers.
""" RecNoArray
@doc """
    JointArray$arrayargs
$arraypreamble data. Returns a tuple of type `type`, whose
values represent `(samplevalue, timestamp, recordingnumber)`. For a description of
each, see [`SampleArray`](@ref), [`TimeArray`](@ref), and [`RecNoArray`](@ref),
respectively.
""" JointArray

### Array interface ###
length(A::OEContArray) = A.contfile.nsample

size(A::OEContArray) = (length(A),)

Base.IndexStyle(::Type{T}) where {T<:OEContArray} = IndexLinear()

@inline function getindex(A::OEContArray, i::Integer)
    @boundscheck checkbounds(A, i)
    _getindex(A, i)
end

@inline function _getindex(A::OEContArray, i::Integer)
    prepare_block!(A, i)
    relidx = sampno_to_offset(i)
    @inbounds data = block_data(A, relidx)
    return convert_data(A, data)
end

### Array helper functions ###
"Load data block if necessary"
@inline function prepare_block!(A::OEContArray, i::Integer)
    blockno = sampno_to_block(i)
    if blockno != A.blockno
        seek_to_block(A.contfile.io, blockno)
        goodread = read_into!(A.contfile.io, A.block, A.check)
        goodread || throw(CorruptedException("Data block $blockno is malformed"))
        A.blockno = blockno
    end
    nothing
end

"Move io to data block"
@inline function seek_to_block(io::IOStream, blockno::Integer)
    blockpos = block_start_pos(blockno)
    if blockpos != position(io)
        seek(io, blockpos)
    end
end

### location functions ###
@inline sampno_to_block(sampno::Integer) = fld(sampno - 1, CONT_REC_N_SAMP) + 1

@inline sampno_to_offset(sampno::Integer) = mod(sampno - 1, CONT_REC_N_SAMP) + 1

@inline block_start_pos(blockno::Integer) =
    (blockno - 1) * CONT_REC_BLOCK_SIZE + HEADER_N_BYTES

# Un-exported utility function for debugging
function pos_to_blockno(pos::Integer)
    if pos < 0
        throw(ArgumentError("Invalid position"))
    elseif pos < HEADER_N_BYTES
        blockno = 0
    else
        blockno = fld(pos - HEADER_N_BYTES, CONT_REC_BLOCK_SIZE) + 1
    end
    return blockno
end

blockno_to_start_sampno(blockno::Integer) = (blockno - 1) * CONT_REC_N_SAMP + 1

### File access and conversion ###
"Read file data block into data block buffer"
@inline function read_into!(io::IOStream, block::DataBlock, check::Bool = true)
    goodread = read_into!(io, block.head)
    goodread || return goodread
    ## Read the body
    nbytes = readbytes!(io, block.body, CONT_REC_BODY_SIZE)
    goodread = nbytes == CONT_REC_BODY_SIZE
    goodread || return goodread
    if check
        goodread = verify_tail!(io, block.tail)
    else
        skip(io, CONT_REC_TAIL_SIZE)
    end
    goodread && convert_block!(block)
    return goodread
end
"Read block header into header buffer"
@inline function read_into!(io::IOStream, head::BlockHeader)
    goodread = true
    try
        head.timestamp = read(io, CONT_REC_TIME_BITTYPE)
        head.nsample = read(io, CONT_REC_N_SAMP_BITTYPE)
        goodread = head.nsample == CONT_REC_N_SAMP
        if goodread
            head.recordingnumber = read(io, CONT_REC_REC_NO_BITTYPE)
        end
    catch exception
        if isa(exception, EOFError)
            goodread = false
        else
            rethrow(exception)
        end
    end
    return goodread
end
read_into!(io::IOStream, head::BlockHeader, ::Bool) = read_into!(io, head)

"Convert the wacky data format in OpenEphys continuous files"
@inline function convert_block!(block::DataBlock)
    ptr = Ptr{CONT_REC_SAMP_BITTYPE}(pointer(block.body))
    if ENDIAN_BOM == 0x04030201
        # Host is little endian: Always true for now
        # Correct for big endianness of this data block
        for idx in 1:CONT_REC_N_SAMP
            unsafe_store!(ptr, ntoh(unsafe_load(ptr, idx)), idx)
        end
    end
    @compat unsafe_copyto!(pointer(block.data), ptr, CONT_REC_N_SAMP)
    return block.data
end

### Methods to access data in buffer ###
@inline @propagate_inbounds block_data(A::SampleArray, rel_idx::Integer) =
    @inbounds return A.block.data[rel_idx]
@inline block_data(A::TimeArray, rel_idx::Integer) =
    A.block.timestamp + rel_idx - 1
@inline block_data(A::RecNoArray, ::Integer) = A.block.recordingnumber

@inline @propagate_inbounds function block_data(A::JointArray, rel_idx::Integer)
    sample = A.block.data[rel_idx]
    timestamp = A.block.head.timestamp + rel_idx - 1
    recno = A.block.head.recordingnumber
    return sample, timestamp, recno
end

### Methods to convert raw values into desired ones ##
convert_data(OE::A, data) where {A<:OEContArray} =
    convert_data(A, OE.contfile.header, data)
function convert_data(
    ::Type{SampleArray{T}}, H::OriginalHeader, data::Integer
) where {T<:AbstractFloat}
    return data * convert(T, H.bitvolts)
end
function convert_data(
    ::Type{SampleArray{T}}, ::OriginalHeader, data::Integer
) where {T<:Integer}
    return convert(T, data)
end
function convert_data(
    ::Type{TimeArray{T}}, H::OriginalHeader, data::Integer
) where {T<:AbstractFloat}
    return convert(T, (data - 1) / H.samplerate) # First sample is at time zero
end
function convert_data(
    ::Type{TimeArray{T}}, ::OriginalHeader, data::Integer
) where {T<:Integer}
    return convert(T, data)
end
function convert_data(
    ::Type{RecNoArray{T}}, ::OriginalHeader, data::Integer
) where {T}
    return convert(T, data)
end
function convert_data(
    ::Type{JointArray{Tuple{S,T,R}}}, H::OriginalHeader, data::Tuple
) where {S<:sampletype,T<:timetype,R<:rectype}
    samp = convert_data(SampleArray{S}, H, data[1])
    timestamp = convert_data(TimeArray{T}, H, data[2])
    recno = convert_data(RecNoArray{R}, H, data[3])
    return samp, timestamp, recno
end

### Verification methods ###
"Verify end of block marker"
function verify_tail!(io::IOStream, tail::Vector{UInt8})
    nbytes = readbytes!(io, tail, CONT_REC_TAIL_SIZE)
    goodread = nbytes == CONT_REC_TAIL_SIZE && tail == CONT_REC_END_MARKER
    return goodread
end

"Check that file could be comprised of header and complete data blocks"
function check_filesize(file::IOStream)
    filesizeok = rem(filesize(file) - HEADER_N_BYTES, CONT_REC_BLOCK_SIZE) == 0
end

### Utility methods ###
function count_blocks(file::IOStream)
    fsize = stat(file).size
    return div(fsize - HEADER_N_BYTES, CONT_REC_BLOCK_SIZE)
end

count_data(numblocks::Integer) = numblocks * CONT_REC_N_SAMP
count_data(file::IOStream) = count_data(count_blocks(file))


text:
module Indirect

end


text:
module Models

using PyCall
using ASE: ASECalculator

"""
`EMT()`

Creates an `ASECalculator` that uses `ase.calculators.emt` to compute
energy and forces. This is very slow and is only included for
demonstration purposes.
"""
function EMT()
   emt = pyimport("ase.calculators.emt")
   return ASECalculator(emt.EMT())
end

"""
`EAM(potfile)`

Creates an `ase.calculators.eam` calculator to compute
energy and forces. This is very slow and is only included for
demonstration purposes.
"""
function EAM(potfile)
   eam = pyimport("ase.calculators.eam")
   return ASECalculator( eam.EAM(potential=potfile) )
end

function QuippyPot(id; kwargs...)
   QuippyPotential = pyimport("quippy.potential").Potential
   return ASECalculator(QuippyPotential(id; kwargs...))
end

function QuippySW()
   QuippyPotential = pyimport("quippy.potential").Potential
   qpot = QuippyPotential("IP SW", param_str = """<SW_params n_types="3" label="PRB_31_plus_H_Ge">
   <comment> Stillinger and Weber, Phys. Rev. B  31 p 5262 (1984), extended for other elements.  Ge and Si-Ge from Ethier and Lewis '92 </comment>
   <per_type_data type="1" atomic_num="1" />
   <per_type_data type="2" atomic_num="14" />
   <per_type_data type="3" atomic_num="32" />
   <per_pair_data atnum_i="1" atnum_j="1" AA="0.0" BB="0.0"
         p="0" q="0" a="1.0" sigma="1.0" eps="0.0" />
   <per_pair_data atnum_i="1" atnum_j="14" AA="8.581214" BB="0.0327827"
         p="4" q="0" a="1.25" sigma="2.537884" eps="2.1672" />
   <per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
         p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />

   <!-- pair data for Ge and Si-Ge -->
   <per_pair_data atnum_i="32" atnum_j="32" AA="7.049556277" BB="0.6022245584"
         p="4" q="0" a="1.80" sigma="2.181" eps="1.92590365783410138248" />
   <per_pair_data atnum_i="14" atnum_j="32" AA="7.049556277" BB="0.6022245584"
         p="4" q="0" a="1.80" sigma="2.138" eps="2.04313391101890694121" />

   <!-- triplet terms: atnum_c is the center atom, neighbours j and k -->
   <per_triplet_data atnum_c="1"  atnum_j="1"  atnum_k="1"
         lambda="21.0" gamma="1.20" eps="2.1675" />
   <per_triplet_data atnum_c="1"  atnum_j="1"  atnum_k="14"
         lambda="21.0" gamma="1.20" eps="2.1675" />
   <per_triplet_data atnum_c="1"  atnum_j="14" atnum_k="14"
         lambda="21.0" gamma="1.20" eps="2.1675" />

   <per_triplet_data atnum_c="14" atnum_j="1"  atnum_k="1"
         lambda="21.0" gamma="1.20" eps="2.1675" />
   <per_triplet_data atnum_c="14" atnum_j="1"  atnum_k="14"
         lambda="21.0" gamma="1.20" eps="2.1675" />
   <per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
         lambda="21.0" gamma="1.20" eps="2.1675" />

   <!-- triplet data for Ge and Si-Ge -->
   <per_triplet_data atnum_c="14" atnum_j="14" atnum_k="32"
      lambda="23.1" gamma="1.20" eps="2.10444772465437788017" />
   <per_triplet_data atnum_c="14" atnum_j="32" atnum_k="32"
      lambda="25.5" gamma="1.20" eps="2.04326828917050691242" />
   <per_triplet_data atnum_c="32" atnum_j="14" atnum_k="14"
      lambda="25.5" gamma="1.20" eps="2.04326828917050691242" />
   <per_triplet_data atnum_c="32" atnum_j="14" atnum_k="32"
      lambda="28.1" gamma="1.20" eps="1.98396169354838709676" />
   <per_triplet_data atnum_c="32" atnum_j="32" atnum_k="32"
      lambda="31.0" gamma="1.20" eps="1.92590365783410138247" />

   </SW_params>""")
   return ASECalculator(qpot)
end

function has_quippy()
   try
      pyimport("quippy")
      return true
   catch
      return false
   end
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1ae29a38-24fb-11eb-33b0-9d171e0cbcd0
begin
	using Revise, Images, ImageView, PhysicalOptics, FFTW, Plots
end

# ╔═╡ 2e643bb6-24fb-11eb-20e7-cf0175edcdb2
md"## Point Spread Function 4f system
Here we simulate the point spread function of a 4f-system. We use the following qualitative procedure:
* Propagate over distance f to the lens
* Multiply with lens transmission
* Propagate another f
* cut out the pupil
* Propagate another f to the second lens
* Multiply with lens transmission
* Propagate to the image plane

Luckily there is a numerical procedure to calculate the following steps with one Fourier transform:
* Propagate over distance f to the lens
* Multiply with lens transmission
* Propagate another f

Therefore we do only:
* Fourier Transform
* cut out the pupil
* Fourier Transform
"

# ╔═╡ 2e48f63a-24fb-11eb-2e00-47522330d613
begin
	L = 10e-3
	N = 2048
	f = 100e-3
	radius = 0.5e-3
	E0 = zeros((N, N))
	E0[1025, 1025] = 1
	E1, L1 = lens_propagate(E0, L, f)
	@show L1
	E1_ = circ(E1, radius, L1)
	E2, L2 = lens_propagate(E1_, L1, f)
	psf_ = abs2.(E2)
	psf = psf_ ./ sum(psf_)
end

# ╔═╡ 23584178-259a-11eb-13c5-9db519daef75
md"## Compare with analytical solution"

# ╔═╡ 2ca65ee2-2531-11eb-07c0-275d6a39d8cf
begin
	psf_analytical = jinc_psf((N, N), L, radius, f; shift=true)
end

# ╔═╡ 13283706-2534-11eb-19c9-e74e50186bae
begin
	psf1 = psf ./ maximum(psf)
	psf2 = psf_analytical ./ maximum(psf_analytical)
	imshow(cat(psf1, psf2, psf1 .- psf2, dims=3))
end

# ╔═╡ 8051b950-25a0-11eb-2e2c-f7986397d24e
begin
	xs = fftpos(L, N)
	x1 = 2 * N ÷5
	x2 = 3 * N ÷5
	slice = (N ÷2 + 1, x1:x2)
	plot(xs[x1:x2], psf1[slice...])
	plot!(xs[x1:x2], psf2[slice...])
end

# ╔═╡ Cell order:
# ╠═1ae29a38-24fb-11eb-33b0-9d171e0cbcd0
# ╟─2e643bb6-24fb-11eb-20e7-cf0175edcdb2
# ╠═2e48f63a-24fb-11eb-2e00-47522330d613
# ╟─23584178-259a-11eb-13c5-9db519daef75
# ╠═2ca65ee2-2531-11eb-07c0-275d6a39d8cf
# ╠═13283706-2534-11eb-19c9-e74e50186bae
# ╠═8051b950-25a0-11eb-2e2c-f7986397d24e

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
module SmartShockFinder

using BSON: @save
using TimerOutputs
using PrettyTables
using UnPack: @unpack
using Distributions: Uniform
using Downloads: download
using Flux
using Flux: onecold, crossentropy, @epochs, Data.DataLoader, params
using HDF5: h5open
using OffsetArrays: OffsetArray, OffsetVector
using StaticArrays: SVector, push
using Statistics: mean
using Random: shuffle!
using Trixi: TreeMesh, leaf_cells, length_at_cell, vandermonde_legendre, multiply_dimensionwise!, legendre_polynomial_and_derivative, n_directions
using Trixi: has_neighbor, has_children, polynomial_interpolation_matrix, has_any_neighbor, gauss_lobatto_nodes_weights

include("network.jl")
include("traindata/traindata.jl")

include("show_output.jl")

include("inputdata/inputdata1d.jl")
include("inputdata/inputdata2d.jl")
include("functions/function1d.jl")
include("functions/function2d.jl")

include("meshes/1d/mesh1d.jl")
include("meshes/2d/mesh2d.jl")

export generate_inputdata1d, generate_inputXi1d, generate_inputdata2d, generate_inputXi2d
export NNPP, NNPPh, NNRH, CNN, print_output, Network
export Traindata_Settings,Model_Settings, train_network
export train_network2dcnn,  generate_traindataset1d, generate_traindataset2d
export generate_validdataset1d, generate_validdataset2d

end


text:
include("CSVInput.jl")
include("GAMSInput.jl")
include("XLSXInput.jl")

export get_col
export get_descriptor
export get_name
export get_range
export get_sheet
export set_col!
export set_descriptor!
export set_name!
export set_range!
export set_sheet!

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
"""
    struct IntegralOperator{T,K,S,V} <: AbstractMatrix{T}

Representation of an integral operator which takes a density ϕ defined on `Y`
and integrates it with `kernel` for all elements `x ∈ X`.
"""
struct IntegralOperator{T,K,S,V} <: AbstractMatrix{T}
    kernel::K
    target_surface::S
    source_surface::V
end

kernel(iop::IntegralOperator) = iop.kernel
target_surface(iop::IntegralOperator) = iop.target_surface
source_surface(iop::IntegralOperator) = iop.source_surface

kernel_type(iop::IntegralOperator) = iop |> kernel |> kernel_type
combined_field_coefficients(iop::IntegralOperator) = iop |> kernel |> combined_field_coefficients

IntegralOperator{T}(k,X,Y) where {T} = IntegralOperator{T,typeof(k),typeof(X),typeof(Y)}(k, X, Y)

function IntegralOperator(k, X, Y=X)
    T = return_type(k)
    assert_concrete_type(T)
    IntegralOperator{T}(k, X, Y)
end

function Base.size(iop::IntegralOperator)
    X = target_surface(iop)
    Y = source_surface(iop)
    (length(dofs(X)), length(dofs(Y)))
end

# kernel_type(iop::IntegralOperator) = kernel_type(iop.kernel)
function Base.getindex(iop::IntegralOperator,i::Integer,j::Integer)
    k = kernel(iop)
    targets = target_surface(iop) |> dofs
    sources = source_surface(iop) |> dofs
    return k(targets[i],sources[j])*weight(sources[j])
end

# convenience constructors
SingleLayerOperator(op::AbstractPDE,X,Y=X)        = IntegralOperator(SingleLayerKernel(op), X, Y)
DoubleLayerOperator(op::AbstractPDE,X,Y=X)        = IntegralOperator(DoubleLayerKernel(op), X, Y)
AdjointDoubleLayerOperator(op::AbstractPDE,X,Y=X) = IntegralOperator(AdjointDoubleLayerKernel(op), X, Y)
HyperSingularOperator(op::AbstractPDE,X,Y=X)      = IntegralOperator(HyperSingularKernel(op), X, Y)

ambient_dimension(iop::IntegralOperator) = ambient_dimension(iop.kernel)

# Applying Laplace's double-layer to a constant will yield either 1 or -1,
# depending on whether the target point is inside or outside the obstacle
function isinside(x::SVector, mesh::NystromMesh)
    N   = ambient_dimension(mesh)
    pde = Laplace(dim=N)
    K   = DoubleLayerKernel(pde)
    u   = sum(dofs(mesh)) do source
         K(x, source) * weight(source)
    end
    u + 0.5 < 0
    # u < 0
end
isinside(x::Tuple,mesh::NystromMesh) = isinside(SVector(x), mesh)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
import .Hyperopt

@params struct HyperoptAlg <: AbstractOptimizer
    sub_alg::AbstractOptimizer
end

"""
    HyperoptOptions: options performing starting point optimization using Hyperopt.jl

- `sub_options`: options for the sub-optimizer.
- `lb`: Lower bound of starting point, if don't specify it, the default value will be `nothing`,
            then will end up be replaced by the lower bound of optimization problem.
- `ub`: Upper bound of starting point, same as above.
- `searchspace_size::Integer`: How many potential starting points we generate.
- `iters::Integer`: Among all generated potential starting points, how many of them will be evaluated.
- `sampler::Hyperopt.Sampler`: An instance of 'Hyperopt.Sampler', which decides search algorithm.
- `ctol`: infeasibility tolerance for accepting a solution as feasible
- `keep_all`: if true, all the solutions of the sub-problems will be saved
"""
@params struct HyperoptOptions
    sub_options
    lb
    ub
    searchspace_size::Integer
    iters::Integer
    sampler::Hyperopt.Sampler
    ctol
    keep_all::Bool
end
function HyperoptOptions(;
    sub_options, lb = nothing, ub = nothing,
    iters = 40, searchspace_size = iters,
    sampler = RandomSampler(), ctol = 1e-5,
    keep_all = true,
)
    return HyperoptOptions(
        sub_options, lb, ub, searchspace_size,
        iters, sampler, ctol, keep_all
    )
end

@params struct HyperoptWorkspace <: Workspace
    sub_workspace::Workspace
    x0::AbstractVector
    options::HyperoptOptions
end

function Workspace(model::VecModel, alg::HyperoptAlg, x0::AbstractVector; options, kwargs...)
    sub_options = options.sub_options isa Function ? options.sub_options(1) : options.sub_options
    return HyperoptWorkspace(
        Workspace(
            model, alg.sub_alg, copy(x0);
            options = sub_options, kwargs...,
        ),
        x0,
        options,
    )
end

"""
    HyperoptResult

When using multiple x0 in [`optimize`](@ref), return this result, including following fields:
- `minimum`: the optimal value among all searches.
- `minmizer`: the optimal solutions found.
- `results`: all the search results.
- `optimal_ind`: the index of the optimal solution in `results`.
"""
@params struct HyperoptResult <: AbstractResult
    minimum
    minimizer
    results
    optimal_ind
end

function optimize!(workspace::HyperoptWorkspace)
    @unpack options, x0, sub_workspace = workspace
    @unpack model = sub_workspace
    @unpack lb, ub, searchspace_size, ctol, keep_all = options
    @unpack sampler, iters = options

    lb = lb === nothing ? getmin(model) : lb
    ub = ub === nothing ? getmax(model) : ub
    @assert all(isfinite, lb) && all(isfinite, ub) throw("Please use finite bounds for the starting point search.")

    @info "Multistart algorithm started."
    @assert searchspace_size >= 1 "searchspace_size must be a positive integer. "

    # Generate a search space
    _sampler = sampler isa Hyperopt.Hyperband ? sampler.inner : sampler
    if _sampler isa Hyperopt.GPSampler
        params, candidates = get_linrange_candidates(lb, ub, searchspace_size)
        objective = (i, x0...) -> begin
            if sampler isa Hyperopt.Hyperband
                @assert options.sub_options isa Function
                sub_options = options.sub_options(i)
                _sub_workspace = @set sub_workspace.options = sub_options
            else
                _sub_workspace = sub_workspace
            end
            _x0 = [x0...]
            reset!(_sub_workspace, _x0)
            return optimize!(_sub_workspace), _x0
        end
    else
        params, candidates = get_sobol_candidates(lb, ub, x0, searchspace_size)
        objective = (i, x0) -> begin
            if sampler isa Hyperopt.Hyperband
                @assert options.sub_options isa Function
                sub_options = options.sub_options(i)
                _sub_workspace = @set sub_workspace.options = sub_options
            else
                _sub_workspace = sub_workspace
            end
            reset!(_sub_workspace, x0)
            return optimize!(_sub_workspace), x0
        end
    end
    _ho = Hyperopt.Hyperoptimizer(
        iterations = iters, sampler = sampler,
        objective = objective, params = params,
        candidates = candidates,
    )
    ho = @set _ho.objective = (args...) -> begin
        if length(args) == 2 && sampler isa Hyperopt.Hyperband
            if _sampler isa Hyperopt.GPSampler
                res, _ = _ho.objective(args[1], args[2]...)
            else
                res, _ = _ho.objective(args[1], args[2])
            end
        else
            res, _ = _ho.objective(args...)
        end
        push!(_ho.results, res.minimum)
        if sampler isa Hyperopt.Hyperband
            return isfeasible(model, res.minimizer, ctol = ctol) ? res.minimum : Inf, res.minimizer
        else
            return res
        end
    end
    if _sampler isa Union{Hyperopt.LHSampler, Hyperopt.CLHSampler, Hyperopt.GPSampler}
        Hyperopt.init!(_sampler, ho)
    end
    if sampler isa Hyperopt.Hyperband
        optimal = Hyperopt.hyperband(ho)
        return HyperoptResult(optimal[1], optimal[3], nothing, nothing)
    else
        hor = map(ho) do (args)
            ho.objective(values(args)...)
        end
        ind = argmin(map(res -> (isfeasible(model, res.minimizer, ctol = ctol) ? res.minimum : Inf), hor))
        optimal = hor[ind]
        if keep_all
            return HyperoptResult(optimal.minimum, optimal.minimizer, identity.(hor), ind)
        else
            return HyperoptResult(optimal.minimum, optimal.minimizer, nothing, nothing)
        end
    end
end

function get_linrange_candidates(lb, ub, n)
    return Tuple(Symbol.(:x0_, 1:length(lb))), Tuple(LinRange.(lb, ub, n))
end

function get_sobol_candidates(lb, ub, x0, n)
    sobol_seq_generator = SobolSeq(lb, ub)
    return (:x0,), ([[x0]; collect(next!(sobol_seq_generator) for _ = 1:n-1)],)
end

function RandomSampler(args...; kwargs...)
    return Hyperopt.RandomSampler(args...; kwargs...)
end
function Hyperband(args...; kwargs...)
    return Hyperopt.Hyperband(args...; kwargs...)
end
function LHSampler(args...; kwargs...)
    return Hyperopt.LHSampler(args...; kwargs...)
end
function CLHSampler(dims = Hyperopt.Continuous())
    return Hyperopt.CLHSampler(dims = [dims])
end
function GPSampler(args...; kwargs...)
    return Hyperopt.GPSampler(Hyperopt.Min, args...; kwargs...)
end


text:
using heuristicopt

function main()
    population = reshape(rand(1:10, 100), (10, 10))
    crossoverRate = 0.5
    mutationRate = 0.1
    iterations = 100
    gaparam = GA.GAParameter(
                         crossoverRate,
                         mutationRate,
                         iterations
                         )

    select = GATool.selectrandom
    crossover = GATool.crossoverOne
    mutate = GATool.mutaterandom
    gafunc = GA.GAFunction(select, crossover, mutate)
    model = GA.GAModel(
        gaparam,
        gafunc,
        objfunc
    )
    predictPopulation, bestFitness = GA.predict(model, population)
    print("bestFitness", bestFitness)
end

function objfunc(individual)
    0
end

main()

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
module CommonLH

using ArgCheck, DocStringExtensions, Formatting, Printf, StatsBase

# Check
export check_float, check_float_array, is_monotone
# Display
export show_text_file, show_string_vector
export MultiIO, print_flush, println_flush
# Keyword arguments
export KwArgs, default_value, has_default, kw_arg
# User input
export ask_for_choice, ask_yes_no
# Vector
export bisecting_indices, find_indices, find_index,
    all_at_least, all_at_most, all_greater, all_less,
    any_at_least, any_at_most, any_greater, any_less,
    any_nan
# Grids
export AbstractGrid, LinearGrid, LinSpacedGrid, PowerSpacedGrid, grid, intervals
# Probability matrices
export validate_prob_matrix, validate_prob_vector, prob_j, prob_k, prob_j_k, prob_k_j, ev_given_j, ev_given_k, scale_prob_array!
# Arrays
export bracket_array!, scale_array!
# Discretize
export discretize, bin_edges_from_percentiles, discretize_given_percentiles, count_bins

include("kwargs.jl")
include("check.jl")
include("display.jl")
include("user_input.jl")
include("vector.jl")
include("grids.jl")
include("probabilities.jl")
include("arrays.jl");
include("stats/discretize.jl");

# export validate, validate_scalar
# include("check.jl")

# include("matrix.jl")

# include("statistics.jl")

# export test_header, test_divider, test_dir
# include("testing.jl")

end # module


text:
### Observation model  ############################################################################

using LinearAlgebra

function observation_weight(pomdp::SingleOCFPOMDP, sp::SingleOCFState, o::SingleOCFObs)

    if ( is_observation_absent(pomdp, o) )
        o_absent = true
    else
        o_absent = false
    end

    if (is_state_absent(pomdp, sp) )
        sp_absent = true
    else
        sp_absent = false
    end

   if ( AutomotivePOMDPs.is_observable_fixed(pomdp.ego_vehicle.state, VehicleState(VecSE2(pomdp.ego_vehicle.state.posG.x+sp.ped_s, pomdp.ego_vehicle.state.posG.y+sp.ped_T, 0.), 0.0), pomdp.env) == false )
   # if ( AutomotiveSensors.occlusion_checker(VehicleState(VecSE2(0., 0., 0.), 0.0), VehicleState(VecSE2(sp.ped_s, sp.ped_T, 0.), 0.0), pomdp.env.obstacles) == false )
        occluded = true
    else
        occluded = false
    end

    if ( o_absent && sp_absent )   # absent
        return 1.0
    end

    if ( o_absent && occluded )     # occluded
        return 1.0
    end

    if ( !o_absent && !sp_absent && !o_absent )       # visible

        std_obs = zeros(Float64, 4, 4)
        std_obs[1,1] = 0.2        # x
        std_obs[2,2] = 0.2        # y
        std_obs[3,3] = 0.2        # theta
        std_obs[4,4] = 0.2        # velocity
        ob_dist = MvNormal([sp.ped_s, sp.ped_T, sp.ped_theta ,sp.ped_v], std_obs)

        (ego_y_state_space, ego_v_state_space) = getEgoDataInStateSpace(pomdp, o.ego_y, o.ego_v)

        obs_cont = SVector(ego_y_state_space, ego_v_state_space, o.ped_s, o.ped_T, o.ped_theta, o.ped_v)
        obs_int, obs_weight = interpolants(pomdp.state_space_grid, obs_cont)

        po = 0.0
        N = length(obs_int)
        for i=1:N
            obs_int_state = pomdp.state_space[obs_int[i]]
            po += obs_weight[i]*pdf(ob_dist, [obs_int_state.ped_s, obs_int_state.ped_T, obs_int_state.ped_theta, obs_int_state.ped_v])
        end
        po = po / N
        return po

    else
        return 0.0
    end

end

function POMDPs.observation(pomdp::SingleOCFPOMDP, a::SingleOCFAction, sp::SingleOCFState)

    states = SingleOCFObs[]
    sizehint!(states, 100);
    probs = Float64[]
    sizehint!(probs, 100);

    push!(states, sp)
    push!(probs, 1)

    return SparseCat(states,probs)

end


text:
# --- SYSTEM PACKAGES ---------------------------------------------------------- #
using MAT
using Dates
using Kegg
using ProgressMeter
# ------------------------------------------------------------------------------ #

# --- LOCAL -------------------------------------------------------------------- #
include("Extensions.jl")
include("Utility.jl")
include("Types.jl")
include("VFFParser.jl")
include("Cobra.jl")
include("QAQC.jl")
# ------------------------------------------------------------------------------ #

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
using Keldysh, Test

@testset "contour" begin
  let c = FullContour(tmax=2.0, β=5.0)
    @test nbranches(c) == 3

    b = map(x -> x.domain, c.branches)
    @test b == (forward_branch, backward_branch, imaginary_branch)

    c = twist(c)
    b = map(x -> x.domain, c.branches)
    @test b == (backward_branch, imaginary_branch, forward_branch)

    for b in (forward_branch, backward_branch, imaginary_branch)
      @test c[b].domain == b
    end

    @test forward_branch ∈ c
    @test backward_branch ∈ c
    @test imaginary_branch ∈ c
  end

  let c = KeldyshContour(tmax=2.0)
    @test nbranches(c) == 2

    b = map(x -> x.domain, c.branches)
    @test b == (forward_branch, backward_branch)

    c = twist(c)
    b = map(x -> x.domain, c.branches)
    @test b == (backward_branch, forward_branch)

    for b in (forward_branch, backward_branch)
      @test c[b].domain == b
    end

    @test forward_branch ∈ c
    @test backward_branch ∈ c
    @test imaginary_branch ∉ c
  end

  let c = ImaginaryContour(β=5.0)
    @test nbranches(c) == 1

    b = map(x -> x.domain, c.branches)
    @test b == (imaginary_branch,)

    for b in (imaginary_branch,)
      @test c[b].domain == b
    end

    @test forward_branch ∉ c
    @test backward_branch ∉ c
    @test imaginary_branch ∈ c
  end
end