text:
"""
    Module Solver
Information about finite element spaces and DG solver
"""
module FiniteElement

import Basis1D
import Basis2D
import Geometry
using DGApplication
using LinearAlgebra
using FELinearAlgebra

export Solver
export dg_rhs!
export nstates, nelements, dofsV, dofsF, nQV, nQF

"""
    struct Solver
All DG-specific spatial discretization/solver information
"""
struct Solver
    " Mesh on which FE space is defined "
    mesh::Geometry.Mesh
    " Application defining the PDE itself we are spatial discretizing "
    app::Application
    " Order of DG method "
    order::Int
    " Reference nodes in 1D. N-D version can be formed by kronecker product "
    x1D
    " Gauss quadrature points in 1D. N-D version can be formed by kronecker product "
    xQ1D
    " Gauss quadrature weights in ND, for use within volume "
    wQV
    " Gauss quadrature weights in (N-1)D, for use along faces "
    wQF
    ###########
    # TODO: should these matrices possibly be stored in a more abstract fashion?
    ###########
    " Global element mass matrix, defined as ∫ Φ_i Φ_j dx "
    Mel
    " Global element gradient matrix, defined as ∫ ∇_lΦ_i Φ_j dx "
    Gels
    " Global element stiffness matrix, defined as ∫ ∇Φ_i ⋅ ∇Φ_j dx"
    Sel
    ###########
    # End TODO
    ###########
    " Interpolation matrix from element nodes to element quadrature points "
    InterpV
    " Interpolation matrix from face nodes to face quadrature points "
    InterpF
    """
    Gradient of basis functions interpolating from element nodes to
    element quadrature points. Can alternatively be viewed as local element
    gradient matrix K, defined as ∇_lΦ_i
    """
    dInterpV
    """
    Gradient of basis functions interpolating from face nodes to face
    quadrature points. Can alternatively be viewed as a local element face
    gradient matrix, defined as ∇_lΦ_i along faces
    """
    dInterpF
    " Interpolation matrix from bilinear map Tk points to element nodes "
    InterpTk
    " Interpolation matrix from bilinear map Tk points to element quadrature points "
    InterpTkQ

    #### Work arrays
    " Solution vector stored on volume quadrature points "
    u_interpV
    " Solution vector stored on face quadrature points "
    u_interpF
    " Array of `SolutionVector`, storing ∇u in each dimension "
    ∇u
    " Single `SolutionVector`, storing ∇u with fastest dimension = dim*nodes"
    ∇u_compute
    " Gradient ∇u stored on volume quadrature points "
    ∇u_interpV
    " Gradient ∇u stored on face quadrature points "
    ∇u_interpF

    function Solver(order::Int, mesh::Geometry.Mesh, app::Application)
        # Initialize node, quadrature, and interpolation matrix information
        (x1D,xQ1D,wQV,wQF,xTk1D,InterpV,InterpF,
            dInterpV,dInterpF,InterpTk,InterpTkQ) =
            precompute_interpolation_matrices(order)
        # Update Mesh to know about solver information
        Geometry.setup_nodes!(mesh,InterpTk,order,compute_global=true)
        Geometry.setup_quads!(mesh,InterpTkQ)
        Geometry.precompute_jacobians!(mesh,xTk1D,xQ1D)
        if is_second_order(app)
            Geometry.precompute_ldg_switch!(mesh)
        end
        # Initialize work arrays
        ns = DGApplication.nstates(app)
        ne = mesh.n_elements
        #u = SolutionVector{Float64}(size(InterpV,2), ns, ne)
        u_interpV = Array{Float64,3}(undef, size(InterpV,1),ns,ne)
        u_interpF = Array{Float64,4}(undef, size(InterpF,1),ns,Geometry.N_FACES,ne)
        if is_second_order(app)
            ∇u = [SolutionVector{Float64}(size(InterpV,2),ns,ne) for i=1:Geometry.DIM]
            ∇u_compute = SolutionVector{Float64}(Geometry.DIM*size(InterpV,2),ns,ne)
            ∇u_interpV = Array{Float64,4}(undef, size(InterpV,1),ns,ne,Geometry.DIM)
            ∇u_interpF = Array{Float64,5}(undef, size(InterpF,1),ns,Geometry.N_FACES,ne,Geometry.DIM)
        else
            ∇u = []
            ∇u_compute = []
            ∇u_interpV = []
            ∇u_interpF = []
        end
        # Initialize local matrices
        Locals = precompute_local_matrices(order,mesh,app,wQV,InterpV,InterpF,dInterpV,dInterpF)

        return new(mesh,app,order, x1D,xQ1D,wQV,wQF, Locals...,
            InterpV,InterpF,dInterpV,dInterpF,InterpTk,InterpTkQ,
            u_interpV,u_interpF, ∇u,∇u_compute,∇u_interpV,∇u_interpF)
    end

end

" Number of states in state vector "
DGApplication.nstates(solver::Solver) = DGApplication.nstates(solver.app)
" Number of elements in mesh "
nelements(solver::Solver) = solver.mesh.n_elements
" Degrees of freedom per element "
dofsV(solver::Solver) = size(solver.InterpV,2)
" Degrees of freedom per element face "
dofsF(solver::Solver) = size(solver.InterpF,2)
" Number of quadrature points per element "
nQV(solver::Solver) = size(solver.InterpV,1)
" Number of quadrature points per element face "
nQF(solver::Solver) = size(solver.InterpF,1)

"""
    precompute_interpolation_matrices(order)
Precompute all interpolation matrices used by DG method.
"""
function precompute_interpolation_matrices(order)
    # Reference nodes
    x1D = Basis1D.chebyshev(order)
    # Gauss quadrature points
    (xQ1D,wQ1D) = Basis1D.gauss_quad(2*order)
    wQV = kron(wQ1D,wQ1D)
    wQF = deepcopy(wQ1D)

    # Volume/Element interpolation
    InterpV = Basis2D.interpolation_matrix2D((x1D,x1D),(xQ1D,xQ1D))
    # Face interpolation
    InterpF = Basis1D.interpolation_matrix1D(x1D, xQ1D)

    # Local gradient matrix = gradient of reference bases on the quad points
    dInterpV = Basis2D.dPhi2D((x1D,x1D),(xQ1D,xQ1D))
    # Local face gradient matrix
    dInterpF = Basis1D.dPhi1D(x1D, xQ1D)

    # Interpolation from points on which bilinear_mapping is defined
    xTk1D = Basis1D.chebyshev(1)
    InterpTk = Basis2D.interpolation_matrix2D((xTk1D,xTk1D),(x1D,x1D))
    InterpTkQ = Basis2D.interpolation_matrix2D((xTk1D,xTk1D),(xQ1D,xQ1D))

    return (x1D,xQ1D,wQV,wQF,xTk1D,InterpV,InterpF,dInterpV,dInterpF,InterpTk,InterpTkQ)
end

"""
    precompute_local_matrices()
Precompute all element-local matrices used by DG method.
"""
function precompute_local_matrices(order, mesh::Geometry.Mesh, app::Application,
    wQV,InterpV,InterpF,dInterpV,dInterpF)
    dofsV = size(InterpV,2)
    dofsF = size(InterpF,2)
    nQV = size(InterpV,1)
    nQF = size(InterpF,1)

    # Mass matrix = ∫ Φ_i Φ_j dx
    Mdata = Array{Float64,4}(undef, dofsV,dofsV,1,mesh.n_elements)
    for iK = 1:mesh.n_elements
        # TODO: Eventually this should be done in matrix-free approach
        scaling = wQV.*vec(mesh.Jkdet[:,iK])
        Mdata[:,:,1,iK] = InterpV'*Diagonal(scaling)*InterpV
    end
    Mel = LocalMatrix(Mdata)

    # Gradient matrix = ∫ ∇ Φ_i Φ_j dx
    if is_second_order(app)
        Gdata = zeros(Float64, Geometry.DIM,dofsV,dofsV,1,mesh.n_elements)
        for iK = 1:mesh.n_elements
            # TODO: Eventually this should be done in matrix-free approach
            Jks = @view mesh.Jk[:,:,:,iK]
            # op_q = wQ * J^{-1}*det(J) == wQ * adj(J)
            op_q = Array{Float64,3}(undef, size(Jks))
            op_q[1,1,:] = wQV.* Jks[2,2,:]
            op_q[1,2,:] = wQV.*-Jks[1,2,:]
            op_q[2,1,:] = wQV.*-Jks[2,1,:]
            op_q[2,2,:] = wQV.* Jks[1,1,:]
            @views for iQ = 1:nQV
                grad = op_q[:,:,iQ]*(dInterpV[iQ,:,:]')
                Gdata[1,:,:,1,iK] .+= InterpV[iQ,:]'*grad[1,:]
                Gdata[2,:,:,1,iK] .+= InterpV[iQ,:]'*grad[2,:]
            end
        end
        Gdata = reshape(Gdata,Geometry.DIM*dofsV,dofsV,1,mesh.n_elements)
    else
        Gdata = zeros(Float64, 0,0,0,0)
    end
    Gels = LocalMatrix(Gdata)

    # Stiffness matrix = ∫ ∇Φ_i ⋅ ∇Φ_j dx
    Sdata = zeros(Float64, dofsV,dofsV,1,mesh.n_elements)
    for iK = 1:mesh.n_elements
        # TODO: Eventually this should be done in matrix-free approach
        Jks = @view mesh.Jk[:,:,:,iK]
        # op_q = wQ * J^{-1}*det(J)*J^{-T} == wQ/det(J)*adj(J)*adj(J)^T
        op_q = Array{Float64,3}(undef, size(Jks))
        op_q[1,1,:] =  wQV./mesh.Jkdet[:,iK].*(Jks[1,2,:].^2+Jks[2,2,:].^2)
        op_q[1,2,:] = -wQV./mesh.Jkdet[:,iK].*(Jks[1,2,:].*Jks[1,1,:]+Jks[2,2,:].*Jks[2,1,:])
        op_q[2,1,:] = op_q[1,2,:]
        op_q[2,2,:] =  wQV./mesh.Jkdet[:,iK].*(Jks[1,1,:].^2+Jks[2,1,:].^2)
        @views for iQ = 1:nQV
            grad = op_q[:,:,iQ]*(dInterpV[iQ,:,:]')
            Sdata[:,:,1,iK] .+= dInterpV[iQ,:,:]*grad
        end
    end
    Sel = LocalMatrix(Sdata)

    return (Mel, Gels, Sel)
end

"""
    dg_rhs!(du, u_curr, all_param, t)
Given u_curr, evaluate the RHS of the system of ODE: u_t = f(u_curr, t)
where f represents the spatial discretization of our application
Note that all_param is a tuple containing (solver,) TODO: what else?
this is equivalent to du = Mel \\ (S+K_c(u)-F_c(u)) in the case of 1st order apps or
∇u = Mel \\ (-Gels*u + F_l(u)),
du = Mel \\ ( S+K_c(u)+K_d(u,∇u)-F_c(u)-F_d(u,∇u) ) in the case of 2nd order apps.
"""
function dg_rhs!(du, u_curr, all_param, t)
    (solver,) = all_param # TODO: TBD - what else?

    interpolate!(solver, u_curr, solver.u_interpV, solver.u_interpF)
    if is_second_order(solver.app)
        # ∇u = -Gs*u
        mul!(solver.∇u_compute,solver.Gels,u_curr,-1.0,0.0)
        convert_to_dim_array!(solver.∇u, solver.∇u_compute)
        # ∇u += F_l(u)
        local_flux!(solver)

        # ∇u = (I⊗M)\∇u
        for l=1:Geometry.DIM
            ldiv!(solver.Mel, solver.∇u[l])
            @views interpolate!(solver, solver.∇u[l],
                solver.∇u_interpV[:,:,:,l], solver.∇u_interpF[:,:,:,:,l])
        end

    end
    # du = s
    source_term!(du, solver, t)
    # du -= F_c(u)
    convect_flux!(du, solver)
    # du += K_c(u)
    convect_volume!(du, solver)
    if is_second_order(solver.app)
        # du -= F_d(u,∇u)
        diffusive_flux!(du, solver)
        # du += K_d(u,∇u)
        diffusive_volume!(du, solver)
    end

    # du = M\du
    ldiv!(solver.Mel, du)
    # Mass matrix might be addable directly to ODE - double check documentation

end

"""
    interpolate!(solver::Solver, curr::SolutionVector, to_interpV, to_interpF)
Interpolate curr from SolutionVector to volume quadrature points and store in
to_interpV. Interpolate curr from SolutionVector to face quadrature points and
store in to_interpF.
"""
function interpolate!(solver::Solver, curr::SolutionVector, to_interpV, to_interpF)
    nstates = true_size(curr, 2)
    # First grab solution on faces and pack into array on_faces
    on_faces = Array{Float64,3}(undef, dofsF(solver),nstates,Geometry.N_FACES)
    for iK = 1:nelements(solver)
        for iF = 1:Geometry.N_FACES
            for iS = 1:nstates
                for iFN = 1:dofsF(solver)
                    on_faces[iFN,iS,iF] = curr.data[solver.mesh.ef2n[iFN,iF],iS,iK]
                end
            end
        end
        # Face interpolation: to_interpF = InterpF*on_faces
        on_face_vecs = reshape(on_faces, dofsF(solver), nstates*Geometry.N_FACES)
        @views to_interpF_vecs = reshape(to_interpF[:,:,:,iK], nQF(solver), nstates*Geometry.N_FACES)
        mul!(to_interpF_vecs, solver.InterpF,on_face_vecs)

        # Volume interpolation: to_interpV = InterpV*curr
        curr_vecs = @view curr.data[:,:,iK]
        to_interpV_vecs = @view to_interpV[:,:,iK]
        mul!(to_interpV_vecs, solver.InterpV,curr_vecs)
    end
end

"""
    convert_to_dim_array!(∇u, ∇u_compute)
Update ∇u (an array of `SolutionVector`s, one for each direction) from the
computed values stored in the single, long `SolutionVector` ∇u_compute
"""
function convert_to_dim_array!(∇u, ∇u_compute)
    (n, ns, ne) = true_size(∇u[1])
    dim = length(∇u)
    (n_times_dim, nsc, nec) = true_size(∇u_compute)
    @assert n*dim == n_times_dim && ns == nsc && ne == nec
    for iK = 1:ne
        for iS = 1:ns
            for iN = 1:n
                for l = 1:dim
                    ∇u[l].data[iN,iS,iK] = ∇u_compute.data[(iN-1)*dim+l,iS,iK]
                end
            end
        end
    end
end

"""
    local_flux!(solver::Solver)
Compute solver.∇u += F_l(u), the local DG flux for use in formulation of
second-order terms. Uses a downwind formulation for flux term.
"""
function local_flux!(solver::Solver)
    uK = zeros(nstates(solver))
    uN = zeros(nstates(solver))
    fstar = Array{Float64,2}(undef, nstates(solver),nQF(solver))
    JWI = Array{Float64,2}(undef, nQF(solver),dofsF(solver))
    flux_contribution = Array{Float64,2}(undef, nstates(solver),dofsF(solver))

    for iK = 1:nelements(solver)
        # There are l equations to handle for local flux term
        for l = 1:Geometry.DIM
            for iF = 1:Geometry.N_FACES
                nK = solver.mesh.e2e[iF,iK]
                normalK = solver.mesh.normals[l,iF,iK]

                if nK < 0 # Boundary condition
                    for iFQ = 1:nQF(solver)
                        uK .= solver.u_interpF[iFQ,:,iF,iK]
                        @views boundary_flux_l!(fstar[:,iFQ], solver.app, uK, nK, normalK)
                    end
                else # Inter-element flux_l
                    nF = solver.mesh.e2f[iF,iK]
                    for iFQ = 1:nQF(solver)
                        uK .= solver.u_interpF[iFQ                    ,:,iF,iK]
                        uN .= solver.u_interpF[solver.mesh.nfq2fq[iFQ],:,nF,nK]
                        @views numerical_flux_l!(fstar[:,iFQ], solver.app, uK, uN, solver.mesh.ldg_switch[iF,iK], normalK)
                    end
                end
                # JWI = Diag(J)*Diag(w)*InterpF
                mul!(JWI, Diagonal(solver.wQF .* solver.mesh.JkFdet[:,iF,iK]), solver.InterpF)
                # flux contribution = fstar*JWI
                mul!(flux_contribution, fstar, JWI)
                # Add up local face contributions into global ∇u array
                for iFN = 1:dofsF(solver)
                    for iS = 1:nstates(solver)
                        solver.∇u[l].data[solver.mesh.ef2n[iFN,iF],iS,iK] += flux_contribution[iS,iFN]
                    end
                end
            end
        end
    end
end

"""
    source_term!(du, solver::Solver, t)
Compute du = source(x,t) within element volumes
"""
function source_term!(du, solver::Solver, t)
    JWI = Array{Float64,2}(undef, nQV(solver),dofsV(solver))
    f = Array{Float64,2}(undef, nstates(solver),nQV(solver))
    for iK = 1:nelements(solver)
        # JWI = Diag(J)*Diag(w)*InterpV
        mul!(JWI, Diagonal(solver.wQV .* solver.mesh.Jkdet[:,iK]), solver.InterpV)

        for iQ = 1:nQV(solver)
            @views source!(f[:,iQ], solver.app, t)
        end
        # du[:,:,iK] = JWI'*f
        @views mul!(du.data[:,:,iK], JWI', f')
    end
end

"""
    convect_flux!(du, solver::Solver)
Compute du -= Fc(u), which is the 1st order (convection) flux term along faces
"""
function convect_flux!(du, solver::Solver)
    uK = zeros(nstates(solver))
    uN = zeros(nstates(solver))
    fstar = Array{Float64,2}(undef, nstates(solver),nQF(solver))
    JWI = Array{Float64,2}(undef, nQF(solver),dofsF(solver))
    flux_contribution = Array{Float64,2}(undef, nstates(solver),dofsF(solver))

    for iK = 1:nelements(solver)
        for iF = 1:Geometry.N_FACES
            nK = solver.mesh.e2e[iF,iK]
            normalK = @view solver.mesh.normals[:,iF,iK]

            if nK < 0 # Boundary condition
                for iFQ = 1:nQF(solver)
                    uK .= solver.u_interpF[iFQ,:,iF,iK]
                    @views boundary_flux_c!(fstar[:,iFQ], solver.app, uK, nK, normalK)
                end
            else # Inter-element flux_c
                nF = solver.mesh.e2f[iF,iK]
                for iFQ = 1:nQF(solver)
                    uK .= solver.u_interpF[iFQ                    ,:,iF,iK]
                    uN .= solver.u_interpF[solver.mesh.nfq2fq[iFQ],:,nF,nK]
                    @views numerical_flux_c!(fstar[:,iFQ], solver.app, uK, uN, normalK)
                end
            end

            # JWI = Diag(J)*Diag(w)*InterpF
            mul!(JWI, Diagonal(solver.wQF .* solver.mesh.JkFdet[:,iF,iK]), solver.InterpF)
            # flux contribution = fstar*JWI
            mul!(flux_contribution, fstar, JWI)
            # Add up local face contributions into global du array
            for iFN = 1:dofsF(solver)
                for iS = 1:nstates(solver)
                    du.data[solver.mesh.ef2n[iFN,iF],iS,iK] -= flux_contribution[iS,iFN]
                end
            end
        end
    end
end

"""
    convect_volume!(du, solver::Solver)
Compute du += Kc(u), which is the 1st order (convection) local volume integral
within elements
"""
function convect_volume!(du, solver::Solver)
    uK = zeros(nstates(solver))
    fc = Array{Float64,3}(undef, nstates(solver), Geometry.DIM, nQV(solver))
    grad = Array{Float64,3}(undef, Geometry.DIM, nstates(solver), nQV(solver))
    # If not using Julia *, should initialize with DIM in slowest dimension

    for iK = 1:nelements(solver)
        for iQ = 1:nQV(solver)
            uK .= solver.u_interpV[iQ,:,iK]
            @views flux_c!(fc[:,:,iQ], solver.app, uK)
        end

        # TODO: Eventually this should be done in matrix-free approach
        Jks = @view solver.mesh.Jk[:,:,:,iK]
        # op_q = wQ * J^{-1}*det(J) == wQ * adj(J)
        op_q = similar(Jks)
        op_q[1,1,:] = solver.wQV.* Jks[2,2,:]
        op_q[1,2,:] = solver.wQV.*-Jks[1,2,:]
        op_q[2,1,:] = solver.wQV.*-Jks[2,1,:]
        op_q[2,2,:] = solver.wQV.* Jks[1,1,:]
        @views for iQ = 1:nQV(solver)
            # grad = op * fc'
            mul!(grad[:,:,iQ], op_q[:,:,iQ], fc[:,:,iQ]')
        end
        # du += dInterpV'*grad'
        @views mul!(du.data[:,:,iK], solver.dInterpV[:,:,1]', grad[1,:,:]', 1.0, 1.0)
        @views mul!(du.data[:,:,iK], solver.dInterpV[:,:,2]', grad[2,:,:]', 1.0, 1.0)
    end
end

"""
    function diffusive_flux!(du, solver)
Compute du += Fd(u,∇u), which is the 2nd order (diffusive) flux term along faces
"""
function diffusive_flux!(du, solver)
    uK = zeros(nstates(solver))
    uN = zeros(nstates(solver))
    DuK = zeros(nstates(solver), Geometry.DIM)
    DuN = zeros(nstates(solver), Geometry.DIM)
    fstar = Array{Float64,2}(undef, nstates(solver),nQF(solver))
    JWI = Array{Float64,2}(undef, nQF(solver),dofsF(solver))
    flux_contribution = Array{Float64,2}(undef, nstates(solver),dofsF(solver))

    for iK = 1:nelements(solver)
        for iF = 1:Geometry.N_FACES
            nK = solver.mesh.e2e[iF,iK]
            normalK = @view solver.mesh.normals[:,iF,iK]

            if nK < 0 # Boundary condition
                for iFQ = 1:nQF(solver)
                    uK .= solver.u_interpF[iFQ,:,iF,iK]
                    DuK .= solver.∇u_interpF[iFQ,:,iF,iK,:]
                    @views boundary_flux_d!(fstar[:,iFQ], solver.app, uK, DuK, nK, normalK)
                end
            else # Inter-element flux_c
                nF = solver.mesh.e2f[iF,iK]
                for iFQ = 1:nQF(solver)
                    uK .= solver.u_interpF[iFQ                    ,:,iF,iK]
                    uN .= solver.u_interpF[solver.mesh.nfq2fq[iFQ],:,nF,nK]
                    DuK .= solver.∇u_interpF[iFQ                    ,:,iF,iK,:]
                    DuN .= solver.∇u_interpF[solver.mesh.nfq2fq[iFQ],:,nF,nK,:]
                    @views numerical_flux_d!(fstar[:,iFQ], solver.app, uK, uN, DuK, DuN, solver.mesh.ldg_switch[iF,iK],normalK)
                end
            end

            # JWI = Diag(J)*Diag(w)*InterpF
            mul!(JWI, Diagonal(solver.wQF .* solver.mesh.JkFdet[:,iF,iK]), solver.InterpF)
            # flux contribution = fstar*JWI
            mul!(flux_contribution, fstar, JWI)
            # Add up local face contributions into global du array
            for iFN = 1:dofsF(solver)
                for iS = 1:nstates(solver)
                    du.data[solver.mesh.ef2n[iFN,iF],iS,iK] += flux_contribution[iS,iFN]
                end
            end
        end
    end
end

"""
    diffusive_volume!(du, solver::Solver)
Compute du -= Kd(u,∇u), which is the 2nd order (diffusive) local volume integral
within elements
"""
function diffusive_volume!(du, solver::Solver)
    uK = zeros(nstates(solver))
    DuK = zeros(nstates(solver),Geometry.DIM)
    fd = Array{Float64,3}(undef, nstates(solver), Geometry.DIM, nQV(solver))
    grad = Array{Float64,3}(undef, Geometry.DIM, nstates(solver), nQV(solver))
    # If not using Julia *, should initialize with DIM in slowest dimension

    for iK = 1:nelements(solver)
        for iQ = 1:nQV(solver)
            uK .= solver.u_interpV[iQ,:,iK]
            DuK .= solver.∇u_interpV[iQ,:,iK,:]
            @views flux_d!(fd[:,:,iQ], solver.app, uK, DuK)
        end

        # TODO: Eventually this should be done in matrix-free approach
        Jks = @view solver.mesh.Jk[:,:,:,iK]
        # op_q = wQ * J^{-1}*det(J) == wQ * adj(J)
        op_q = similar(Jks)
        op_q[1,1,:] = solver.wQV.* Jks[2,2,:]
        op_q[1,2,:] = solver.wQV.*-Jks[1,2,:]
        op_q[2,1,:] = solver.wQV.*-Jks[2,1,:]
        op_q[2,2,:] = solver.wQV.* Jks[1,1,:]
        @views for iQ = 1:nQV(solver)
            # grad = op * fd'
            mul!(grad[:,:,iQ], op_q[:,:,iQ], fd[:,:,iQ]')
        end
        # du -= dInterpV'*grad'
        @views mul!(du.data[:,:,iK], solver.dInterpV[:,:,1]', grad[1,:,:]', -1.0, 1.0)
        @views mul!(du.data[:,:,iK], solver.dInterpV[:,:,2]', grad[2,:,:]', -1.0, 1.0)
    end
end

include("FEImplicit.jl")

end # module FiniteElement

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
@defcomp Discontinuity begin

  region = Index()
  y_year=Parameter(index=[time], unit="year")
  y_year_0 = Parameter(unit="year")

  rand_discontinuity = Parameter(unit="unitless", default=.5)

  irefeqdis_eqdiscimpact=Variable(index=[region], unit="%")
  wincf_weightsfactor=Parameter(index=[region], unit="unitless", default=[1, 0.8, 0.8, 0.4, 0.8, 0.8, 0.6, 0.6])
  wdis_gdplostdisc=Parameter(unit="%", default=15.)

  igdpeqdis_eqdiscimpact=Variable(index=[time,region], unit="%")
  rgdp_per_cap_NonMarketRemainGDP=Parameter(index=[time,region], unit="\$/person")
  GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
  ipow_incomeexponent=Parameter(unit="unitless", default=-0.13333333333333333)

  igdp_realizeddiscimpact=Variable(index=[time,region], unit="%")
  occurdis_occurrencedummy=Variable(index=[time], unit="unitless")
  expfdis_discdecay=Variable(index=[time], unit="unitless")

  distau_discontinuityexponent=Parameter(unit="unitless", default=90.)

  idis_lossfromdisc=Variable(index=[time], unit="degreeC")
  tdis_tolerabilitydisc=Parameter(unit="degreeC", default=3.)
  rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
  pdis_probability=Parameter(unit="%/degreeC", default=20.)

  isatg_saturationmodification = Parameter(unit="unitless")
  isat_satdiscimpact=Variable(index=[time,region], unit="%")

  isat_per_cap_DiscImpactperCapinclSaturation=Variable(index=[time,region], unit="\$/person")
  rcons_per_cap_DiscRemainConsumption=Variable(index=[time, region], unit = "\$/person")
  rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit = "\$/person")

    function run_timestep(p, v, d, t)

        v.idis_lossfromdisc[t] = max(0, p.rt_g_globaltemperature[t] - p.tdis_tolerabilitydisc)

        if is_first(t)
            if v.idis_lossfromdisc[t]*(p.pdis_probability/100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t]=exp(-(p.y_year[t] - p.y_year_0)/p.distau_discontinuityexponent)
        else
            if v.idis_lossfromdisc[t]*(p.pdis_probability/100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            elseif v.occurdis_occurrencedummy[t-1] == 1
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t]=exp(-(p.y_year[t] - p.y_year[t-1])/p.distau_discontinuityexponent)
        end

        for r in d.region
            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor[r]*p.wdis_gdplostdisc

            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

            if is_first(t)
                v.igdp_realizeddiscimpact[t,r]=v.occurdis_occurrencedummy[t]*(1-v.expfdis_discdecay[t])*v.igdpeqdis_eqdiscimpact[t,r]
            else
                v.igdp_realizeddiscimpact[t,r]=v.igdp_realizeddiscimpact[t-1,r]+v.occurdis_occurrencedummy[t]*(1-v.expfdis_discdecay[t])*(v.igdpeqdis_eqdiscimpact[t,r]-v.igdp_realizeddiscimpact[t-1,r])
            end

            if v.igdp_realizeddiscimpact[t,r] < p.isatg_saturationmodification
                v.isat_satdiscimpact[t,r] = v.igdp_realizeddiscimpact[t,r]
            else
                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100-p.isatg_saturationmodification)*((v.igdp_realizeddiscimpact[t,r]-p.isatg_saturationmodification)/((100-p.isatg_saturationmodification)+(v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
            end
            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r]/100)*p.rgdp_per_cap_NonMarketRemainGDP[t,r]
            v.rcons_per_cap_DiscRemainConsumption[t,r] = p.rcons_per_cap_NonMarketRemainConsumption[t,r] - v.isat_per_cap_DiscImpactperCapinclSaturation[t,r]
        end
    end
end

ext: jl


text:
######################################################################################################
# Functions for implementing our sampling rules as well as generic way to compute f and ∇f           #
######################################################################################################
using JuMP;
import Tulip;

function hμ_in_lambda(hμ, hi, K)
    for i=1:K
        if ((i!=hi) && (hμ[hi]-hμ[i])<=eps())
            return false;
        end
    end
    return true;
end
function hμ_in_lambda_threshold(reward, K, threshold)
    for i=1:K
        if abs(reward[i]-threshold)<=eps()
            return false;
        end
    end
    return true;
end
function is_complete_square(n)
    p = floor(Int, sqrt(n));
    return p*p == n;
end
function pseudo_lipschitz(hμ, arms, L)
    K = length(arms);
    maxR = maximum([abs(hμ[i]-hμ[j])/(L*abs(arms[i]-arms[j])) for i=1:K for j=(i+1):K]);
    return maximum([maxR*L, L]);
end
function is_lipschitz(hμ, arms, L)
    R = pseudo_lipschitz(hμ, arms, L);
    return (R <= L+eps());
end

"""
Required by our method: solving the LP formulation of mixed-strategy 2-player zero-sum game
"""
function solveZeroSumGame(M_payoff, K, n_row)
    m = Model(with_optimizer(Tulip.Optimizer));
    @variable(m, x[1:K] >= 0)
    @variable(m, w)
    for j in 1:n_row
        @constraint(m, sum(M_payoff[j][k]*x[k] for k=1:K) >= w)
    end
    @constraint(m, sum(x[i] for i=1:K) == 1)
    @objective(m, Max, w)
    optimize!(m);
    f_success = termination_status(m);
    z = JuMP.value.(x);
    return z;
end

"""
Standard BAI: computing f and ∇f by our Proposition 1
"""
function compute_f_∇f_standard_bai(hw, hμ, ξ, hi, r, K)
    μbar = [[(hw[i]*ξ[i]+hw[j]*ξ[j])/(hw[i]+hw[j]) for j=1:K] for i=1:K];
    suboptimal = [i for i=1:K if i!=hi];
    # construct ∇f
    ∇f = [[0.0 for j=1:K] for i=1:K];
    for j in suboptimal
        ∇f[j][hi] = d(getexpfam(pep,hi),hμ[hi],μbar[hi][j]);
        ∇f[j][j] = d(getexpfam(pep,j),hμ[j],μbar[hi][j]);
    end
    # construct f
    f = [hw'∇f[j] for j in suboptimal];
    fmin = minimum(f);
    if r > eps()
        fidx = [j for (idxj,j) in enumerate(suboptimal) if (f[idxj]<fmin+r)]
    elseif abs(r)<eps()
        fidx = [suboptimal[argmin(f)]];
    else
        fidx = suboptimal;
    end
    return f, ∇f, fidx;
end

"""
Linear BAI: computing f and ∇f by our Proposition 1
"""
# Envelope theorem simplifies the computation of the confusing parameters
function alt_min_linear_bai(hw, hμ, arms, Vxinv)
    K = length(arms); dim = length(hμ);
    hr = [hμ'arms[k] for k=1:K];
    hi = argmax(hr);
    suboptimal = [i for i=1:K if i!=hi];
    # construct ∇f
    λ = zeros(dim,K);
    for k in suboptimal
        direction = arms[hi]-arms[k];
        λ[:,k] = hμ - (direction'hμ / ((direction')*Vxinv*direction)) * Vxinv*(direction);
    end
    ∇f = [[0.0 for i=1:K] for j=1:K];
    for k in suboptimal
        for i=1:K
            ∇f[k][i] = ((arms[i]')*(hμ-λ[:,k]))^2 / 2;
        end
    end
    # construct f
    f = [hw'∇f[j] for j in suboptimal];
    return minimum(f), f, ∇f, suboptimal;
end
function compute_f_∇f_linear_bai(hw, hμ, r, arms, Vxinv)
    fmin, f, ∇f, suboptimal = alt_min_linear_bai(hw, hμ, arms, Vxinv);
    if r > eps()
        fidx = [j for (idxj,j) in enumerate(suboptimal) if (f[idxj]<fmin+r)]
    elseif abs(r)<eps()
        fidx = [suboptimal[argmin(f)]];
    else
        fidx = suboptimal;
    end
    return f, ∇f, fidx;
end

"""
Linear Threshold: computing f and ∇f by our Proposition 1
"""
# Envelope theorem simplifies the computation of the confusing parameters
function alt_min_linear_threshold(hw, hμ, arms, Vxinv)
    K = length(pep.arms); dim = length(hμ);
    # construct ∇f
    λ = zeros(dim,K);
    for k=1:K
        x = pep.threshold - (hμ')*pep.arms[k]; ak = pep.arms[k];
        λ[:,k] = hμ + sign(x) * (x/((ak')*Vxinv*ak)) * Vxinv * ak;
    end
    ∇f = [[0.0 for i=1:K] for j=1:K];
    for k=1:K
        for i=1:K
            ∇f[k][i] = ((pep.arms[i]')*(hμ-λ[:,k]))^2 / 2;
        end
    end
    # construct f
    f = [hw'∇f[j] for j=1:K];
    return minimum(f), f, ∇f;
end
function compute_f_∇f_linear_threshold(pep, hw, hμ, r, Vxinv)
    fmin, f, ∇f = alt_min_linear_threshold(hw, hμ, pep.arms, Vxinv);
    if r > eps()
        fidx = [j for j=1:length(pep.arms) if (f[j]<fmin+r)]
    else
        fidx = [argmin(f)];
    end
    return f, ∇f, fidx;
end

"""
Lipschitz BAI: computing f and ∇f by our Proposition 1
"""
function confusing_parameter_function(θ, j, hw, hμ, hi, arms, L)
    K = length(arms);
    return sum([2*hw[k]*maximum([θ-L*abs(arms[k]-arms[j])-hμ[k], 0]) for k=1:K]) - sum([2*hw[k]*maximum([hμ[k]-θ-L*abs(arms[k]-arms[hi]), 0]) for k=1:K]);
end
# Envelope theorem simplifies the computation of the confusing parameters
function alt_min_lipschitz(hw, hμ, arms, L)
    K = length(arms); hi = argmax(hμ);
    suboptimal = [i for i=1:K if i!=hi];
    # confusing parameters
    Θ = zeros(K);
    for j in suboptimal
        Θ[j] = binary_search(θ -> confusing_parameter_function(θ, j, hw, hμ, hi, arms, L), hμ[j], hμ[hi]);
    end
    # construct ∇f
    λ = [[0.0 for i=1:K] for j=1:K];
    for j in suboptimal
        for k=1:K
            λ[j][k] = minimum([maximum([Θ[j]-L*abs(arms[k]-arms[j]), hμ[k]]), Θ[j]+L*abs(arms[k]-arms[hi])]);
        end
    end
    ∇f = [[0.0 for i=1:K] for j=1:K];
    for j in suboptimal
        for k=1:K
            ∇f[j][k] = (hμ[k]-λ[j][k])^2/2.0;
        end
    end
    # construct f
    f = [hw'∇f[j] for j in suboptimal];
    return minimum(f), f, ∇f, suboptimal;
end
function compute_f_∇f_lipschitz_bai(hw, hμ, r, arms, L)
    fmin, f, ∇f, suboptimal = alt_min_lipschitz(hw, hμ, arms, L);
    if r > eps()
        fidx = [j for (idxj,j) in enumerate(suboptimal) if (f[idxj]<fmin+r)]
    elseif abs(r)<eps()
        fidx = [suboptimal[argmin(f)]];
    else
        fidx = suboptimal;
    end
    return f, ∇f, fidx;
end


text:
using Turing: VarName, insdelim
using Turing.VarReplay: is_inside

using Test

varname(s::Symbol)  = nothing, s
varname(expr::Expr) = begin
  # Initialize an expression block to store the code for creating uid
  local sym
  @assert expr.head == :ref "expr needs to be an indexing expression, e.g. :(x[1])"
  indexing_ex = Expr(:block)
  # Add the initialization statement for uid
  push!(indexing_ex.args, quote indexing_list = [] end)
  # Initialize a local container for parsing and add the expr to it
  to_eval = []; pushfirst!(to_eval, expr)
  # Parse the expression and creating the code for creating uid
  find_head = false
  while length(to_eval) > 0
    evaling = popfirst!(to_eval)   # get the current expression to deal with
    if isa(evaling, Expr) && evaling.head == :ref && ~find_head
      # Add all the indexing arguments to the left
      pushfirst!(to_eval, "[", insdelim(evaling.args[2:end])..., "]")
      # Add first argument depending on its type
      # If it is an expression, it means it's a nested array calling
      # Otherwise it's the symbol for the calling
      if isa(evaling.args[1], Expr)
        pushfirst!(to_eval, evaling.args[1])
      else
        # push!(indexing_ex.args, quote pushfirst!(indexing_list, $(string(evaling.args[1]))) end)
        push!(indexing_ex.args, quote sym = Symbol($(string(evaling.args[1]))) end) # store symbol in runtime
        find_head = true
        sym = evaling.args[1] # store symbol in compilation time
      end
    else
      # Evaluting the concrete value of the indexing variable
      push!(indexing_ex.args, quote push!(indexing_list, string($evaling)) end)
    end
  end
  push!(indexing_ex.args, quote indexing = reduce(*, indexing_list) end)
  indexing_ex, sym
end

genvn(sym::Symbol) = VarName(gensym(), sym, "", 1)
genvn(expr::Expr) = begin
  ex, sym = varname(expr)
  VarName(gensym(), sym, eval(ex), 1)
end

space = Set([:x, :y, :(z[1])])
vn1 = genvn(:x)
vn2 = genvn(:y)
vn3 = genvn(:(x[1]))
vn4 = genvn(:(z[1][1]))
vn5 = genvn(:(z[2]))
vn6 = genvn(:z)

@test is_inside(vn1, space)
@test is_inside(vn2, space)
@test is_inside(vn3, space)
@test is_inside(vn4, space)
@test ~is_inside(vn5, space)
@test ~is_inside(vn6, space)


text:
using FluxMNIST
using Flux: throttle
using ArgParse
using Dates

function savemodel(model::FluxMNIST.Model{T}, outputdir::AbstractString) where {T<:AbstractFloat}
    if !isdir(outputdir)
        mkpath(outputdir)
    end
    ts = Dates.format(now(), dateformat"yyyymmddHHMMSS")
    filename = "model-flux-f$(sizeof(T)*8)_$(ts).bson"
    dst = joinpath(outputdir, filename)
    FluxMNIST.savemodel(model, dst)
end

main(args::AbstractString) = main(split(args))
function main(args::Vector{<:AbstractString}=[])
    s = ArgParseSettings()
    s.description="FluxMNIST.jl sample for Float32"
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--batchsize"; arg_type=Int; default=1000; help="minibatch size")
        ("--epochs"; arg_type=Int; default=10; help="number of epochs for training")
        ("--output"; arg_type=String; default="./"; help="output dir")
        ("--fp64"; action=:store_true; help="use `Float64` (default: `Float32`)")
    end
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    parsed_args = parse_args(args, s; as_symbols=true)
    batchsize = parsed_args[:batchsize]
    epochs = parsed_args[:epochs]
    outputdir = parsed_args[:output]
    FType = parsed_args[:fp64] ? Float64 : Float32

    model = FluxMNIST.Model{FType}()
    traindata, tX, tY = FluxMNIST.loadMNIST(FType, batchsize)
    accuracy = FluxMNIST.Accuracy(model.m)
    evalcb = throttle(() -> @show(accuracy(tX, tY)), 10)
    @time FluxMNIST.train!(model, traindata; epochs=epochs, cb=evalcb)
    @show(accuracy(tX, tY))
    # FluxMNIST.savemodel(model)
    savemodel(model, outputdir)
end

abspath(PROGRAM_FILE) == @__FILE__() && main(ARGS)


text:
# Advent of code 2021, day 4

function parseinput()
    paragraphs = split(strip(read("AoCdata/AoC_2021_day4.txt", String)), "\n\n")
    draworder = [parse(Int, num) for num in split(strip(first(paragraphs)), ",")]
    matrices, picks = Matrix{Int}[], Matrix{Int}[]
    for mattxt in paragraphs[(begin + 1):end]
        push!(matrices, reshape([parse(Int, s) for s in split(strip(mattxt))]', (5, 5)))
        push!(picks, zeros(Int, 5, 5))
    end
    return draworder, matrices, picks
end

function part1and2()
    win(mat) = any(sum(mat; dims=1) .== 5) || any(sum(mat; dims=2) .== 5)
    draw, mats, picks = parseinput()
    totalwinners, part1, winners, pastwinners = 0, 0, Int[], Int[]
    for d in draw
        for n in eachindex(mats)
            if (idx = findfirst(j -> j == d, mats[n])) != nothing
                picks[n][idx] = 1
            end
        end
        winners = filter(
            x -> !(x in pastwinners), findall(i -> win(picks[i]), eachindex(mats))
        )
        isempty(winners) && continue
        push!(pastwinners, winners...)
        winmat, winpick = mats[winners[end]], picks[winners[end]]
        unpicksum = sum([winmat[n] for n in eachindex(winmat) if winpick[n] == 0])
        if part1 == 0
            part1 = winners[1]
            println("Part 1: ", unpicksum * d)
        elseif length(pastwinners) == length(mats)  # all have won
            lastwinner = last(pastwinners)
            println("Last winner is $lastwinner: ", unpicksum * d)
            break
        end
    end
end

part1and2()


text:
# private math functions

"""
A helper function for `ldexpk`

First note that `r = (q >> n) << n` clears the lowest n bits of q, i.e. returns 2^n where n is the
largest integer such that q >= 2^n

For numbers q less than 2^m the following code does the same as the above snippet
    `r = ( (q>>v + q) >> n - q>>v ) << n`
For numbers larger than or equal to 2^v this subtracts 2^n from q for q>>n times.

The function returns q(input) := q(output) + offset*r

In the code for ldexpk we actually use
    `m = ( (m>>n + m) >> n -  m>>m ) << (n-2)`.
So that x has to be multplied by u four times `x = x*u*u*u*u` to put the value  of the offset
exponent amount back in.
"""
@inline function _split_exponent(q, n, v, offset)
    m = q >> v
    m = (((m + q) >> n) - m) << (n - offset)
    q = q - (m << offset)
    m, q
end
@inline split_exponent(::Type{Float64}, q::Int) = _split_exponent(q, UInt(9), UInt(31), UInt(2))
@inline split_exponent(::Type{Float32}, q::Int) = _split_exponent(q, UInt(6), UInt(31), UInt(2))

"""
    ldexpk(a, n)

Computes `a × 2^n`.
"""
@inline function ldexpk(x::T, q::Int) where {T<:Union{Float32,Float64}}
    bias = exponent_bias(T)
    emax = exponent_raw_max(T)
    m, q = split_exponent(T, q)
    m += bias
    m = ifelse(m < 0, 0, m)
    m = ifelse(m > emax, emax, m)
    q += bias
    u = integer2float(T, m)
    x = x * u * u * u * u
    u = integer2float(T, q)
    x * u
end

@inline function ldexp2k(x::T, e::Int) where {T<:Union{Float32,Float64}}
    x * pow2i(T, e >> 1) * pow2i(T, e - (e >> 1))
end

@inline function ldexp3k(x::T, e::Int) where {T<:Union{Float32,Float64}}
    reinterpret(T, reinterpret(Unsigned, x) + (Int64(e) << significand_bits(T)) % uinttype(T))
end

# threshold values for `ilogbk`
const threshold_exponent(::Type{Float64}) = 300
const threshold_exponent(::Type{Float32}) = 64

"""
    ilogbk(x) -> Int

Returns the integral part of the logarithm of `|x|`, using 2 as base for the logarithm; in other
words this returns the binary exponent of `x` so that

    x = significand × 2^exponent

where `significand ∈ [1, 2)`.
"""
@inline function ilogbk(d::T) where {T<:Union{Float32,Float64}}
    m = d < T(2)^-threshold_exponent(T)
    d = ifelse(m, d * T(2)^threshold_exponent(T), d)
    q = float2integer(d) & exponent_raw_max(T)
    q = ifelse(m, q - (threshold_exponent(T) + exponent_bias(T)), q - exponent_bias(T))
end

# similar to ilogbk, but argument has to be a normalized float value
@inline function ilogb2k(d::T) where {T<:Union{Float32,Float64}}
    (float2integer(d) & exponent_raw_max(T)) - exponent_bias(T)
end

let
global atan2k_fast
global atan2k

c20d =  1.06298484191448746607415e-05
c19d = -0.000125620649967286867384336
c18d =  0.00070557664296393412389774
c17d = -0.00251865614498713360352999
c16d =  0.00646262899036991172313504
c15d = -0.0128281333663399031014274
c14d =  0.0208024799924145797902497
c13d = -0.0289002344784740315686289
c12d =  0.0359785005035104590853656
c11d = -0.041848579703592507506027
c10d =  0.0470843011653283988193763
c9d  = -0.0524914210588448421068719
c8d  =  0.0587946590969581003860434
c7d  = -0.0666620884778795497194182
c6d  =  0.0769225330296203768654095
c5d  = -0.0909090442773387574781907
c4d  =  0.111111108376896236538123
c3d  = -0.142857142756268568062339
c2d  =  0.199999999997977351284817
c1d =  -0.333333333333317605173818

c9f = -0.00176397908944636583328247f0
c8f =  0.0107900900766253471374512f0
c7f = -0.0309564601629972457885742f0
c6f =  0.0577365085482597351074219f0
c5f = -0.0838950723409652709960938f0
c4f =  0.109463557600975036621094f0
c3f = -0.142626821994781494140625f0
c2f =  0.199983194470405578613281f0
c1f = -0.333332866430282592773438f0

global @inline atan2k_fast_kernel(x::Float64) = @horner x c1d c2d c3d c4d c5d c6d c7d c8d c9d c10d c11d c12d c13d c14d c15d c16d c17d c18d c19d c20d
global @inline atan2k_fast_kernel(x::Float32) = @horner x c1f c2f c3f c4f c5f c6f c7f c8f c9f

@inline function atan2k_fast(y::T, x::T) where {T<:Union{Float32,Float64}}
    q = 0
    if x < 0
        x = -x
        q = -2
    end
    if y > x
        t = x; x = y
        y = -t
        q += 1
    end
    s = y / x
    t = s * s
    u = atan2k_fast_kernel(t)
    t = u * t * s + s
    q * T(PI_2) + t
end

global @inline atan2k_kernel(x::Double{Float64}) = @horner x.hi c1d c2d c3d c4d c5d c6d c7d c8d c9d c10d c11d c12d c13d c14d c15d c16d c17d c18d c19d c20d
global @inline atan2k_kernel(x::Double{Float32}) = dadd(c1f, x.hi * (@horner x.hi c2f c3f c4f c5f c6f c7f c8f c9f))

@inline function atan2k(y::Double{T}, x::Double{T}) where {T<:Union{Float32,Float64}}
    q = 0
    if x < 0
        x = -x
        q = -2
    end
    if y > x
        t = x; x = y
        y = -t
        q += 1
    end

    s = ddiv(y, x)
    t = dsqu(s)
    t = dnormalize(t)

    u = atan2k_kernel(t)

    t = dmul(t, u)
    t = dmul(s, dadd(T(1.0), t))
    T <: Float64 && abs(s.hi) < 1e-200 && (t = s)
    t = dadd(dmul(T(q), MDPI2(T)), t)
    return t
end
end

const under_expk(::Type{Float64}) = -1000.0
const under_expk(::Type{Float32}) = -104f0

@inline function expk_kernel(x::Float64)
    c10 = 2.51069683420950419527139e-08
    c9  = 2.76286166770270649116855e-07
    c8  = 2.75572496725023574143864e-06
    c7  = 2.48014973989819794114153e-05
    c6  = 0.000198412698809069797676111
    c5  = 0.0013888888939977128960529
    c4  = 0.00833333333332371417601081
    c3  = 0.0416666666665409524128449
    c2  = 0.166666666666666740681535
    c1  = 0.500000000000000999200722
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9 c10
end

@inline function  expk_kernel(x::Float32)
    c5 = 0.00136324646882712841033936f0
    c4 = 0.00836596917361021041870117f0
    c3 = 0.0416710823774337768554688f0
    c2 = 0.166665524244308471679688f0
    c1 = 0.499999850988388061523438f0
    return @horner x c1 c2 c3 c4 c5
end

@inline function expk(d::Double{T}) where {T<:Union{Float32,Float64}}
    q = round(T(d) * T(MLN2E))
    qi = unsafe_trunc(Int, q)

    s = dadd(d, -q * L2U(T))
    s = dadd(s, -q * L2L(T))

    s = dnormalize(s)

    u = expk_kernel(T(s))

    t = dadd(s, dmul(dsqu(s), u))
    t = dadd(T(1.0), t)
    u = ldexpk(T(t), qi)

    (d.hi < under_expk(T)) && (u = T(0.0))
    return u
end

@inline function expk2_kernel(x::Double{Float64})
    c11 = 0.1602472219709932072e-9
    c10  = 0.2092255183563157007e-8
    c9  = 0.2505230023782644465e-7
    c8  = 0.2755724800902135303e-6
    c7  = 0.2755731892386044373e-5
    c6  = 0.2480158735605815065e-4
    c5  = 0.1984126984148071858e-3
    c4  = 0.1388888888886763255e-2
    c3  = 0.8333333333333347095e-2
    c2  = 0.4166666666666669905e-1
    c1 = 0.1666666666666666574e0
    u = @horner x.hi c2 c3 c4 c5 c6 c7 c8 c9 c10 c11
    return dadd(dmul(x, u), c1)
end

@inline function  expk2_kernel(x::Double{Float32})
    c5 = 0.1980960224f-3
    c4 = 0.1394256484f-2
    c3 = 0.8333456703f-2
    c2 = 0.4166637361f-1
    c1 = 0.166666659414234244790680580464f0
    u = @horner x.hi c2 c3 c4 c5
    return dadd(dmul(x, u), c1)
end

@inline function expk2(d::Double{T}) where {T<:Union{Float32,Float64}}
    q = round(T(d) * T(MLN2E))
    qi = unsafe_trunc(Int, q)

    s = dadd(d, -q * L2U(T))
    s = dadd(s, -q * L2L(T))

    t = expk2_kernel(s)
    t = dadd(dmul(s, t), T(0.5))
    t = dadd(s, dmul(dsqu(s), t))
    t = dadd(T(1.0), t)

    t = Double(ldexp2k(t.hi, qi), ldexp2k(t.lo, qi))

    (d.hi < under_expk(T)) && (t = Double(T(0.0)))
    return t
end

@inline function logk2_kernel(x::Float64)
    c8 = 0.13860436390467167910856
    c7 = 0.131699838841615374240845
    c6 = 0.153914168346271945653214
    c5 = 0.181816523941564611721589
    c4 = 0.22222224632662035403996
    c3 = 0.285714285511134091777308
    c2 = 0.400000000000914013309483
    c1 = 0.666666666666664853302393
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8
end

@inline function logk2_kernel(x::Float32)
    c4 = 0.240320354700088500976562f0
    c3 = 0.285112679004669189453125f0
    c2 = 0.400007992982864379882812f0
    c1 = 0.666666686534881591796875f0
    return @horner x c1 c2 c3 c4
end

@inline function logk2(d::Double{T}) where {T<:Union{Float32,Float64}}
    e  = ilogbk(d.hi * T(1.0/0.75))
    m  = scale(d, pow2i(T, -e))

    x  = ddiv(dsub2(m, T(1.0)), dadd2(m, T(1.0)))
    x2 = dsqu(x)

    t  = logk2_kernel(x2.hi)

    s = dmul(MDLN2(T), T(e))
    s = dadd(s, scale(x, T(2.0)))
    s = dadd(s, dmul(dmul(x2, x), t))
    return s
end

@inline function logk_kernel(x::Double{Float64})
    c10 = 0.116255524079935043668677
    c9 = 0.103239680901072952701192
    c8 = 0.117754809412463995466069
    c7 = 0.13332981086846273921509
    c6 = 0.153846227114512262845736
    c5 = 0.181818180850050775676507
    c4 = 0.222222222230083560345903
    c3 = 0.285714285714249172087875
    c2 = 0.400000000000000077715612
    c1 = Double(0.666666666666666629659233, 3.80554962542412056336616e-17)
    dadd2(dmul(x, @horner x.hi c2 c3 c4 c5 c6 c7 c8 c9 c10), c1)
end

@inline function logk_kernel(x::Double{Float32})
    c4 = 0.240320354700088500976562f0
    c3 = 0.285112679004669189453125f0
    c2 = 0.400007992982864379882812f0
    c1 = Double(0.66666662693023681640625f0, 3.69183861259614332084311f-9)
    dadd2(dmul(x, @horner x.hi c2 c3 c4), c1)
end

@inline function logk(d::T) where {T<:Union{Float32,Float64}}
    o = d < floatmin(T)
    o && (d *= T(Int64(1) << 32) * T(Int64(1) << 32))

    e  = ilogb2k(d * T(1.0/0.75))
    m  = ldexp3k(d, -e)

    o && (e -= 64)

    x  = ddiv(dsub2(m, T(1.0)), dadd2(T(1.0), m))
    x2 = dsqu(x)

    t  = logk_kernel(x2)

    s = dmul(MDLN2(T), T(e))
    s = dadd(s, scale(x, T(2.0)))
    s = dadd(s, dmul(dmul(x2, x), t))
    return s
end


text:
function STMEquations(arg0::CR3BPSystem)
    return STMEquations((CR3BPSystem,), arg0)
end

function compute_derivatives(obj::STMEquations, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "computeDerivatives", Vector{jdouble}, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_name(obj::STMEquations)
    return jcall(obj, "getName", JString, ())
end

function get_state_transition_matrix(obj::STMEquations, arg0::SpacecraftState)
    return jcall(obj, "getStateTransitionMatrix", RealMatrix, (SpacecraftState,), arg0)
end

function init(obj::AdditionalEquations, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function set_initial_phi(obj::STMEquations, arg0::SpacecraftState)
    return jcall(obj, "setInitialPhi", SpacecraftState, (SpacecraftState,), arg0)
end


text:
#=
Section Allocation

   Suppose you have $n$ students in a class who need to be assigned to
      $m$ discussion sections. Each student needs to be assigned to exactly
      one section. Each discussion section should have between 6 and 10
      students. Suppose an $n \times m$ preference matrix $P$ is given,
      where $P_{ij}$ gives student $i$'s ranking for section $j$
      (1 would mean it is the student's top choice, 10,000 or a large number
      would mean the student can not attend that section).

   The goal will be to get an allocation matrix $X$, where $X_{ij} = 1$
      if student $i$ is assigned to section $j$ and $0$ otherwise.

  source: https://github.com/JuliaOpt/Convex.jl/blob/master/examples/
                             section_allocation/section_allocation.ipynb
=#

using Convex, GLPKMathProgInterface

include("data.jl"); # data.jl has our preference matrix, P

X = Variable(size(P), :Bin)

#=
   We want every student to be assigned to exactly one section.
   So, every row must have exactly one non-zero entry.
   In other words, the sum of all the columns for every row is 1.
   We also want each section to have between 6 and 10 students,
       so the sum of all the rows for every column should be between
       these constraints = [sum(X, 2) == 1, sum(X, 1) <= 10, sum(X, 1) >= 6]

   Our objective is simple sum(X .* P), which can be more efficiently
       represented as vec(X)' * vec(P).
   Since each entry of X is either 0 or 1, this is basically summing up
       the rankings of students that were assigned to them.
   If all students got their first choice, this value will be the number
       of students since the ranking of the first choice is 1.
=#

p = minimize(vec(X)' * vec(P), constraints)

solve!(p, GLPKSolverMIP())


text:
 # function to allocate data using linked lists

  function linkedlist(box, DATA)                        # Box = points, side = box's side, cutoff
    Nc = trunc(Int64 ,DATA.side / DATA.cutoff)          # Number of cells
    first_atom = zeros(Int64,Nc, Nc)                    # vector that contain the first atom of each cell
    next_atom  = zeros(Int64, DATA.N)                   # vector that that the value of a atom position correponds to the next atom's position

    for iat in 1:DATA.N
      icell = trunc(Int64, box[iat][1]/DATA.cutoff) + 1      #  classification of each particle in a specific cell
      jcell = trunc(Int64, box[iat][2]/DATA.cutoff) + 1      #

      if icell==0
        icell = icell + 1
      elseif icell< 0
        while icell <= 0
          icell += 1
         end
      end

      if jcell ==0
        jcell = jcell + 1
      elseif jcell < 0
        while jcell <= 0
          jcell = jcell + 1
        end
      end

      next_atom[iat] = first_atom[icell, jcell]         #    saving the next atom
      first_atom[icell, jcell] = iat                    #    saving the first atom
    end

    return Nc, first_atom, next_atom
  end


text:
include(joinpath("linReg", "chosenVariables_test.jl"))
include(joinpath("linReg", "research_test.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
using Test
using Dex
using HTTP
using JSON

function createconfig(workdir::String)
    cfg = """issuer: http://127.0.0.1:9999/dex
storage:
  type: sqlite3
  config:
    file: $workdir/conf/dex.db
web:
  http: 127.0.0.1:8080
logger:
  level: "debug"
  format: "text"
staticClients:
- id: local
  redirectURIs:
    - 'http://127.0.0.1:8888/auth/login'
  name: 'test client'
  secret: '000000000000000000000000'
oauth2:
  skipApprovalScreen: true
enablePasswordDB: true
staticPasswords:
- email: "admin@example.com"
  # bcrypt hash of the string "password"
  hash: "\$2a\$10\$2b2cU8CPhOTaGrs1HRQuAueS7JTT5ZHsHSzYiFPm1leZck7Mc8T4W"
  username: "admin"
  userID: "08a8684b-db88-4b73-90a9-3cd1661f5466"
"""
    cfgfile = joinpath(workdir, "test.cfg")

    open(cfgfile, "w") do f
        println(f, cfg)
    end

    cfgfile
end

function test_dex_config()
    body = ""

    @info("waiting for dex to come up")
    # try for 10 secs
    while isempty(body)
        sleep(2)
        try
            resp = HTTP.get("http://127.0.0.1:8080/dex/.well-known/openid-configuration")
            body = String(resp.body)
        catch ex
            @info("dex not ready yet...")
        end
    end

    isempty(body) && error("dex did not come up")

    spec = JSON.parse(body)
    @info("testing dex openid configuration")
    @test spec["issuer"] == "http://127.0.0.1:9999/dex"
    @test spec["userinfo_endpoint"] == "http://127.0.0.1:9999/dex/userinfo"
    @test spec["token_endpoint"] == "http://127.0.0.1:9999/dex/token"
    @test spec["jwks_uri"] == "http://127.0.0.1:9999/dex/keys"
    nothing
end

function test()
    workdir = mktempdir()
    mkpath(workdir)
    cfgfile = createconfig(workdir)
    dex = DexCtx(workdir)

    @info("setting up Dex", workdir)
    setup(dex, cfgfile)
    @test isfile(Dex.conffile(dex))

    @info("starting Dex")
    start(dex)
    sleep(2)
    @test isfile(Dex.pidfile(dex))
    @test isrunning(dex)

    test_dex_config()

    @info("restarting Dex")
    restart(dex; delay_seconds=2)
    sleep(2)
    @test isfile(Dex.pidfile(dex))
    @test isrunning(dex)

    test_dex_config()

    @info("stopping Dex")
    stop(dex)
    sleep(2)
    @test !isfile(Dex.pidfile(dex))
    @test !isrunning(dex)
    @test isfile(Dex.logfile(dex))

    @test_throws Exception setup(dex, cfgfile)
    @test nothing === setup(dex, cfgfile; force=true)

    @info("cleaning up")
    rm(workdir; recursive=true, force=true)
    @info("done")

    nothing
end

test()


text:
module TestReduce

using FLoops
using Test
using Transducers: SplitBy

function test_sum()
    @floop for (x, y) in zip(1:3, 1:2:6)
        a = x + y
        b = x - y
        @reduce(s += a, t += b)
    end
    @test (s, t) === (15, -3)
end

function test_sum_with_init()
    @floop for (x, y) in zip(1:3, 1:2:6)
        a = x + y
        b = x - y
        @reduce(s = 0im + a, t = 0im + b)
    end
    @test (s, t) === (15 + 0im, -3 + 0im)
end

function test_sum_nested_loop()
    @floop for x in 1:3, y in 1:2:6
        a = x + y
        b = x - y
        @reduce(s += a, t += b)
    end
    @test (s, t) === (45, -9)
end

function test_findmax()
    @floop for (i, v) in pairs([0, 1, 3, 2])
        y = 2v
        @reduce() do (ymax = -Inf; y), (imax = 0; i)
            if isless(ymax, y)
                ymax = y
                imax = i
            end
        end
    end
    @test (ymax, imax) == (6, 3)
end

function test_findminmax()
    @floop for (i, v) in pairs([0, 1, 3, 2])
        y = 2v
        @reduce() do (ymax = -Inf; y), (imax = 0; i)
            if isless(ymax, y)
                ymax = y
                imax = i
            end
        end
        @reduce() do (ymin = Inf; y), (imin = 0; i)
            if isless(y, ymin)
                ymin = y
                imin = i
            end
        end
    end
    @test (ymax, imax) == (6, 3)
    @test (ymin, imin) == (0, 1)
end

function test_break()
    @floop for x in 1:10
        @reduce(s += x)
        x == 3 && break
    end
    @test s == 6
end

function test_findmax_with_filtering()
    @floop for (i, v) in pairs([0, 1, 4, 3, 1])
        if isodd(v)
            @reduce() do (vmax = -Inf; v), (imax = 0; i)
                if isless(vmax, v)
                    vmax = v
                    imax = i
                end
            end
        end
    end
    @test (vmax, imax) == (3, 4)
end

function test_nonunique_arguments()
    @testset "`@reduce(a += x, b *= x)`" begin
        xs = 1:10
        @floop for x in xs
            @reduce(a += x, b *= x)
        end
        @test (sum(xs), prod(xs)) == (a, b)
    end

    @testset "`@reduce(a = 0 + x, b = 1 * x)`" begin
        xs = 1:10
        @floop for x in xs
            @reduce(a = 0 + x, b = 1 * x)
        end
        @test (sum(xs), prod(xs)) == (a, b)
    end

    @testset "`@reduce a = 0 + x; @reduce b = 1 * x`" begin
        xs = 1:10
        @floop for x in xs
            @reduce a = 0 + x
            @reduce b = 1 * x
        end
        @test (sum(xs), prod(xs)) == (a, b)
    end
end

function floop_with_init_binops(xs, ex = nothing)
    @floop ex for x in xs
        @reduce(s = 0 + x, c = 0 + 1)
    end
    return f() = (s, c)
end

function two_floops(xs, ex = nothing)
    @floop ex for x in xs
        @reduce(s = 0 + x, c = 0 + 1)
    end
    @floop ex for x in xs
        @reduce(v = 0 + (x - (s / c))^2)
    end
    return v / c
end

function test_empty()
    @test !FLoops.has_boxed_variables(floop_with_init_binops(1:10))
    @test two_floops(1:8) == 5.25
end

function sum_arrays(arrays, ex = nothing)
    @floop ex for x in arrays
        # @reduce(s = zero(y) .+ y)  # TODO
        @reduce() do (s = zero(x); x)
            s .+= x
        end
    end
    return s
end

function test_at_reduce_init_scope()
    @test sum_arrays([[1], [2], [3]]) == [sum(1:3)]
    @test sum_arrays([[1], [2], [3]], SequentialEx()) == [sum(1:3)]
    @testset "let" begin
        local y = 0
        @floop for x in 1:10
            let y = 123
            end
            @reduce(s = y + x)
        end
        @test s == sum(1:10)
    end
end

function maximum_partition_length(f, xs, ex = nothing)
    @floop ex for chunk in SplitBy(f)(xs)
        @reduce m = max(typemin(Int), length(chunk))
    end
    return m
end

function test_SplitBy()
    @test maximum_partition_length(isodd, 1:10) == 1
    @test maximum_partition_length(isodd, 1:10, SequentialEx()) == 1
    @test maximum_partition_length(==(7), 1:10) == 6
    @test maximum_partition_length(==(7), 1:10, SequentialEx()) == 6
    xs = rand(0:9, 1000)
    @test maximum_partition_length(iszero, xs) ==
          maximum_partition_length(iszero, xs, SequentialEx())
end

function sum_halved_arrays(arrays, ex = nothing)
    @floop ex for x in arrays
        @init y = zero(x)
        y .= x .÷ 2
        @reduce(s = 0 + sum(y))
    end
    return s
end

function two_inits(arrays, ex = nothing)
    @floop ex for x in arrays
        @init begin
            y = zero(x)
            z = similar(y)
        end
        y .= x .÷ 2
        z .= 2 .* y
        r = sum(a * b for (a, b) in zip(y, z))
        @reduce(s = 0 + r)
    end
    return s
end

function two_inits2(arrays, ex = nothing)
    @floop ex for x in arrays
        @init y = zero(x)
        @init z = similar(y)
        y .= x .÷ 2
        z .= 2 .* y
        r = sum(a * b for (a, b) in zip(y, z))
        @reduce(s = 0 + r)
    end
    return s
end

function test_at_init()
    arrays = [[1, 2], [3, 4], [5, 6], [7, 8]]
    desired = sum(sum(x .÷ 2) for x in arrays)
    @test sum_halved_arrays(arrays) == desired
    @test sum_halved_arrays(arrays, SequentialEx()) == desired

    desired = sum((c ÷ 2) * (2 * (c ÷ 2)) for x in arrays for c in x)
    @test two_inits(arrays) == desired
    @test two_inits(arrays, SequentialEx()) == desired
    @test two_inits2(arrays) == desired
    @test two_inits2(arrays, SequentialEx()) == desired
end

function test_just_at_init()
    n = 10
    dest = zeros(Int, n)
    @floop for (i, x) in pairs(1:n)
        @init y = zeros(Int, 3)
        y .= (x, 2x, 3x)
        dest[i] = y[1]
    end
    @test dest == 1:n
end

function probe_init(xs, ex = nothing)
    ninit = Threads.Atomic{Int}(0)
    @floop ex for x in xs
        @init y = begin
            Threads.atomic_add!(ninit, 1)
            nothing
        end
    end
    return ninit[]
end

function test_at_init_called_once()
    @testset "default" begin
        @test probe_init(1:(Threads.nthreads() * 100)) == Threads.nthreads()
    end
    @testset "SequentialEx" begin
        @test probe_init(1:(Threads.nthreads() * 100), SequentialEx()) == 1
    end
end

function init_is_private(xs, ex = nothing)
    @floop ex for x in xs
        @init p = []
        push!(p, x)
        y = pop!(p)
        @reduce s += y
    end
    return (sum = s, isdefined_p = @isdefined(p))
end

function test_init_is_private()
    @test init_is_private(1:10) == (sum = sum(1:10), isdefined_p = false)
    @test init_is_private(1:10, SequentialEx()) == (sum = sum(1:10), isdefined_p = false)
end

function test_unprocessed_at_reduce()
    err = try
        @reduce(s += y, p *= y)
        nothing
    catch err
        err
    end
    @test err isa FLoops.ReduceOpSpec
    @test occursin("used outside `@floop`", sprint(showerror, err))
end

function test_unprocessed_at_init()
    err = try
        @init x = 0
        nothing
    catch err
        err
    end
    @test err isa FLoops.InitSpec
    @test occursin("used outside `@floop`", sprint(showerror, err))
end

function test_invalid_at_init()
    @testset "toplevel" begin
        err = try
            @eval @init(a)
            nothing
        catch err
            err
        end
        @test err isa Exception
        @test occursin("requires an assignment", sprint(showerror, err))
    end
    @testset "non assignment (+=)" begin
        err = try
            @eval @init a += 1
            nothing
        catch err
            err
        end
        @test err isa Exception
        @test occursin("requires an assignment", sprint(showerror, err))
    end
    @testset "non assignment (constant)" begin
        err = try
            @eval @init begin
                a = 1
                1
            end
            nothing
        catch err
            err
        end
        @test err isa Exception
        @test occursin("requires an assignment", sprint(showerror, err))
    end
end

function test_duplicated_accumulators()
    err = try
        @eval @macroexpand @floop for x in xs
           @reduce(y += x, y += x)
        end
        nothing
    catch err
        err
    end
    @test err isa Exception
    @test occursin("`y` used 2 times", sprint(showerror, err))
end

end  # module


text:
using Bonsai, URIs, StructTypes
using JSON3
using StructTypes: @Struct

app = App()
app.docs = "/docs"

@Struct struct Id
    id::Int
end

"""
limit - How many items to return at one time (max 100)
offset - Offset to start returning pets from
"""
struct Limit
    limit::Int
    offset::Union{Int,Missing}
end

@Struct Limit

@Struct struct Pet
    id::Int
    name::String
    tag::String
end

"""
code - a numerical error code
"""
struct Error
    code::Int
    type::String
end

@Struct Error

@Struct struct Next
    x_next::String
end

pets = Dict(
    1 => Pet(
        1,
        "bobby",
        "dog"
    )
)

app.post("/pets/{id}") do stream
    pet = Bonsai.read(stream, Body(Pet))
    pets[pet.id] = pet
    @info "created pet" pet = pet
    Bonsai.write(stream, "ok")
end

app.get("/pets/{id}") do stream
    id = Bonsai.read(stream, PathParams(Id)).id
    if haskey(pets, id)
        pet::Pet = pets[id]
        Bonsai.write(stream, pet, ResponseCodes.Ok())
    else
        err = Error(
            1,
            "Pet not found :("
        )
        Bonsai.write(stream, err, ResponseCodes.NotFound())
    end
end

app.get("/pets") do stream
    query = Bonsai.read(stream, Query(Limit))
    l = Pet[]

    for (i, pet) in values(pets)

        if !ismissing(query.offset) & i < query.offset
            continue
        end

        push!(l, pet)
        if i > query.limit
            break
        end
    end

    Bonsai.write(stream, Headers(Next("/pets?limit=$(query.limit+1)&offset=$(query.offset)")))
    Bonsai.write(stream, Body(l))
end

app.get("**") do stream, next
    try
        next(stream)
    catch e
        if e isa NoHandler
            Bonsai.write(stream, Body("Not found"), ResponseCodes.NotFound())
        else
            @error "error" e = e
            bt = catch_backtrace()
            showerror(stderr, e, bt)
            Bonsai.write(stream, Body(repr(e)), ResponseCodes.InternalServerError())
        end
    end
end

app.get(app.docs) do stream
    @info "docs"
    html = Bonsai.create_docs_html(app)
    Bonsai.write(stream, Body(html))
    Bonsai.write(stream, Headers(content_type="text/html; charset=UTF-8"))
end

app.get(app.docs * ".json") do stream
    Bonsai.write(stream, Body(OpenAPI(app)))
end

# JSON3.write("openapi.json", OpenAPI(app))

start(app, port=10004)

stop(app)

html = Bonsai.create_docs_html(app)
print(html)

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
module ShellFunctions

using Parameters: @with_kw

include("options.jl") # provides Options module
using .Options

include("echo.jl")
include("head.jl")
include("cat.jl")

using DelimitedFiles
include("awk.jl")

end # module


text:
using InteractiveLearning
using Test

@testset "InteractiveLearning.jl" begin
    # Write your tests here.
end


text:
module Signatures

export searchname, signature, fullyqualifiedname, MethodSignature

using JavaCall.JNI

struct MethodSignature
    returntype::Any
    parametertypes::Vector{Any}
end

searchname(name::String) = replace(name, "." => "/")
searchname(name::Symbol) = searchname(string(name))

fullyqualifiedname(name::String) = string("L", searchname(name), ";")

signature(::Type{jboolean}) = "Z"
signature(::Type{jbyte}) = "B"
signature(::Type{jchar}) = "C"
signature(::Type{jshort}) = "S"
signature(::Type{jint}) = "I"
signature(::Type{jlong}) = "J"
signature(::Type{jfloat}) = "F"
signature(::Type{jdouble})= "D"
signature(::Type{jvoid}) = "V"

signature(::Type{Array{T,N}}) where {T,N} = string("[" ^ N, signature(T))

signature(sym::Symbol) = fullyqualifiedname(string(sym))

signature(s::MethodSignature) =
    string("(", map(signature, s.parametertypes)..., ")", signature(s.returntype))

end


text:
# sample sequences of ergms with different parameters' values from R package ergm
# and test the PseudoLikelihoodScoreDrivenERGM filter
using Utilities,AReg,StaticNets,DynNets , JLD,MLBase,StatsBase,CSV, RCall
 using PyCall; pygui(); using PyPlot
 using ErgmRcall

#load simulations results divided in chunks
T = 500
      tot_Nsample = 200
      N_each_chunk = 10
      N = 50
      UM = [-3 , -0.5]
      B_Re  = 0.95
      A_Re  = 0.05
      indTvPar = trues(2)
      indTargPar = falses(2)
      Nsample = N_each_chunk
      # Save all Sampled Data
      load_fold = "./data/estimatesTest/sdergmTest/distrib_static_pars/"
      n_chunk = 1
      data = load(load_fold*"sdergm_dgp__$(N)_T_$(T)_Sample_$(Nsample)_UM_$(UM)_intTV_$(Int.(indTvPar[:]))_intTarg_$(Int.(indTargPar[:]))_N_each_chunk_$(N_each_chunk)_chunk_$(n_chunk).jld" )
      obsT_all = data["obsT_all"][1:end-1]
      estPar_all = data["estPar_all"][1:end-1]
      for n_chunk = 2:Int(tot_Nsample/N_each_chunk)
            data = load(load_fold*"sdergm_dgp__$(N)_T_$(T)_Sample_$(Nsample)_UM_$(UM)_intTV_$(Int.(indTvPar[:]))_intTarg_$(Int.(indTargPar[:]))_N_each_chunk_$(N_each_chunk)_chunk_$(n_chunk).jld" )
            global obsT_all = vcat(obsT_all, data["obsT_all"][1:end-1])
            global estPar_all = vcat(estPar_all, data["estPar_all"][1:end-1])
      end

      tot_Nsample = size(estPar_all)[1]
      W_est_all = reduce(vcat, [[estPar_all[i][1][1]  estPar_all[i][2][1]] for i=1:tot_Nsample])
      B_est_all = reduce(vcat, [[estPar_all[i][1][2]  estPar_all[i][2][2]] for i=1:tot_Nsample])
      A_est_all = reduce(vcat, [[estPar_all[i][1][3]  estPar_all[i][2][3]] for i=1:tot_Nsample])
      W_dgp = UM.*(1-B_Re)
      #plt.close()
      fig = figure( figsize=(10,10)) # Create a new blank figure
      indPar = 1
      subplot(321) # Create the 1st axis of a 2x2 arrax of axes
      grid("on") # Create a grid on the axis
      plt.hist(W_est_all[:,indPar])
      legend("W estimates") # Give the most recent axis a titl
      plt.axvline(x=W_dgp[indPar], color=[1,0,0])
      subplot(323) # Create the 1st axis of a 2x2 arrax of axes
      grid("on") # Create a grid on the axis
      plt.hist(B_est_all[:,indPar])
      legend("B estimates") # Give the most recent axis a titl
      plt.axvline(x=B_Re, color=[1,0,0])
      subplot(325) # Create the 4th axis of a 2x2 arrax of axes
      grid("on") # Create a grid on the axis
      plt.hist(A_est_all[:,1])
      plt.axvline(x=A_Re, color=[1,0,0])
      indPar = 2
      subplot(322) # Create the 1st axis of a 2x2 arrax of axes
      grid("on") # Create a grid on the axis
      plt.hist(W_est_all[:,indPar])
      legend("A estimates") # Give the most recent axis a titl
      plt.axvline(x=W_dgp[indPar], color=[1,0,0])
      subplot(324) # Create the 1st axis of a 2x2 arrax of axes
      grid("on") # Create a grid on the axis
      plt.hist(B_est_all[:,indPar])
      plt.axvline(x=B_Re, color=[1,0,0])
      subplot(326) # Create the 4th axis of a 2x2 arrax of axes
      plt.title("A estimates") # Give the most recent axis a title
      grid("on") # Create a grid on the axis
      plt.hist(A_est_all[:,indPar])
      plt.axvline(x=A_Re, color=[1,0,0])
      suptitle("T = $(T), N= $(N), Unc Means = $(UM), targeting = $(indTargPar)")

# Compute the hessian of the loglikelihood for each sample
using Statistics, ForwardDiff, LinearAlgebra
      Nsample = tot_Nsample
      cov_hat_all = fill( zeros(6,6),Nsample)
      hess_hat_all = fill( zeros(6,6),Nsample)
      for n=1:Nsample
            print(n)
            model = DynNets.SdErgmDirBinGlobalPseudo(obsT_all[1],DynNets.fooGasPar,indTvPar,"")
            tmpParStat,~ = DynNets.estimate(model;indTvPar = BitArray(undef,2),changeStats_T =obsT_all[n][1:5])
            ftot_0 = [tmpParStat[1][1], tmpParStat[2][1]]
            function likeFun(vecReGasParAll::Array{<:Real,1})
                  oneInADterms  = (maxLargeVal + vecReGasParAll[1])/maxLargeVal
                  foo,loglikelValue = DynNets.score_driven_filter_or_dgp(model,vecReGasParAll,indTvPar;
                                                  obsT = obsT_all[n],ftot_0 = ftot_0 .* oneInADterms)
                  #println(vecReGasPar)
                  return - loglikelValue
            end
            vecAllParGasHat = [W_est_all[n,1], B_est_all[n,1], A_est_all[n,1], W_est_all[n,2], B_est_all[n,2], A_est_all[n,2]]
            likeFun(vecAllParGasHat)
            print( ForwardDiff.gradient(likeFun,vecAllParGasHat))
            hess_hat =  ForwardDiff.hessian(likeFun,vecAllParGasHat)
            cov_hat = round.(pinv(hess_hat),digits = 8)
            cov_hat_all[n] = cov_hat
            hess_hat_all[n] = hess_hat
      end

# Compute the covariance of different estimates
par_est_all = [W_est_all[:,1] B_est_all[:,1] A_est_all[:,1] W_est_all[:,2] B_est_all[:,2] A_est_all[:,2]]
      cov_hat_all2 = cov(par_est_all)
      a =  [cov_hat_all[n][diagind(cov_hat_all[n])] for n =1:Nsample]
      cov_hat_all_diag = permutedims(reshape(hcat(a...), (length(a[1]), length(a))))

      fig, axs = plt.subplots(3,2, figsize=(10,10)) # Create a new blank figure
      for indPar = 1:3
       # Create the 1st axis of a 2x2 arrax of axes
      grid("on") # Create a grid on the axis
      axs[indPar,1].hist(cov_hat_all_diag[cov_hat_all_diag[:,indPar].>0,indPar])
      axs[indPar,1].axvline(x=cov_hat_all2[indPar,indPar], color=[1,0,0])
      grid("on") # Create a grid on the axis
      axs[indPar,2].hist(cov_hat_all_diag[cov_hat_all_diag[:,indPar+3].>0,indPar+3])
      axs[indPar,2].axvline(x=cov_hat_all2[indPar+3,indPar+3], color=[1,0,0])
      end
      suptitle("T = $(T), N= $(N), Unc Means = $(UM), targeting = $(indTargPar)")

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using JLD2,CSV
#Load CsV files and save as julia data
T = 59

assetHeader = "names"
for a=2:21
    assetHeader = [assetHeader;"$(a)"]
end

matrT = [ zeros(Int,2,2) for t=1:T]
equityT = [ zeros(Int,1) for t=1:T]
IDsT = [ zeros(Int,1) for t=1:T]
for t=1:T
    matfilename = "/home/Domenico/Dropbox/Network_Ensembles/data/Real_networks_raw_data/bipartite/dataArtAssessing/US_com_banks/matr_$(t).csv"
    matframe = CSV.read(matfilename;header = assetHeader)
    mat = Array{Int,2}(matframe)
    matrT[t] = mat[:,2:end]
    IDsT[t] = mat[:,1]
    eqfilename = "/home/Domenico/Dropbox/Network_Ensembles/data/Real_networks_raw_data/bipartite/dataArtAssessing/US_com_banks/equities_$(t).csv"
    equframe = CSV.read(eqfilename; header = ["Equity"])
    eqvec = Array{Int,2}(equframe)
    equityT[t] = squeeze(eqvec,2)
end

save_path = "/home/Domenico/Dropbox/Network_Ensembles/data/Real_networks_raw_data/bipartite/dataArtAssessing/USComBanksMatsAndEquities.jld"
@save(save_path,equityT,IDsT,matrT)
##

#


text:
using Test
using StaticArrays
using Random
using Nystrom
using LinearAlgebra
Random.seed!(1)

@testset "Density test" begin
    Geometry.clear_entities!()
    pde   = Helmholtz(;dim=3,k=1)
    Ω     = ParametricSurfaces.Sphere()
    Γ     = boundary(Ω) |> Geometry.Domain
    M     = ParametricSurfaces.meshgen(Γ,(10,10))
    mesh  = NystromMesh(view(M,Γ),order=1)
    σ     = Density(target->norm(coords(target)),mesh)
    @test eltype(σ) == Float64
    σ     = Density(mesh) do target
        x = coords(target)
        exp(im*2*norm(x))
    end
    @test eltype(σ) == ComplexF64
end

@testset "TangentialDensity test" begin
    Geometry.clear_entities!()
    Ω   = ParametricSurfaces.Sphere(;radius=1) |> Geometry.Domain
    Γ   = boundary(Ω)
    M   = ParametricSurfaces.meshgen(Γ,(2,2))
    mesh = NystromMesh(view(M,Γ),order=2)
    pde = Elastostatic(;dim=3,μ=2,λ=3)
    T = Nystrom.default_density_eltype(pde)
    xout = SVector(3,3,3)
    c = rand(T)

    σ = Density(mesh) do target  # density defined with a tangential field
        x = coords(target)
        cross(normal(target), SingleLayerKernel(pde)(xout,target)*c)
    end
    tan_σ = TangentialDensity(σ)
    @test Density(tan_σ) ≈ σ

    ncross_σ = ncross(σ)
    ncross_tan_σ = ncross(tan_σ)
    @test TangentialDensity(ncross_σ) ≈ ncross_tan_σ
    @test Density(ncross_tan_σ) ≈ ncross_σ
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
export exact_free_energy

function bitarray(v::Vector{T}, num_bit::Int)::BitArray{2} where T<:Number
    #ba = BitArray{2}(0, 0)
    ba = BitArray(undef, 0, 0)
    ba.len = 64*length(v)
    ba.chunks = UInt64.(v)
    ba.dims = (64, length(v))
    view(ba, 1:num_bit, :)
end

function exact_free_energy(K)
    nbits = size(K, 1)
    configs = bitarray(collect(0:(1<<nbits-1)), nbits)
    energy = sum(configs .* (K*configs), dims=1)
    -log( sum( exp.(-energy)))
end


text:
#Testbed for synthetic spike recovery from example background
using BioBackgroundModels, BioMotifInference, Random, Distributed, Distributions, Serialization
Random.seed!(786)
#CONSTANTS
no_obs=500
obsl=100:200

const folder_path="/bench/PhD/NGS_binaries/BBM/refined_folders"

report_folders=deserialize(folder_path)

bhmm_vec=[BHMM([0.40027008799648645, 0.1997894691941005, 0.006556173444055311, 0.3381222117369811, 0.053302496489949405, 0.0019595611384158307], [0.9910385961237116 0.001055084535731191 0.0030117366072612726 0.0025443692116936516 1.875366316590218e-83 0.002350213521605879; 0.003592562092379169 0.6212087435679351 1.4862557155404555e-91 0.3738436265294334 0.001355067810251799 5.674610926019336e-174; 0.03314627345971247 0.00026886477654034587 0.9628685354878737 3.677787936178394e-5 0.00367954839651186 5.0e-324; 0.0002903769304330447 0.3737391779744011 1.584645416437449e-35 0.6225645652084915 0.002839228767872827 0.0005666511188010834; 0.0005590565306708014 0.0064195824547637215 0.004061509889975232 0.01962586251822172 0.9687912312979926 0.0005427573083774576; 0.0174697691356862 2.0648278441425858e-20 8.008498087945918e-218 0.004151185943888868 0.0013251598095032858 0.9770538851109234], [Categorical([0.2565257692018922, 0.24506145230281004, 0.2649485969649606, 0.23346418153033835]), Categorical([0.052491168994042534, 0.0987036876775239, 0.2383241923003866, 0.610480951028046]), Categorical([0.08371833024203595, 0.07204591736665274, 0.3854013598343639, 0.4588343925569475]), Categorical([0.5954709485303519, 0.2497223085607694, 0.09528006039357192, 0.0595266825153071]), Categorical([0.4744103622376032, 0.04038264172760258, 0.04588736323302444, 0.4393196328017719]), Categorical([0.41173384083242826, 0.44270375576246973, 0.030612125665857307, 0.11495027773924743])])]

bhmm_dist=Categorical([1.])

struc_sig_1=[.1 .7 .1 .1
           .1 .1 .1 .7
           .1 .7 .1 .1]
periodicity=8
struc_frac_obs=.75

tata_box=[.05 .05 .05 .85
          .85 .05 .05 .05
          .05 .05 .05 .85
          .85 .05 .05 .05
          .425 .075 .075 .425
          .85 .05 .05 .05
          .425 .075 .075 .425]
motif_frac_obs=.7
motif_recur_range=1:4

@info "Constructing synthetic sample set..."
obs1, bg_scores1, hmm_truth1, spike_truth1 = synthetic_sample(no_obs,obsl,bhmm_vec,bhmm_dist,[struc_sig_1,tata_box],[(true,(struc_frac_obs,periodicity)),(false,(motif_frac_obs,motif_recur_range))])

e1 = "BMISpikeTest"

#JOB CONSTANTS
const ensemble_size = 750
const no_sources = 4
const source_min_bases = 3
const source_max_bases = 12
const source_length_range= source_min_bases:source_max_bases
const mixing_prior = .5
const models_to_permute = ensemble_size * 3
funcvec=full_perm_funcvec
push!(funcvec, BioMotifInference.permute_source)
push!(funcvec, BioMotifInference.permute_source)
args=[Vector{Tuple{Symbol,Any}}() for i in 1:length(funcvec)]
args[end-1]=[(:weight_shift_freq,0.),(:length_change_freq,1.),(:length_perm_range,1:1)]
args[end]=[(:weight_shift_freq,.1),(:length_change_freq,0.),(:weight_shift_dist,Uniform(.00001,.01))]

instruct = Permute_Instruct(funcvec, ones(length(funcvec))./length(funcvec),models_to_permute,100; args=args)

@info "Assembling source priors..."
prior_array= Vector{Matrix{Float64}}()
source_priors = BioMotifInference.assemble_source_priors(no_sources, prior_array)

@info "Assembling ensemble..."
isfile(e1*"/ens") ? (ens1 = deserialize(e1*"/ens")) :
    (ens1 = IPM_Ensemble(e1, ensemble_size, source_priors, (falses(0,0), mixing_prior), bg_scores1, obs1, source_length_range))

@info "Converging ensemble..."
logZ1 = BioMotifInference.converge_ensemble!(ens1, instruct, backup=(true,250), wk_disp=false, tuning_disp=false, ens_disp=true, conv_plot=false, src_disp=true, lh_disp=false, liwi_disp=false)


text:
using Gumbo # can not be loaded inside of function => use script
using DataFrames

## hardcoded number of constituents
nAss = 505

url = "http://en.wikipedia.org/wiki/List_of_S%26P_500_companies"
htmlContent = download(url) |>
open |>
readall |>
parsehtml

## tableHeaders = [kk.root[2][3][4][4][7][1][1][1][1][1].text,
##                 kk.root[2][3][4][4][7][1][1][2][1].text,
##                 kk.root[2][3][4][4][7][1][1][3][1][1].text,
##                 kk.root[2][3][4][4][7][1][1][4][1][1].text,
##                 kk.root[2][3][4][4][7][1][1][5][1].text,
##                 kk.root[2][3][4][4][7][1][1][6][1].text,
##                 kk.root[2][3][4][4][7][1][1][7][1].text,
##                 kk.root[2][3][4][4][7][1][1][8][1][1].text]

function accessTableEntry(htmlContent::HTMLDocument, ii::Int)
    # old wikipedia page settings
    # return htmlContent.root[2][3][4][4][7][1][ii+1][1][1][1].text
    tickSymb = htmlContent.root[2][3][5][4][7][1][ii+1][1][1][1].text
    sector = htmlContent.root[2][3][5][4][7][1][ii+1][4][1].text
    subindustry = htmlContent.root[2][3][5][4][7][1][ii+1][5][1].text
    return (tickSymb, sector, subindustry)
end

companies = UTF8String[]
sectors = UTF8String[]
subindustries = UTF8String[]
for ii=1:nAss
    try
        tickSymb, sector, subindustry = accessTableEntry(htmlContent, ii)
        push!(companies, tickSymb)
        push!(sectors, sector)
        push!(subindustries, subindustry)
    catch
    end
end

pathToStore = joinpath(Pkg.dir("EconDatasets"), "data",
                       "SP500TickerSymbols.csv")

pathToStoreConsituents = joinpath(Pkg.dir("EconDatasets"), "data",
                                  "SP500Industries.csv")

df = DataFrame(TickerSymbol = companies, Sector = sectors,
               SubIndustry = subindustries)

writecsv(pathToStore, companies)
writetable(pathToStoreConsituents, df)


text:
# Abstract type for lensing operators
abstract type LenseOp <: ImplicitOp{Basis,Spin,Pix} end

const FΦTuple = FieldTuple{<:Basis,<:NamedTuple{(:f,:ϕ)}}

#
# δfϕₛ_δfϕₜ is an operator which computes the lensing jacobian between time s
# and time t,
#
# [δfₛ/δfₜ  δfₛ/δϕ;
#  δϕ/δfₜ  δϕ/δϕ]
#
# The operator acts on a FieldTuple containing (f,ϕ).
#
# Individual lensing algorithms need to implement how to apply this operator,
# inverse, and/or transpose. Different algorithms may need the field at time s,
# time t, or both, hence both are stored here
#
# Note, the bottom row is trivially [0,1], but included to make the Jacobian
# square and easier to reason about
#
struct δfϕₛ_δfϕₜ{s, t, L<:LenseOp, Fₛ<:Field, Fₜ<:Field} <: ImplicitOp{Basis,Spin,Pix}
    L::L
    fₛ::Fₛ
    fₜ::Fₜ
    δfϕₛ_δfϕₜ{s,t}(l::L,fₛ::Fₛ,fₜ::Fₜ) where {s,t,L,Fₛ,Fₜ} = new{s,t,L,Fₛ,Fₜ}(l,fₛ,fₜ)
end
# convenience constructors which use f̃ to mean f_(t=1) and/or f to mean f_(t=0)
δf̃ϕ_δfϕₜ(L,f̃,fₜ,::Val{t}) where {t} = δfϕₛ_δfϕₜ{1.,t}(L,f̃,fₜ)
δf̃ϕ_δfϕₜ(L,f̃,fₜ,::Val{1}) = IdentityOp
δfϕ_δfϕₜ(L,f,fₜ,::Val{t}) where {t} = δfϕₛ_δfϕₜ{0.,t}(L,f,fₜ)
δfϕ_δfϕₜ(L,f,fₜ,::Val{0}) = IdentityOp
δfϕ_δf̃ϕ(L,f,f̃) = δfϕₛ_δfϕₜ{0,1}(L,f,f̃)
δf̃ϕ_δfϕ(L,f̃,f) = δfϕₛ_δfϕₜ{1,0}(L,f̃,f)
# inverse Jacobians are the same as switching time t and s
\(J::δfϕₛ_δfϕₜ{s,t}, f::Field) where {s,t} = δfϕₛ_δfϕₜ{t,s}(J.L,J.fₜ,J.fₛ) * f

# operator for [δ/δϕ L(ϕ)*f]
# (this just involves picking out one block of δf̃ϕ_δfϕ)
δLf_δϕ(f, ϕ, ::Type{L}=LenseFlow) where {L} = δLf_δϕ(f, L(ϕ))
δLf_δϕ(f, L::LenseOp) = FuncOp(
    op  = g -> (δf̃ϕ_δfϕ(L,f,f)  * FΦTuple(g,zero(getϕ(L))))[2],
    opᴴ = g -> (δf̃ϕ_δfϕ(L,f,f)' * FΦTuple(g,zero(getϕ(L))))[2]
)

# some syntactic sugar for making lensing operators that lense from time t1 to t2
struct →{t1,t2} end
→(t1::Real,t2::Real) = →{float(t1),float(t2)}()
getindex(L::LenseOp, ::→{t,t}) where {t} = 1
getindex(L::LenseOp, ::→{0,1}) = L
getindex(L::LenseOp, i::→)  = _getindex(L,i)
_getindex(L::LenseOp, ::→{t1,t2}) where {t1,t2} = error("Lensing from time $t1 to $t2 with $(typeof(L)) is not implemented.")

struct NoLensing <: LenseOp end
NoLensing(ϕ) = NoLensing()
cache(::NoLensing, ::Field) = NoLensing()
cache!(::NoLensing, ::Field) = NoLensing()
*(::NoLensing, f::Field) = f
\(::NoLensing, f::Field) = f
adjoint(L::NoLensing) = L
_getindex(L::NoLensing, i::→) = L
δfϕₛ_δfϕₜ{t₀,t₁}(L::NoLensing,::Any,::Any) where {t₀,t₁} = IdentityOp


text:
using GZip

"""
    make_ref(src, ref)
---
With plink file `src.{bed,bim,fam}`, this function:
1. imputes the few missing genotypes in this dataset
2. phasing is also done by `beagle.jar`
3. write results to `ref.{bed,bim,fam}`
"""
function make_ref(src, ref; window=300)
    title("Make imputation reference")
    item("Check parameters")
    if !isfile(src * ".bed")
        warning("$src.bed doesn't exist")
        return
    end
    pathto, target = splitdir(ref)
    if length(pathto) == 0
        warning("Creating reference dataset in the current directory")
    else
        isdir(pathto) || mkdir(pathto)
    end
    done()

    item("Create reference")
    tmp = mktempdir(".")
    run(`$plink --sheep --bfile $src --recode vcf-iid bgz --out $tmp/mid`)
    run(pipeline(`java -jar $beagle
                    window=$window
                    gt=$tmp/mid.vcf.gz
                    ne=$nsNe
                    out=$ref`,
                 `grep Window`))
    rm(tmp, recursive=true, force=true)
    done()

    item("Convert to plink format")
    run(`$plink --sheep --make-bed --vcf $ref.vcf.gz --out $ref`)
    done()
end

"""
    make_ref(ref, fmp, plk...)
---
Make a reference `ref`, using linkage map file `fmp`.
This function will
- exclude the extra SNP in `plk...`, but not in `lmp`,
- merge the result `plk...`
- impute the *few* (usually are) missing SNP and create a reference file.
"""
function make_ref(ref, fmp, plk...)
    title("Merging several platforms into one reference")
    lmp = read_map(fmp)
    sst = Set(lmp.snp)          # SNP set
    dir, _ = splitdir(ref)      # write all the temp files and results in the ref dir.
    isdir(dir) || mkpath(dir)
    open("$dir/list", "w") do io
        for p in plk
            println(io, "$p-s")
        end
    end

    # check if a SNP needs flip.  I just remove them (308 SNP)
    alsnp = Dict{String, Set{Char}}()
    for p in plk
        bim = read_bim("$p.bim")
        for (_, snp, _, _, a1, a2) in eachrow(bim)
            if haskey(alsnp, snp)
                push!(alsnp[snp], a1)
            else
                alsnp[snp] = Set(a1)
            end
            push!(alsnp[snp], a2)
        end
    end
    excl = Set{String}()
    for (snp, val) in alsnp
        if length(val) > 2
            push!(excl, snp)
        end
    end

    for p in plk
        open("$dir/keep.snp", "w") do io
            for line in eachline("$p.bim")
                snp = split(line)[2]
                if snp ∈ sst && snp ∉ excl
                    println(io, snp)
                end
            end
        end
        run(`$plink --sheep --bfile $p
		--extract $dir/keep.snp
		--make-bed
		--out $p-s`)
    end
    run(`$plink --sheep --merge-list $dir/list --out $dir/1`)
    run(`$plink --sheep --bfile $dir/1 --list-duplicate-vars suppress-first --out $dir/1`)
    run(`$plink --sheep --bfile $dir/1 --exclude $dir/1.dupvar --make-bed --out $ref`)
    #make_ref("$dir/2", ref)
end

"""
    function impute(cur, ref, out; nthread=7)
---
Impute `cur`.{bed,bim,fam} with reference/training set `ref`.{bed,bim,fam}, and put
the result to `out`.{bed,bim,fam}.

**Note**:

- SNP only in `cur` will be removed.
- beagle was assigned 7 threads by default to avoid some weird errors.
"""
function impute(cur, ref, out; window=300)
    title("Impute $cur.{bed,bim,fam} with training set $ref")
    item("Check parameters")
    if !isfile(cur * ".bed")
        warning("$cur.bed doesn't exist")
        return
    end
    if !isfile("$ref.vcf.gz")
        warning("Reference $ref.vcf.gz is not ready")
        return
    end
    pathto, target = splitdir(out)
    if(length(pathto) == 0)
        warning("Creating target in the current directory")
    else
        isdir(pathto) || mkdir(pathto)
    end
    done()

    item("Prepare $cur and remove extra SNP in $cur")
    tmp = mktempdir(".")
    refsnp = begin
        df = read_bim("$ref.bim")
        Set(df.snp)
    end
    cursnp = begin
        df = read_bim("$cur.bim")
        Set(df.snp)
    end
    exclude = setdiff(cursnp, refsnp)
    write("$tmp/exclude.snp", join(exclude, '\n'), '\n')
    run(`$plink
            --sheep
            --bfile $cur
            --exclude $tmp/exclude.snp
            --recode vcf-iid bgz
            --out $tmp/mid`)
    done()

    item("Impute $cur with $ref")
    run(pipeline(`java -jar $beagle
                    window=$window
                    ne=$nsNe
                    ref=$ref.vcf.gz
                    gt=$tmp/mid.vcf.gz
                    out=$out`,
                 `grep Window`))
    done()

    item("Convert result to plink format")
    run(`$plink --sheep --make-bed --vcf $out.vcf.gz --out $out`)
    done()

    rm(tmp, recursive=true, force=true)
end

"""
    function merge_vcfs(tgt, refs...)
---
Merge several VCF files of the same platform,
or imputed to the same platform,
to create one reference.

This is a very dirty function.  The refs must be of the same platform.
"""
function merge_vcfs(tgt, refs...)
    out = gzopen("$tgt.vcf.gz", "w")
    gz = Any[]
    header = String[]

    for r in refs
        z = gzopen(r)
        push!(gz, z)
    end

    for line in eachline(gz[1])
        if line[2] == '#'
            println(out, line)
        else
            print(out, line)
            break
        end
    end

    for i in 2:length(gz)
        for line in eachline(gz[i])
            if line[2] != '#'
                print(out, '\t', join(split(line)[10:end], '\t'))
                break
            end
        end
    end
    print(out, '\n')

    for line in eachline(gz[1])
        print(out, line)
        for i in 2:length(gz)
            l = readline(gz[i])
            print(out, '\t', join(split(l)[10:end], '\t'))
        end
        print(out, '\n')
    end

    for z in gz
        close(z)
    end
    close(out)

    run(`$plink --sheep --make-bed --vcf $tgt.vcf.gz --out $tgt`)
end


text:
###############################################################################
#
# Update the fluorescence in the leaf
#
###############################################################################
"""
    leaf_fluorescence!(
                fluo_set::FluoParaSet{FT},
                leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Compute fluorescence yield, Kr, Ks, and Kp for leaf, given
- `fluo_set` [`FluoParaSet`](@ref) type parameter set
- `leaf` [`Leaf`](@ref) struct
"""
function leaf_fluorescence!(
            fluo_set::FluoParaSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    @unpack Ag, Kd, Kf, Kp_max, maxPSII = leaf;
    @unpack Kr1, Kr2, Kr3 = fluo_set;

    # Actual effective ETR:
    leaf.Ja  = max(0, Ag / leaf.e2c);

    # Effective photochemical yield:
    if leaf.Ja <= 0
        _φ   = maxPSII;
    else
        _φ   = maxPSII*leaf.Ja/leaf.J_pot;
    end

    leaf.φ   = min(1/maxPSII, _φ);
    # degree of light saturation: 'x' (van der Tol e.Ap. 2014)
    x        = max(0,  1-leaf.φ/maxPSII);

    # Max PSII rate constant, x_α = exp(log(x)*Kr2);
    x_α      = x ^ Kr2;
    leaf.Kr  = Kr1 * (1 + Kr3)* x_α / (Kr3 + x_α);
    leaf.Kp  = max(0, leaf.φ*(Kf+Kd+leaf.Kr)/(1-leaf.φ));

    leaf.Fo  = Kf / (Kf+Kp_max+Kd                );
    leaf.Fo′ = Kf / (Kf+Kp_max+Kd+leaf.Kr+leaf.Ks);
    leaf.Fm  = Kf / (Kf       +Kd                );
    leaf.Fm′ = Kf / (Kf       +Kd+leaf.Kr+leaf.Ks);
    leaf.ϕs  = leaf.Fm′ * (1 - leaf.φ);

    #
    #
    # TODO should these change with Ks?
    #
    #
    # leaf.eta  = leaf.ϕs/leaf.Fo
    # don't need this anymore
    # better to use ϕs directly for SIF as Fo is not always fqe=0.01
    leaf.qQ  = 1 - (leaf.ϕs-leaf.Fo′) / (leaf.Fm-leaf.Fo′);
    leaf.qE  = 1 - (leaf.Fm-leaf.Fo′) / (leaf.Fm′-leaf.Fo);
    leaf.NPQ = leaf.Kr / (Kf + Kd);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Test
using ConvOptDL.Utils

@testset "Utils Test (batched_outer)" begin
    A = hcat([ones(4)*i for i in 1:5]...)
    B = hcat([ones(3)*i for i in 1:5]...)
    refC = ones(4, 3, 5)
    for i in 1:5
        refC[:, :, i] .*= i^2
    end
    @test batched_outer(A, B) == refC
end

@testset "Utils Test (batched_kronecker)" begin
    A = [
        1 2;
        3 4;
    ]
    B = ones(3, 2)
    ans_ = repeat(repeat(A, inner=(3, 2)), outer=(1, 1, 5))
    @test batched_kronecker(repeat(A, outer=(1, 1, 5)), repeat(B, outer=(1, 1, 5))) == ans_
end


text:
### Vector Algorithms

import Base: sum, min, max, minimum, maximum, countnz, any, all, sort,
                union, find, cumsum, diff, findmax, findmin, sortperm

export sortIndex, sortByKey, diff2, minidx, maxidx

# Reduction

for (op,fn) in ((:sum, :af_sum_all), (:product, :af_product_all),
                (:maximum, :af_max_all), (:minimum, af_min_all))

    @eval function ($op){T<:Real}(a::AFArray{T})
        real = Base.Ref{Cdouble}(0)
        imag = Base.Ref{Cdouble}(0)
        eval($(quot(fn)))(real, imag, a)
        real[]
    end

    @eval function ($op){T<:Complex}(a::AFArray{T})
        real = Base.Ref{Cdouble}(0)
        imag = Base.Ref{Cdouble}(0)
        eval($(quot(fn)))(real, imag, a)
        complex(real[], imag[])
    end

end

function countnz{T}(a::AFArray{T})
    real = Base.Ref{Cdouble}(0)
    imag = Base.Ref{Cdouble}(0)
    af_count_all(real, imag, a)
    Int(real[])
end

function countnz{T}(a::AFArray{T}, dim::Integer)
    dim = dim - 1
    out = new_ptr()
    af_count(out, a, dim)
    AFArray{Cuint}(out[])
end

for (op,fn) in ((:any, :af_any_true_all),(:all, :af_all_true_all))

    @eval function ($op){T}(a::AFArray{T})
        real = Base.Ref{Cdouble}(0)
        imag = Base.Ref{Cdouble}(0)
        eval($fn)(real, imag, a)
        Bool(real[])
    end

end

for (op,fn) in ((:any, :af_any_true), (:all, :af_all_true))

    @eval function ($op){T}(a::AFArray{T}, dim::Integer)
        dim = dim - 1
        out = new_ptr()
        eval($fn)(out, a, dim)
        AFArray{Bool}(out[])
    end

end

for (op, fn) in ((:sum, :af_sum), (:product, :af_product),
                (:maximum, :af_max), (:minimum, :af_min))

    @eval function ($op){T}(a::AFArray{T}, dim::Integer)
        dim = dim - 1
        out = new_ptr()
        eval($(quot(fn)))(out, a, dim)
        AFArray{T}(out[])
    end

end

# Sorting

function sort{T}(a::AFArray{T}, dim::Integer = 1; rev = false)
    if dim == 2
        error("ArrayFire doesn't support sorting along this dimension")
    end
    out = new_ptr()
    af_sort(out, a, Cuint(dim-1), !rev)
    AFArray{T}(out[])
end

function sortperm{T}(a::AFArray{T}, dim::Integer = 1; rev = false)
    if dim == 2
        error("ArrayFire doesn't support sorting along this dimension")
    end
    out = new_ptr()
    indices = new_ptr()
    af_sort_index(out, indices, a, Cuint(dim-1), !rev)
    AFArray{Int32}(indices[]) + 1
end

function sortByKey{S,T}(keys::AFArray{S}, values::AFArray{T}, dim::Integer = 1; rev = false)
    if dim == 2
        error("ArrayFire doesn't support sorting along this dimension")
    end
    out_keys = new_ptr()
    out_values = new_ptr()
    af_sort_by_key(out_keys, out_values, keys, values, Cuint(dim - 1), !rev)
    AFArray{S}(out_keys[]), AFArray{T}(out_values[])
end

# Set Operations

function union{T,S}(a::AFArray{T}, b::AFArray{S}; is_unique = false)
    out = new_ptr()
    af_set_union(out, a, b, is_unique)
    AFArray{af_promote(T,S)}(out[])
end

function unique{T}(a::AFArray{T}; is_sorted = false)
    out = new_ptr()
    af_set_unique(out, a, is_sorted)
    AFArray{T}(out[])
end

function setdiff{T,S}(a::AFArray{T}, b::AFArray{S}; is_unique = false)
    out = new_ptr()
    af_set_intersect(out, a, b, is_unique)
    AFArray{af_promote(T,S)}(out[])
end

# Inclusive Scan Operations

function cumsum{T}(a::AFArray{T}, dim::Integer = 1)
    out = new_ptr()
    af_accum(out, a, Cint(dim-1))
    AFArray{T}(out[])
end

function find(a::AFArray)
    out = new_ptr()
    af_where(out, a)
    AFArray{backend_eltype(out[])}(out[]) + 1
end

# Numerical Differentiation

function diff{T}(a::AFArray{T}, dim::Integer = 1)
    out = new_ptr()
    af_diff1(out, a, Cint(dim-1))
    AFArray{T}(out[])
end

function diff2{T}(a::AFArray{T}, dim::Integer = 1)
    out = new_ptr()
    af_diff2(out, a, Cint(dim-1))
    AFArray{T}(out[])
end

# Max, Min with Indices

for (op, fn) in ((:findmax, :af_imax_all), (:findmin, :af_imin_all))

    @eval function ($op){T<:Real}(a::AFArray{T})
        real = Base.Ref{Cdouble}(0)
        imag = Base.Ref{Cdouble}(0)
        idx = Base.Ref{Cuint}(0)
        eval($fn)(real, imag, idx, a)
        real[], Int(idx[]) + 1
    end

    @eval function ($op){T<:Real}(a::AFArray{Complex{T}})
        real = Base.Ref{Cdouble}(0)
        imag = Base.Ref{Cdouble}(0)
        idx = Base.Ref{Cuint}(0)
        eval($fn)(real, imag, idx, a)
        complex(real[], imag[]), Int(idx[]) + 1
    end

end

for (op, fn) in ((:minidx, :af_imin), (:maxidx, :af_imax))

    @eval function ($op){T}(a::AFArray{T}, dim::Int)
        out = new_ptr()
        idx = new_ptr()
        eval($fn)(out, idx, a, dim-1)
        AFArray{T}(out[]),
        AFArray{backend_eltype(idx[])}(idx[]) + 1
    end

end


text:
using GtkBuilderAid
using Gtk
using Cairo
using Base.Test

# Force loading a few things into memory
# Need to do this to make the builder work in early versions
dummy = @GtkWindow("DUMMY")
dummy_box = @GtkBox(:v)
push!(dummy_box, @GtkLabel("dummy"))
dummy_button_box = @GtkButtonBox(:h)
push!(dummy_box, dummy_button_box)
push!(dummy_button_box, @GtkButton("dummy"))
destroy(dummy)

include("aidbuild.jl")
include("misc.jl")


text:
using DifferentialEquations, SparseArrays, BandedMatrices, Plots
using Optim,ApproxFun
include("cheby.jl")

const Λ = Float64(800.0)
const h = Float64(6.4)
const ρmax = Float64(0.03)
const S = Chebyshev(0.00..ρmax)
const ρ = points(S, gridnum)

flow(ρbar, Vp, Vpp, t) =
    (
        (3 / sqrt(Vp + exp(2 * t)) + 1 / sqrt(Vp + 2 * Vpp * ρbar + exp(2 * t))) *
        exp(5 * t)
    ) / (12 * pi^2)

flow2(ρbar, Vp, Vpp, t) =
    exp(6 * t) *
    (2 * Vp + 3 * Vpp * ρbar + 2 * exp(2 * t)) *
    (Vp + 2 * Vpp * ρbar + exp(2 * t))^-1 *
    (Vp + exp(2 * t))^-1 *
    (16 * pi^2)^-1

coeffa = Tnt * u0
coeffb = chderive(Float64(0.0), ρmax, coeffa, gridnum)
coeffc = chderive(Float64(0.0), ρmax, coeffb, gridnum)
Vp1 = Tn * coeffb
Vpp1 = Tn * coeffc
#display(plot(Vpp))

du .= flow.(ρ, Vp, Vpp, t)

function eqn(du, u, p, t)
    coeffa1 = Tnt * u
    coeffb1 = chderive(Float64(0.0), ρmax, coeffa1, gridnum)
    coeffc1 = chderive(Float64(0.0), ρmax, coeffb1, gridnum)
    Vp = Tn * coeffb1
    Vpp = Tn * coeffc1
    #display(plot(Vpp))
    du .= flow.(ρ, Vp, Vpp, t)
    println(t)
    #println(coeffa1[end-3:end])
end

function eqn2(du, u, p, t)
    V = Fun(S, ApproxFun.transform(S, u))
    du .= flow2.(ρ, V'.(ρ), V''.(ρ), t)
    println(t)
end

# coeffa = rand(gridnum)
# coeffb = rand(gridnum)
# coeffc = rand(gridnum)
# V1 = Fun(S, ApproxFun.transform(S, u0))
# du .= flow.(ρ, V'.(ρ), V''.(ρ), t)
# println(t)

λ * (x - Float64(0.013))^2
2 *λ * (0.0 - Float64(0.014))

const λ = Float64(16.0)
inifun1(x) = λ * (x - Float64(0.0121))^2
const m2_Λ = Float64(0.001)
const λ4 = Float64(0.1 / 6)
inifun2(x) = m2_Λ * x + λ4 * x^2

dt0=0.0
u0 = inifun1.(ρ)

prob = ODEProblem(eqn, u0, (Float64(0.0), dt0+Float64(log(1 / Λ))))
# alg = AutoTsit5(Rodas5())
# alg2 = AutoTsit5(TRBDF2())#-3.63
# alg3 = AutoTsit5(Rodas4P())
# alg4 = AutoTsit5(KenCarp3())#-3.6
# alg5 = AutoTsit5(QNDF1())
#
#
# alg6=  AutoVern9(ROS34PW2())
# alg7=  AutoVern9(TRBDF2())
#
# alg8=  AutoVern9(Rodas5())

alg7=  AutoTsit5(ROS34PW3())

using LinearAlgebra
BLAS.set_num_threads(1)

#choseable :
#Rodas5():7.77s
#Rodas4P():5.58s
#ROS34PW3():11.97s
#ROS34PW1a():12.88s
#ROS34PW1b():13.00s

@time sol = solve(
    prob,
    AutoVern7(ROS34PW3()),
    adaptive = true,
    reltol = 1e-6,
    abstol = 1e-6,
    save_everystep = true,
    dt=0.0000001,
    dtmax=0.001
)

c=0.35*10^6 *Λ^-3
Vout0 = Fun(S, ApproxFun.transform(S, u0))
Vout1 = Fun(S, ApproxFun.transform(S, (sol.u)[end]))
Vout2 = Fun(S, ApproxFun.transform(S, (sol.u)[end]-c*sqrt.(2*ρ)))

plot(ρ[1:end],u0[1:end],seriestype=:scatter)

plot(ρ[53:end],[(sol.u)[end][53:end]],seriestype=:scatter)

plot(ρ[123:end-20],[(sol.u)[end][123:end-20]],seriestype=:scatter)
plot!(ρ[123:end-20],[(sol.u)[end][123:end-20]],seriestype=:scatter)

plot(ρ,[(sol.u)[end]],seriestype=:scatter)

plot(Vout1,0.0,0.02)

sqrt(2 * (optimize(Vout1, 0.0, ρmax).minimizer) * Λ^2)

phi0=sqrt(2 * (optimize(Vout2, 0.0, ρmax).minimizer) * Λ^2)
rho0=optimize(Vout2, 0.0, ρmax).minimizer

sqrt((Vout1'(rho0)+2*rho0*Vout1''(rho0)) *Λ^2)

sqrt(Vout1'(rho0)*Λ^2)

sqrt(2*minpos*m2fun(minpos)*Λ^2)

testu=Tnt*(sol.u)[end]
Fun(testu)

max_stars_repo_path: Numerical/O(N) Solver.jl

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_path: Numerical/O(N) Solver.jl

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_path: Numerical/O(N) Solver.jl

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection