text:
using SpecialFunctions # for gamma function reasons
using LinearAlgebra    # for Gauss quadrature
using OffsetArrays     # for notational clarity
using SparseArrays     # for data structures

using GLMakie          # plotting
using LaTeXStrings     # plotting

"""
unimesh1D(xmin, xmax, K)
# Description
    Generates a uniform 1D mesh
# Arguments
    xmin: smallest value of array
    xmax: largest values of array
    K: number of elements in an array
# Return Values: VX, EtoV
    VX: vertex values | an Array of size K+1
    EtoV: element to node connectivity | a Matrix of size Kx2
# Example
xmin = -1
xmax =  1
K    =  4
VX, EtoV = unimesh1D(xmin, xmax, K)
"""
function unimesh1D(xmin, xmax, K)
    VX = collect(0:K) ./ K .* (xmax - xmin) .+ xmin
    EtoV = Int.(ones(K, 2))
    for i = 1:K
        EtoV[i,1] = Int(i)
        EtoV[i,2] = Int(i+1)
    end
    return VX, EtoV
end

# Mathy aliases
const Γ = gamma

# Coefficients in the Jacobi polynomial recurrence relations.
aᴾ(α, β, n) = 2/(2n+α+β) * √(n * (n+α+β) * (n+α) * (n+β) / (2n+α+β-1) / (2n+α+β+1))
bᴾ(α, β, n) = -(α^2 - β^2) / (2n+α+β) / (2n+α+β+2)

# code checked against the matlab code
"""
jacobi(x, α, β, n)
# Description
- Evaluates the jacobi polynomial at the point x
# Arguments
- `x`: point at which you will evaluate the jacobi polynomial
- `α`: first parameter for Jacobi polynomials
- `β`: second parameter for Jacobi polynomials
- `n` : order
# Return
-  `y`: the value of the of the Jacobi polynomial
"""
function jacobi(x, α, β, n::Int)
    Pᵅᵝ = n <= 1 ? OffsetArray(zeros(2), 0:1) : OffsetArray(zeros(n+1), 0:n)
    Pᵅᵝ[0] = √(2.0^-(α+β+1) * Γ(α+β+2) / Γ(α+1) / Γ(β+1))
    Pᵅᵝ[1] = Pᵅᵝ[0]/2 * √((α+β+3) / (α+1) / (β+1)) * ((α+β+2)*x + α - β)
    for n′ in 1:n-1
        Pᵅᵝ[n′+1] = ((x - bᴾ(α,β,n′)) * Pᵅᵝ[n′] - aᴾ(α,β,n′) * Pᵅᵝ[n′-1]) / aᴾ(α, β, n′+1)
    end
    return Pᵅᵝ[n]
end

"""
djacobi(x, α, β, n)
# Description
- Evaluates the derivative of the jacobi polynomial at the point x
# Arguments
- `x`: point at which you will evaluate the derivative of the jacobi polynomial
- `α`: first parameter for Jacobi polynomials
- `β`: second parameter for Jacobi polynomials
- `n` : order
# Return
-  `y`: the derivative of the of the Jacobi polynomial
"""
djacobi(x, α, β, n::Int) = √(n * (n+α+β+1)) * jacobi(x, α+1, β+1, n-1)

"""
vandermonde(x, α, β, N)
# Description
    Return vandermonde matrix of order N at the values x
    Allocates a little bit of memory
# Arguments
-   `x`: points at which to evaluate the Jacobi polynomials
-   `α`: first parameter for Jacobi polynomials
-   `β`: second parameter for Jacobi polynomials
-   `N`: maximum order of Jacobi polynomial to include
# Return Values
-   `v`: vandermonde matrix
# Example
    See LegendreTests.jl
"""
function vandermonde(x, α, β, N)
    # compute first two coefficients
    γ0 = 2^(α + β + 1) * factorial(α) * factorial(β) / ((α + β + 1) * factorial(α + β))
    γ1 = (α + 1) * (β + 1) / (α + β + 3) * γ0

    # create view to assign values
    v = zeros(length(x), N+1)
    v1 = view(v, :, 1)
    @. v1 = 1 / sqrt(γ0)

    # explicitly compute second coefficient
    if N == 0
        return v
    end

    v2 = view(v, :, 2)
    @. v2 = ( (α + β + 2) * x/2 + (α - β)/2) / sqrt(γ1)

    if N == 1
        return v
    end

    aʲ = 2 / (2 + α + β) * sqrt((α+1) * (β+1) / (α + β + 3))

    for i in 3:(N+1)
        # get views for ith, i-1th, and i-2th columns
        vi = view(v, :, i)
        vM1 = view(v, :, i-1)
        vM2 = view(v, :, i-2)

        # compute new a and b values
        h1 = 2 * (i-2) + α + β
        aⁱ = 2 / (h1 + 2) * sqrt((i-1) * (i-1 + α + β) * (i-1 + α) * (i-1 + β) / ((h1 + 1) * (h1 + 3)))
        bⁱ = - (α^2 - β^2) / (h1 * (h1 + 2))

        # compute coefficients for ith column
        @. vi = 1 / aⁱ * (-aʲ * vM2 + (x - bⁱ) * vM1)

        # save a coefficient for next iteration
        aʲ = aⁱ
    end

    return v
end

"""
dvandermonde(x, α, β, N)
# Description
    Return the gradient of the vandermonde matrix of order N at the values x
    Allocates a little bit of memory
# Arguments
-   `x`: points at which to evaluate the Jacobi polynomials
-   `α`: first parameter for Jacobi polynomials
-   `β`: second paramater for Jacobi polynomials
-   `N`: maximum order of Jacobi polynomial to include
# Return Values
-   `vr`: gradient of vandermonde matrix
# Example
    See LegendreTests.jl
"""
function dvandermonde(x, α, β, N)
    # create empty matrix (also handles first set of derivatives)
    vr = zeros(length(x), N+1)

    if N == 0
        return vr
    end

    # set values using vandermonde matrix
    v = vandermonde(x, α+1, β+1, N)
    for i in 1:N
        vi = view(v, :, i)
        vrP1 = view(vr, :, i+1)
        @. vrP1 = sqrt(i * (α + β + i+1)) * vi
    end

    return vr
end

"""
dmatrix(x, α, β, N)
# Description
    Return the differentiation matrix of order N at the values x
    Allocates too much memory
# Arguments
-   `x`: points at which to evaluate the Jacobi polynomials
-   `α`: first parameter for Jacobi polynomials
-   `β`: second paramater for Jacobi polynomials
-   `N`: maximum order of Jacobi polynomial to include
# Return Values
-   `D`: the differentiation matrix
# Example
    See LegendreTests.jl
"""
function dmatrix(x, α, β, N)
    # calculate vandermonde matrix and grad of vandermonde matrix
    vr = dvandermonde(x, α, β, N)
    v  =  vandermonde(x, α, β, N)

    # calculate values using D = vr * v^-1
    d = vr / v

    return d
end

"""
lift1D(V, y)
for computing fluxes
helps compute a surface integral of a quantity
note that the parentheses are necessary to prevent too much multiplcation
the E function takes the surface integrals are presents it
with respect to the full space inside an element
the entire operator represents how fluxes flow
into the interior of an element
"""
function lift1D(V)
    m,n = size(V)

    E = zeros(m , 2)
    E[1,1] = 1.0
    E[m,2] = 1.0

    return V * (transpose(V) * E)
end

"""
lift1D_v2(V, y)
for computing fluxes
nodal form
helps compute a surface integral of a quantity
note that the parantheses are necessary to prevent too much multiplcation
the E function takes the surface integrals are presents it
with respect to the full space inside an element
the entire operator represents how fluxes flow
into the interior of an element
"""
function lift1D_v2(V)
    m,n = size(V)

    E = zeros(m , 2)
    E[1,1] = 1.0
    E[m,2] = 1.0

    return E
end

"""
jacobiGQ(α, β, N)
# Description
    Guass Quadrature points and weights for the Jacobi Polynomial (α,β)
# Input
α, β: Jacobi polynomial descriptors
N:    order of quadrature points
# Return: x,w
x: quadrature points | array of size N+1
w: quadrature weights | array of size N+1
#Example
α = 0
β = 0
N = 4
x, w = jacobiGQ(α, β, N)
"""
function jacobiGQ(α, β, N)
    N == 0 && return [(α-β) / (α+β+2)], [2]

    # Form symmetric matrix from recurrence.
    dv = OffsetArray(zeros(N+1), 0:N)  # diagonal vector
    ev = OffsetArray(zeros(N+1), 0:N)  # sub/super-diagonal vector

    for n in 0:N
        dv[n] = bᴾ(α, β, n)
        ev[n] = aᴾ(α, β, n)
    end

    # Create full matrix combining the two.
    # Need to pass arrays that are not offset.
    J = SymTridiagonal(dv[0:N], ev[1:N])
    (α + β) ≈ 0 && (J[1, 1] = 0)

    # Compute quadrature points and weights by eigenvalue solve.
    x, V = eigen(J)
    w = @. V[1, :]^2 * 2^(α+β+1) / (α+β+1)
    @. w *= factorial(α) * factorial(β) / factorial(α+β)

    return x, w
end

"""
jacobiGL(α, β, N)
# Description
    Guass Labatto quadrature points for the Jacobi Polynomial (α,β)
    The quadrature weights are computed as well (but not returned)
# Arguments
- `α, β`: Jacobi polynomial descriptors
- `N`:    order of quadrature
# Return: x
- `x`: quadrature points  | array of size N+1
# Examples
```julia-repl
julia> x = jacobiGL(0, 0, 4)
5-element Array{Float64,1}:
 -1.0
 -0.6546536707079759
  4.440892098500626e-16
  0.6546536707079771
  1.0
```
"""
function jacobiGL(α, β, N)
    N == 0 && error("What are you doing? Go back to finite volume land.")
    N == 1 && return [-1, 1]

    x = zeros(N+1)
    x[1], x[N+1] = -1, 1

    x_GQ, _ = jacobiGQ(α+1, β+1, N-2)
    x[2:N] .= x_GQ

    return x
end

abstract type AbstractMesh end

"""
gridvalues1D(xmin, xmax, K)

# Description

    Generates physical gridpoints with each element

# Arguments

    VX: vertex values | an Array of size K+1

    EtoV: element to node connectivity | a Matrix of size Kx2

    r: LGL nodes in reference element | an array

# Return Values: x

    x: physical coordinates of solution

# Example (uses ../utils.jl as well)

xmin = 0
xmax = 2π
K = 4
# call functions
VX, EtoV = unimesh1D(xmin, xmax, K)
r = jacobiGL(0, 0, 4)
x = gridvalues1D(VX, EtoV, r)
# x[:,1] is the physical coordinates within the first element
# for plotting
f(x) = sin(x)
plot(x, f.(x))
# scatter(x,f.(x)) tends to work better
"""
function gridvalues1D(VX, EtoV, r)
    # get low and high edges
    va = view(EtoV, :, 1)
    vb = view(EtoV, :, 2)

    # compute physical coordinates of the grid points
    x = ones(length(r),1) * (VX[va]') .+ 0.5 .* (r .+ 1 ) * ((VX[vb] - VX[va])')
    return x
end

"""
facemask1D(r)

# Description

    creates face mask

# Arguments

    r: GL points

# Return Values: x

    fmask1: standard facemask
    fmask2: alternate form

# Example | ../utils.jl

r = jacobiGL(0, 0, 4)
fmask = fmask1D(r)

"""
function fmask1D(r)
    # check if index is left or right edge
    fm1 = @. abs(r+1) < eps(1.0);
    fm2 = @. abs(r-1) < eps(1.0);
    fmask1 = (fm1,fm2)

    # alternate form
    tmp = collect(1:length(r))
    fmask2  = [tmp[fm1]; tmp[fm2]]

    fmask = (fmask1, fmask2)
    return fmask
end

"""
edgevalues1D(fmask, x)

# Description

    calculates edge values

# Arguments

    fmask: face mask for GL edges

    x:  physical coordinates of solution on each element

# Return Values: x

    fx: face values of x

# Example | ../utils.jl

r = jacobiGL(0, 0, 4)
x = gridvalues1D(VX, EtoV, r)
fmask = fmask1D(r)[1]
fx = edgevalues1D(fmask,x)

# the locations of the edges in element 1 is fx[:, 1]

"""
function edgevalues1D(fmask, x)
    # compute x values at selected indices
    fx1 = x[fmask[1],:]
    fx2 = x[fmask[2],:]

    # return list of physical edge positions
    fx = [fx1; fx2]
    return fx
end

"""
normals1D(K)

# Description

    calculates face normals

# Arguments

    K: number of elements

# Return Values: normals

    normals: face normals along each grid

# Example

"""
function normals1D(K)
    normals  = ones(2,K)
    @. normals[1,:] *= -1
    return normals
end

"""
geometric_factors(x, Dʳ)

# Description

    computes the geometric factors for local mappings of 1D elements

# Arguments

    x: physical coordinates of solution for each element

    Dʳ:

# Return Values: rx, J

    rx: inverse jacobian

    J: jacobian (in 1D a scalar)

# Example

"""
function geometric_factors(x, Dʳ)
    J = Dʳ * x
    rx = 1 ./ J # for 1D
    return rx, J
end

"""
connect1D(EtoV)

# Description

    builds global connectivity arrays for 1D

# Arguments

    EtoV: element to node connectivity | a Matrix of size Kx2

# Return Values: EtoE, EtoF

    EtoE: element to element connectivity
    EtoF: element to face connectivity

# Example

"""
function connect1D(EtoV)
    nFaces = 2 # for 1d elements

    # Find number of elements and vertices
    K = size(EtoV,1)
    total_faces = nFaces * K
    Nv = K+1

    # list of local face to local vertex connections
    vn = [1, 2]

    # build global face to vertex array
    FtoV = Int.(spzeros(total_faces, Nv))
    let sk = 1
        for k = 1:K
            for faces = 1:nFaces
                FtoV[sk, EtoV[k, vn[faces]]] = 1;
                sk += 1
            end
        end
    end

    # build global face to face array
    FtoF = FtoV * (FtoV') - sparse(I, total_faces, total_faces)

    # find all face to face connections
    # check
    #(faces1, faces2) = findnz(FtoF)
    faces1, faces2 = findnz(FtoF .== 1)

    # convert global face number to element and face numbers
    element1 = @. floor(Int, (faces1 - 1) / nFaces) + 1
    face1    = @. Int( mod( (faces1 - 1),  nFaces) + 1)
    element2 = @. floor(Int, (faces2 - 1) / nFaces) + 1
    face2    = @. Int( mod( (faces2 - 1),  nFaces) + 1)

    # Rearrange into Nelement x Nfaces sized arrays
    ind = diag( LinearIndices(ones(K, nFaces))[element1,face1] ) # this line is a terrible idea.
    EtoE = collect(1:K) * ones(1, nFaces)
    EtoF = ones(K, 1) * (collect(1:nFaces)')
    EtoE[ind] = copy(element2);
    EtoF[ind] = face2;
    return EtoE, EtoF
end

"""
buildmaps1D(K, nGL, nFP, nFaces, fmask, EtoE, EtoF, x)
# Description

    connectivity matrices for element to elements and elements to face

# Arguments

-   `K`: number of elements
-   `nGL`: number of points within an element (polynomial degree + 1)
-   `nFP`: 1
-   `nFaces`: 2
-   `fmask`: an element by element mask to extract edge values
-   `EtoE`: element to element connectivity
-   `EtoF`: element to face connectivity
-   `x`: Guass lobatto points

# Return Values: vmapM, vmapP, vmapB, mapB, mapI, mapO, vmapI, vmapO

-   `vmapM`: vertex indices, (used for interior u values)
-   `vmapP`: vertex indices, (used for exterior u values)
-   `vmapB`: vertex indices, corresponding to boundaries
-   `mapB`: use to extract vmapB from vmapM
-   `mapI`: Index of left boundary condition
-   `mapO`: Index of right boundary condition

# Example | uses ../utils.jl

K = 3
n = 3; α = 0; β = 0; xmin = 0; xmax = 2π;
nGL = n + 1
nFP = 1
nFaces = 2

r = jacobiGL(α, β, n)

VX, EtoV = unimesh1D(xmin, xmax, K)
EtoE, EtoF = connect1D(EtoV)
x = gridvalues1D(VX, EtoV, r)
fx = edgevalues1D(r,x)

vmapM, vmapP, vmapB, mapB, mapI, mapO, vmapI, vmapO = buildmaps1D(K, nGL, nFP, nFaces, fmask, EtoE, EtoF, x)
"""
function buildmaps1D(K, nGL, nFP, nFaces, fmask, EtoE, EtoF, x)
    # number volume nodes consecutively
    nodeids = reshape(collect(1:(K*nGL)), nGL, K)
    vmapM = zeros(nFP, nFaces, K)
    vmapP = zeros(nFP, nFaces, K)
    # find index of face nodes wrt volume node ordering
    for k1 in 1:K
        for f1 in 1:nFaces
            vmapM[:, f1, k1] = nodeids[fmask[f1], k1]
        end
    end

    for k1 = 1:K
        for f1 = 1:nFaces
            # find neighbor
            k2 = Int.( EtoE[k1, f1])
            f2 = Int.( EtoF[k1, f1])

            # find volume node numbers of left and right nodes
            vidM = Int.( vmapM[:, f1, k1])
            vidP = Int.( vmapM[:, f2, k2])

            x1 = x[vidM]
            x2 = x[vidP]

            # compute distance matrix
            # need to figure out this part
            D = @. (x1 - x2)^2
            m = length(x1)
            for j = 1:m
                if D[j] < eps(1.0)*10^5
                    vmapP[j, f1, k1] = vidP[j]
                end
            end
        end
    end

    # reshape arrays
    vmapP = Int.( reshape(vmapP, length(vmapP)) )
    vmapM = Int.( reshape(vmapM, length(vmapM)) )

    # Create list of boundary nodes
    mapB = Int.( collect(1:length(vmapP))[vmapP .== vmapM] )
    vmapB = Int.( vmapM[mapB] )

    # inflow and outflow maps
    mapI = 1
    mapO = K * nFaces
    vmapI = 1
    vmapO = K*nGL
    return vmapM, vmapP, vmapB, mapB, mapI, mapO, vmapI, vmapO
end

"""
make_periodic1D!(vmapP, u)

# Description

    makes the grid periodic by modifying vmapP.
    Assumes that the first node is connected to the last.

# Arguments

    vmapP: exterior vertex map
    u: vertex vector

# Return Values: none

# Example

"""
function make_periodic1D!(vmapP, u)
    vmapP[1] = length(u)
    vmapP[end] = 1

    return nothing
end

struct Mesh{T,S,U,W} <: AbstractMesh
    # inputs
    K::S
    n::S

    # face stuff
    nFP::S
    nFaces::S

    # GL points
    r::U
    x::T

    # vertex maps
    vmapM::W
    vmapP::W
    vmapB::W
    mapB::W

    # inflow/outflow maps
    mapI::S
    mapO::S
    vmapI::S
    vmapO::S

    # structures for computation
    D::T
    M::T
    Mi::T
    lift::T
    rx::T
    normals::T
    fscale::T
end

"""
mesh(K, n, xmin, xmax)

# Description

    outer_constructor mesh struct

# Arguments

- `K`: number of elements
- `n`: polynomial order
- `xmin`: lower bound
- `xmax`: upper bound

# Return Values: x
    return grid values

"""
function Mesh(KK, nn, xmin, xmax; periodic = false)
    # initialize parameters
    K = KK
    α = 0; β = 0;
    n = nn

    # number of vertices
    nGL = n+1
    nFP = 1
    nFaces = 2

    # compute Gauss Lobatto grid
    r = jacobiGL(α, β, n)

    # build grid
    VX, EtoV = unimesh1D(xmin, xmax, K)

    # build coordinates of all the nodes
    x = gridvalues1D(VX, EtoV, r)

    # build connectivity matrix
    EtoE, EtoF = connect1D(EtoV)

    # build face masks
    fmask1,fmask2 = fmask1D(r)

    # build connectivity maps
    vmapM,vmapP,vmapB,mapB, mapI,mapO,vmapI,vmapO = buildmaps1D(K, nGL,nFP,nFaces, fmask1, EtoE,EtoF, x)

    # build differentiation matrix
    D = dmatrix(r, α, β, n)

    # build surface integral terms
    V = vandermonde(r, α, β, n)
    lift = lift1D(V)

    # build mass matrix and inverse of mass matrix
    Mi = V * V'
    M = inv(Mi)

    # calculate geometric factors
    rx,J = geometric_factors(x, D)

    # build surface normals
    normals = normals1D(K)

    # build inverse metric at the surface
    fscale = 1 ./ J[fmask2,:]

    # hack for periodicity
    if periodic
        vmapP[1] = vmapM[end]
        vmapP[end] = vmapM[1]
    end

    return Mesh{typeof(x),typeof(K),typeof(r),typeof(vmapP)}(K, n, nFP, nFaces, r, x, vmapM, vmapP, vmapB, mapB, mapI, mapO, vmapI, vmapO, D, M, Mi, lift, rx, normals, fscale)
end

function compute_volume_terms(data::AbstractArray, mesh::Mesh)
    q = mesh.D * data
    @. q *= mesh.rx
    return q
end

function compute_surface_terms(data::AbstractArray, 𝒢::Mesh)
    diffs = reshape( (data[𝒢.vmapM] - data[𝒢.vmapP]), (𝒢.nFP * 𝒢.nFaces, 𝒢.K ))
    # Include factor of 2 for the weak-strong form
    @. diffs *= 1.0 / 2.0
    # Compute Lift Operator
    lifted = - 𝒢.lift * (𝒢.fscale .* 𝒢.normals .* diffs)
    return lifted
end

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
module GenomicIndexes

import CodecBGZF

include("bgzfindex.jl")
include("bai.jl")
include("tabix.jl")
include("csi.jl")

end # module

ext: jl


text:
using CompScienceMeshes, BEAST

width, height, h = 1.0, 0.5, 0.05
G1 = meshrectangle(width, height, h)
G2 = CompScienceMeshes.rotate(G1, 0.5π * x̂)
G3 = CompScienceMeshes.rotate(G1, 1.0π * x̂)
G = weld(G1, G2, G3)

X = raviartthomas(G)

κ = 1.0
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n × E) × n
T = Maxwell3D.singlelayer(wavenumber=κ)

@hilbertspace j; @hilbertspace k
efie = @discretise T[k,j]==e[k]  j∈X k∈X
u = gmres(efie)

Θ, Φ = range(0.0,stop=π,length=100), 0.0
ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
farfield = potential(MWFarField3D(κ*im), ffpoints, u, X)

using Plots
using LinearAlgebra
plot(Θ,norm.(farfield))


text:
using Documenter, Pilot

makedocs(;
    modules=[Pilot],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/yngtodd/Pilot.jl/blob/{commit}{path}#L{line}",
    sitename="Pilot.jl",
    authors="yngtodd",
#    assets=String[],
)

deploydocs(;
    repo="github.com/yngtodd/Pilot.jl",
)


text:
using StaticCompiler
using StaticTools
using LoopVectorization

@inline function mul!(C::MallocArray, A::MallocArray, B::MallocArray)
    @turbo for n ∈ indices((C,B), 2), m ∈ indices((C,A), 1)
        Cmn = zero(eltype(C))
        for k ∈ indices((A,B), (2,1))
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
    return C
end

function loopvec_matrix(argc::Int, argv::Ptr{Ptr{UInt8}})
    argc == 3 || return printf(stderrp(), c"Incorrect number of command-line arguments\n")
    rows = parse(Int64, argv, 2)            # First command-line argument
    cols = parse(Int64, argv, 3)            # Second command-line argument

    # LHS
    A = MallocArray{Float64}(undef, rows, cols)
    @turbo for i ∈ axes(A, 1)
        for j ∈ axes(A, 2)
           A[i,j] = i*j
        end
    end

    # RHS
    B = MallocArray{Float64}(undef, cols, rows)
    @turbo for i ∈ axes(B, 1)
        for j ∈ axes(B, 2)
           B[i,j] = i*j
        end
    end

    # # Matrix multiplication
    C = MallocArray{Float64}(undef, cols, cols)
    mul!(C, B, A)

    # Print to stdout
    printf(C)
    # Also print to file
    fp = fopen(c"table.tsv",c"w")
    printf(fp, C)
    fclose(fp)
    # Clean up matrices
    free(A)
    free(B)
    free(C)
end

# Attempt to compile
path = compile_executable(loopvec_matrix, (Int64, Ptr{Ptr{UInt8}}), "./")


text:
module TestConstraint

using JuMP
using LinearAlgebra
using Test

include(joinpath(@__DIR__, "utilities.jl"))
include(joinpath(@__DIR__, "JuMPExtension.jl"))

function _test_constraint_name_util(constraint, name, F::Type, S::Type)
    @test name == @inferred JuMP.name(constraint)
    model = constraint.model
    @test constraint.index == JuMP.constraint_by_name(model, name).index
    if !(model isa JuMPExtension.MyModel)
        @test constraint.index ==
              JuMP.constraint_by_name(model, name, F, S).index
    end
end

function test_SingleVariable_constraints(ModelType, ::Any)
    m = ModelType()
    @variable(m, x)

    # x <= 10.0 doesn't translate to a SingleVariable constraint because
    # the LHS is first subtracted to form x - 10.0 <= 0.
    @constraint(m, cref, x in MOI.LessThan(10.0))
    _test_constraint_name_util(
        cref,
        "cref",
        JuMP.VariableRef,
        MOI.LessThan{Float64},
    )
    c = JuMP.constraint_object(cref)
    @test c.func == x
    @test c.set == MOI.LessThan(10.0)

    @variable(m, y[1:2])
    @constraint(m, cref2[i = 1:2], y[i] in MOI.LessThan(float(i)))
    _test_constraint_name_util(
        cref2[1],
        "cref2[1]",
        JuMP.VariableRef,
        MOI.LessThan{Float64},
    )
    c = JuMP.constraint_object(cref2[1])
    @test c.func == y[1]
    @test c.set == MOI.LessThan(1.0)
end

function test_VectorOfVariables_constraints(ModelType, ::Any)
    m = ModelType()
    @variable(m, x[1:2])

    cref = @constraint(m, x in MOI.Zeros(2))
    c = JuMP.constraint_object(cref)
    @test c.func == x
    @test c.set == MOI.Zeros(2)

    cref = @constraint(m, [x[2], x[1]] in MOI.Zeros(2))
    c = JuMP.constraint_object(cref)
    @test c.func == [x[2], x[1]]
    @test c.set == MOI.Zeros(2)

    @test_throws DimensionMismatch @constraint(m, [x[2], x[1]] in MOI.Zeros(3))
end

function test_AffExpr_scalar_constraints(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)

    cref = @constraint(model, 2x <= 10)
    @test "" == @inferred JuMP.name(cref)
    JuMP.set_name(cref, "c")
    _test_constraint_name_util(cref, "c", JuMP.AffExpr, MOI.LessThan{Float64})

    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, 2x)
    @test c.set == MOI.LessThan(10.0)

    cref = @constraint(model, 3x + 1 ≥ 10)
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, 3x)
    @test c.set == MOI.GreaterThan(9.0)

    cref = @constraint(model, 1 == -x)
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, 1.0x)
    @test c.set == MOI.EqualTo(-1.0)

    cref = @constraint(model, 2 == 1)
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, zero(JuMP.AffExpr))
    @test c.set == MOI.EqualTo(-1.0)
end

function test_AffExpr_vectorized_constraints(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)
    err = ErrorException(
        "In `@constraint(model, [x, 2x] == [1 - x, 3])`: Unexpected vector in scalar constraint. Did you mean to use the dot comparison operators like .==, .<=, and .>= instead?",
    )
    @test_throws_strip err @constraint(model, [x, 2x] == [1 - x, 3])
    @test_macro_throws ErrorException begin
        @constraint(model, [x == 1 - x, 2x == 3])
    end
    cref = @constraint(model, [x, 2x] .== [1 - x, 3])
    c = JuMP.constraint_object.(cref)
    @test JuMP.isequal_canonical(c[1].func, 2.0x)
    @test c[1].set == MOI.EqualTo(1.0)
    @test JuMP.isequal_canonical(c[2].func, 2.0x)
    @test c[2].set == MOI.EqualTo(3.0)
end

function test_AffExpr_vector_constraints(ModelType, ::Any)
    model = ModelType()
    cref = @constraint(model, [1, 2] in MOI.Zeros(2))
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func[1], zero(JuMP.AffExpr) + 1)
    @test JuMP.isequal_canonical(c.func[2], zero(JuMP.AffExpr) + 2)
    @test c.set == MOI.Zeros(2)
    @test c.shape isa JuMP.VectorShape
    @test_throws DimensionMismatch @constraint(model, [1, 2] in MOI.Zeros(3))
end

function test_delete_constraints(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)
    constraint_ref = @constraint(model, 2x <= 1)
    @test JuMP.is_valid(model, constraint_ref)
    JuMP.delete(model, constraint_ref)
    @test !JuMP.is_valid(model, constraint_ref)
    second_model = ModelType()
    @test_throws Exception JuMP.delete(second_model, constraint_ref)
end

function test_batch_delete_constraints(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:9])
    cons = [@constraint(model, sum(x[1:2:9]) <= 3)]
    push!(cons, @constraint(model, sum(x[2:2:8]) <= 2))
    push!(cons, @constraint(model, sum(x[1:3:9]) <= 1))
    @test all(JuMP.is_valid.(model, cons))
    JuMP.delete(model, cons[[1, 3]])
    @test all((!JuMP.is_valid).(model, cons[[1, 3]]))
    @test JuMP.is_valid(model, cons[2])
    second_model = ModelType()
    @test_throws Exception JuMP.delete(second_model, cons[[1, 3]])
    @test_throws Exception JuMP.delete(second_model, [cons[2]])
end

function test_two_sided_constraints(ModelType, ::Any)
    m = ModelType()
    @variable(m, x)
    @variable(m, y)

    @constraint(m, cref, 1.0 <= x + y + 1.0 <= 2.0)
    _test_constraint_name_util(
        cref,
        "cref",
        JuMP.AffExpr,
        MOI.Interval{Float64},
    )

    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, x + y)
    @test c.set == MOI.Interval(0.0, 1.0)

    cref = @constraint(m, 2x - y + 2.0 ∈ MOI.Interval(-1.0, 1.0))
    @test "" == @inferred JuMP.name(cref)

    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, 2x - y)
    @test c.set == MOI.Interval(-3.0, -1.0)
end

function test_broadcasted_constraint_eq(ModelType, ::Any)
    m = ModelType()
    @variable(m, x[1:2])

    A = [1.0 2.0; 3.0 4.0]
    b = [4.0, 5.0]

    cref = @constraint(m, A * x .== b)
    @test (2,) == @inferred size(cref)

    c1 = JuMP.constraint_object(cref[1])
    @test JuMP.isequal_canonical(c1.func, x[1] + 2x[2])
    @test c1.set == MOI.EqualTo(4.0)
    c2 = JuMP.constraint_object(cref[2])
    @test JuMP.isequal_canonical(c2.func, 3x[1] + 4x[2])
    @test c2.set == MOI.EqualTo(5.0)
end

function test_broadcasted_constraint_leq(ModelType, ::Any)
    m = ModelType()
    @variable(m, x[1:2, 1:2])

    UB = [1.0 2.0; 3.0 4.0]

    cref = @constraint(m, x .+ 1 .<= UB)
    @test (2, 2) == @inferred size(cref)
    for i in 1:2
        for j in 1:2
            c = JuMP.constraint_object(cref[i, j])
            @test JuMP.isequal_canonical(c.func, x[i, j] + 0)
            @test c.set == MOI.LessThan(UB[i, j] - 1)
        end
    end
end

function test_broadcasted_two_sided_constraint(ModelType, ::Any)
    m = ModelType()
    @variable(m, x[1:2])
    @variable(m, y[1:2])
    l = [1.0, 2.0]
    u = [3.0, 4.0]

    cref = @constraint(m, l .<= x + y .+ 1 .<= u)
    @test (2,) == @inferred size(cref)

    for i in 1:2
        c = JuMP.constraint_object(cref[i])
        @test JuMP.isequal_canonical(c.func, x[i] + y[i])
        @test c.set == MOI.Interval(l[i] - 1, u[i] - 1)
    end
end

function test_broadcasted_constraint_with_indices(ModelType, ::Any)
    m = ModelType()
    @variable m x[1:2]
    @constraint m cref1[i = 2:4] x .== [i, i + 1]
    ConstraintRefType = eltype(cref1[2])
    @test cref1 isa JuMP.Containers.DenseAxisArray{Vector{ConstraintRefType}}
    @constraint m cref2[i = 1:3, j = 1:4] x .≤ [i + j, i - j]
    ConstraintRefType = eltype(cref2[1])
    @test cref2 isa Matrix{Vector{ConstraintRefType}}
    @variable m y[1:2, 1:2]
    @constraint m cref3[i = 1:2] y[i, :] .== 1
    ConstraintRefType = eltype(cref3[1])
    @test cref3 isa Vector{Vector{ConstraintRefType}}
end

function test_quadexpr_constraints(ModelType, ::Any)
    model = ModelType()
    @variable(model, x)
    @variable(model, y)

    cref = @constraint(model, x^2 + x <= 1)
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, x^2 + x)
    @test c.set == MOI.LessThan(1.0)

    cref = @constraint(model, y * x - 1.0 == 0.0)
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func, x * y)
    @test c.set == MOI.EqualTo(1.0)

    cref = @constraint(
        model,
        [2x - 4x * y + 3x^2 - 1, -3y + 2x * y - 2x^2 + 1] in SecondOrderCone()
    )
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func[1], -1 + 3x^2 - 4x * y + 2x)
    @test JuMP.isequal_canonical(c.func[2], 1 - 2x^2 + 2x * y - 3y)
    @test c.set == MOI.SecondOrderCone(2)
end

function test_syntax_error_constraint(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:2])
    err = ErrorException(
        "In `@constraint(model, [3, x] in SecondOrderCone())`: unable to " *
        "add the constraint because we don't recognize $([3, x]) as a " *
        "valid JuMP function.",
    )
    @test_throws_strip err @constraint(model, [3, x] in SecondOrderCone())
end

function test_indicator_constraint(ModelType, ::Any)
    model = ModelType()
    @variable(model, a, Bin)
    @variable(model, b, Bin)
    @variable(model, x)
    @variable(model, y)
    for cref in [
        @constraint(model, a => {x + 2y <= 1}),
        @constraint(model, a ⇒ {x + 2y ≤ 1})
    ]
        c = JuMP.constraint_object(cref)
        @test c.func == [a, x + 2y]
        @test c.set == MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(1.0))
    end
    for cref in [
        @constraint(model, !b => {2x + y <= 1})
        @constraint(model, ¬b ⇒ {2x + y ≤ 1})
        # This returns a vector of constraints that is concatenated.
        @constraint(model, ![b, b] .=> {[2x + y, 2x + y] .≤ 1})
    ]
        c = JuMP.constraint_object(cref)
        @test c.func == [b, 2x + y]
        @test c.set == MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.LessThan(1.0))
    end
    err = ErrorException(
        "In `@constraint(model, !(a, b) => {x <= 1})`: Invalid binary variable expression `!(a, b)` for indicator constraint.",
    )
    @test_macro_throws err @constraint(model, !(a, b) => {x <= 1})
    err = ErrorException(
        "In `@constraint(model, a => x)`: Invalid right-hand side `x` of indicator constraint. Expected constraint surrounded by `{` and `}`.",
    )
    @test_macro_throws err @constraint(model, a => x)
    err = ErrorException(
        "In `@constraint(model, a => x <= 1)`: Invalid right-hand side `x <= 1` of indicator constraint. Expected constraint surrounded by `{` and `}`.",
    )
    @test_macro_throws err @constraint(model, a => x <= 1)
    err = ErrorException(
        "In `@constraint(model, a => {x <= 1, x >= 0})`: Invalid right-hand side `{x <= 1, x >= 0}` of indicator constraint. Expected constraint surrounded by `{` and `}`.",
    )
    @test_macro_throws err @constraint(model, a => {x <= 1, x >= 0})
    err = ErrorException(
        "In `@constraint(model, [a, b] .=> {x + y <= 1})`: Inconsistent use of `.` in symbols to indicate vectorization.",
    )
    @test_macro_throws err @constraint(model, [a, b] .=> {x + y <= 1})
end

function test_SDP_constraint(ModelType, VariableRefType)
    m = ModelType()
    @variable(m, x)
    @variable(m, y)
    @variable(m, z)
    @variable(m, w)

    cref = @constraint(m, [x y; z w] in PSDCone())
    c = JuMP.constraint_object(cref)
    @test c.func == [x, z, y, w]
    @test c.set == MOI.PositiveSemidefiniteConeSquare(2)
    @test c.shape isa JuMP.SquareMatrixShape

    @constraint(m, sym_ref, Symmetric([x 1; 1 -y] - [1 x; x -2]) in PSDCone())
    _test_constraint_name_util(
        sym_ref,
        "sym_ref",
        Vector{AffExpr},
        MOI.PositiveSemidefiniteConeTriangle,
    )
    c = JuMP.constraint_object(sym_ref)
    @test JuMP.isequal_canonical(c.func[1], x - 1)
    @test JuMP.isequal_canonical(c.func[2], 1 - x)
    @test JuMP.isequal_canonical(c.func[3], 2 - y)
    @test c.set == MOI.PositiveSemidefiniteConeTriangle(2)
    @test c.shape isa JuMP.SymmetricMatrixShape

    @SDconstraint(m, cref, [x 1; 1 -y] ⪰ [1 x; x -2])
    _test_constraint_name_util(
        cref,
        "cref",
        Vector{AffExpr},
        MOI.PositiveSemidefiniteConeSquare,
    )
    c = JuMP.constraint_object(cref)
    @test JuMP.isequal_canonical(c.func[1], x - 1)
    @test JuMP.isequal_canonical(c.func[2], 1 - x)
    @test JuMP.isequal_canonical(c.func[3], 1 - x)
    @test JuMP.isequal_canonical(c.func[4], 2 - y)
    @test c.set == MOI.PositiveSemidefiniteConeSquare(2)
    @test c.shape isa JuMP.SquareMatrixShape

    @SDconstraint(m, iref[i = 1:2], 0 ⪯ [x+i x+y; x+y -y])
    for i in 1:2
        _test_constraint_name_util(
            iref[i],
            "iref[$i]",
            Vector{AffExpr},
            MOI.PositiveSemidefiniteConeSquare,
        )
        c = JuMP.constraint_object(iref[i])
        @test JuMP.isequal_canonical(c.func[1], x + i)
        @test JuMP.isequal_canonical(c.func[2], x + y)
        @test JuMP.isequal_canonical(c.func[3], x + y)
        @test JuMP.isequal_canonical(c.func[4], -y)
        @test c.set == MOI.PositiveSemidefiniteConeSquare(2)
        @test c.shape isa JuMP.SquareMatrixShape
    end

    @SDconstraint(m, con_d, 0 ⪯ Diagonal([x, y]))
    c = JuMP.constraint_object(con_d)
    @test c.func isa Vector{GenericAffExpr{Float64,VariableRefType}}
    @test JuMP.isequal_canonical(c.func[1], 1x)
    @test iszero(c.func[2])
    @test iszero(c.func[3])
    @test JuMP.isequal_canonical(c.func[4], 1y)
    @test c.set == MOI.PositiveSemidefiniteConeSquare(2)

    @SDconstraint(m, con_d_sym, 0 ⪯ Symmetric(Diagonal([x, y])))
    c = JuMP.constraint_object(con_d_sym)
    @test c.func isa Vector{GenericAffExpr{Float64,VariableRefType}}
    @test JuMP.isequal_canonical(c.func[1], 1x)
    @test iszero(c.func[2])
    @test JuMP.isequal_canonical(c.func[3], 1y)
    @test c.set == MOI.PositiveSemidefiniteConeTriangle(2)

    @SDconstraint(m, con_td, Tridiagonal([z], [x, y], [w]) ⪰ 0)
    c = JuMP.constraint_object(con_td)
    @test c.func isa Vector{GenericAffExpr{Float64,VariableRefType}}
    @test JuMP.isequal_canonical(c.func[1], 1x)
    @test JuMP.isequal_canonical(c.func[2], 1z)
    @test JuMP.isequal_canonical(c.func[3], 1w)
    @test JuMP.isequal_canonical(c.func[4], 1y)
    @test c.set == MOI.PositiveSemidefiniteConeSquare(2)

    @SDconstraint(m, con_td_sym, Symmetric(Tridiagonal([z], [x, y], [w])) ⪰ 0)
    c = JuMP.constraint_object(con_td_sym)
    @test c.func isa Vector{GenericAffExpr{Float64,VariableRefType}}
    @test JuMP.isequal_canonical(c.func[1], 1x)
    @test JuMP.isequal_canonical(c.func[2], 1w)
    @test JuMP.isequal_canonical(c.func[3], 1y)
    @test c.set == MOI.PositiveSemidefiniteConeTriangle(2)

    @SDconstraint(m, con_ut, UpperTriangular([x y; z w]) ⪰ 0)
    c = JuMP.constraint_object(con_ut)
    @test c.func isa Vector{GenericAffExpr{Float64,VariableRefType}}
    @test JuMP.isequal_canonical(c.func[1], 1x)
    @test iszero(c.func[2])
    @test JuMP.isequal_canonical(c.func[3], 1y)
    @test JuMP.isequal_canonical(c.func[4], 1w)
    @test c.set == MOI.PositiveSemidefiniteConeSquare(2)

    @SDconstraint(m, con_lt, 0 ⪯ LowerTriangular([x y; z w]))
    c = JuMP.constraint_object(con_lt)
    @test c.func isa Vector{GenericAffExpr{Float64,VariableRefType}}
    @test JuMP.isequal_canonical(c.func[1], 1x)
    @test JuMP.isequal_canonical(c.func[2], 1z)
    @test iszero(c.func[3])
    @test JuMP.isequal_canonical(c.func[4], 1w)
    @test c.set == MOI.PositiveSemidefiniteConeSquare(2)

    # Julia changed how it reports keyword arguments between 1.3 and 1.4!
    err = if VERSION < v"1.4"
        ErrorException(
            "function build_constraint does not accept keyword arguments",
        )
    else
        MethodError
    end
    @test_throws(
        err,
        @SDconstraint(m, [x 1; 1 -y] ⪰ [1 x; x -2], unknown_kw = 1)
    )
    # Invalid sense == in SDP constraint
    @test_macro_throws ErrorException @SDconstraint(
        m,
        [x 1; 1 -y] == [1 x; x -2]
    )
end

function _test_constraint_name_util(ModelType, VariableRefType)
    model = ModelType()
    @variable(model, x)
    @constraint(model, con, x^2 == 1)
    _test_constraint_name_util(
        con,
        "con",
        GenericQuadExpr{Float64,VariableRefType},
        MOI.EqualTo{Float64},
    )
    JuMP.set_name(con, "kon")
    @test JuMP.constraint_by_name(model, "con") isa Nothing
    _test_constraint_name_util(
        con,
        "kon",
        GenericQuadExpr{Float64,VariableRefType},
        MOI.EqualTo{Float64},
    )
    y = @constraint(model, kon, [x^2, x] in SecondOrderCone())
    err(name) = ErrorException("Multiple constraints have the name $name.")
    @test_throws err("kon") JuMP.constraint_by_name(model, "kon")
    JuMP.set_name(kon, "con")
    _test_constraint_name_util(
        con,
        "kon",
        GenericQuadExpr{Float64,VariableRefType},
        MOI.EqualTo{Float64},
    )
    _test_constraint_name_util(
        kon,
        "con",
        Vector{GenericQuadExpr{Float64,VariableRefType}},
        MOI.SecondOrderCone,
    )
    JuMP.set_name(con, "con")
    @test_throws err("con") JuMP.constraint_by_name(model, "con")
    @test JuMP.constraint_by_name(model, "kon") isa Nothing
    JuMP.set_name(kon, "kon")
    _test_constraint_name_util(
        con,
        "con",
        GenericQuadExpr{Float64,VariableRefType},
        MOI.EqualTo{Float64},
    )
    return _test_constraint_name_util(
        kon,
        "kon",
        Vector{GenericQuadExpr{Float64,VariableRefType}},
        MOI.SecondOrderCone,
    )
end

function test_PSD_constraint_errors(ModelType, ::Any)
    model = ModelType()
    @variable(model, X[1:2, 1:2])
    err = ErrorException(
        "In `@constraint(model, X in MOI.PositiveSemidefiniteConeSquare(2))`:" *
        " instead of `MathOptInterface.PositiveSemidefiniteConeSquare(2)`," *
        " use `JuMP.PSDCone()`.",
    )
    @test_throws_strip(
        err,
        @constraint(model, X in MOI.PositiveSemidefiniteConeSquare(2))
    )
    err = ErrorException(
        "In `@constraint(model, X in MOI.PositiveSemidefiniteConeTriangle(2))`:" *
        " instead of `MathOptInterface.PositiveSemidefiniteConeTriangle(2)`," *
        " use `JuMP.PSDCone()`.",
    )
    @test_throws_strip(
        err,
        @constraint(model, X in MOI.PositiveSemidefiniteConeTriangle(2))
    )
end

function test_matrix_constraint_errors(ModelType, VariableRefType)
    model = ModelType()
    @variable(model, X[1:2, 1:2])
    err = ErrorException(
        "In `@constraint(model, X in MOI.SecondOrderCone(4))`: unexpected " *
        "matrix in vector constraint. Do you need to flatten the matrix " *
        "into a vector using `vec()`?",
    )
    # Note: this should apply to any MOI.AbstractVectorSet. We just pick
    # SecondOrderCone for convenience.
    @test_throws_strip(err, @constraint(model, X in MOI.SecondOrderCone(4)))
end

function test_nonsensical_SDP_constraint(ModelType, ::Any)
    m = ModelType()
    @test_throws_strip(
        ErrorException(
            "In `@variable(m, unequal[1:5, 1:6], PSD)`: Symmetric variables must be square. Got size (5, 6).",
        ),
        @variable(m, unequal[1:5, 1:6], PSD)
    )
    # Some of these errors happen at compile time, so we can't use @test_throws
    @test_throws MethodError @variable(m, notone[1:5, 2:6], PSD)
    @test_throws MethodError @variable(m, oneD[1:5], PSD)
    @test_throws MethodError @variable(m, threeD[1:5, 1:5, 1:5], PSD)
    @test_throws MethodError @variable(m, psd[2] <= rand(2, 2), PSD)

    @test_throws_strip(
        ErrorException(
            "In `@variable(m, -(ones(3, 4)) <= foo[1:4, 1:4] <= ones(4, 4), PSD)`: Non-symmetric bounds, integrality or starting values for symmetric variable.",
        ),
        @variable(m, -ones(3, 4) <= foo[1:4, 1:4] <= ones(4, 4), PSD)
    )
    @test_throws_strip(
        ErrorException(
            "In `@variable(m, -(ones(3, 4)) <= foo[1:4, 1:4] <= ones(4, 4), Symmetric)`: Non-symmetric bounds, integrality or starting values for symmetric variable.",
        ),
        @variable(m, -ones(3, 4) <= foo[1:4, 1:4] <= ones(4, 4), Symmetric)
    )
    @test_throws_strip(
        ErrorException(
            "In `@variable(m, -(ones(4, 4)) <= foo[1:4, 1:4] <= ones(4, 5), Symmetric)`: Non-symmetric bounds, integrality or starting values for symmetric variable.",
        ),
        @variable(m, -ones(4, 4) <= foo[1:4, 1:4] <= ones(4, 5), Symmetric)
    )
    @test_throws_strip(
        ErrorException(
            "In `@variable(m, -(rand(5, 5)) <= nonsymmetric[1:5, 1:5] <= rand(5, 5), Symmetric)`: Non-symmetric bounds, integrality or starting values for symmetric variable.",
        ),
        @variable(
            m,
            -rand(5, 5) <= nonsymmetric[1:5, 1:5] <= rand(5, 5),
            Symmetric
        )
    )
end

function test_sum_constraint(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:3, 1:3])
    @variable(model, y)
    C = [1 2 3; 4 5 6; 7 8 9]

    @test_expression sum(C[i, j] * x[i, j] for i in 1:2, j in 2:3)
    @test_expression sum(C[i, j] * x[i, j] for i = 1:3, j in 1:3 if i != j) - y
    @test JuMP.isequal_canonical(
        @expression(model, sum(C[i, j] * x[i, j] for i in 1:3, j in 1:i)),
        sum(C[i, j] * x[i, j] for i in 1:3 for j in 1:i),
    )
    @test_expression sum(C[i, j] * x[i, j] for i in 1:3 for j in 1:i)
    @test_expression sum(C[i, j] * x[i, j] for i = 1:3 if true for j in 1:i)
    @test_expression sum(
        C[i, j] * x[i, j] for i = 1:3 if true for j = 1:i if true
    )
    @test_expression sum(0 * x[i, 1] for i in 1:3)
    @test_expression sum(0 * x[i, 1] + y for i in 1:3)
    @test_expression sum(0 * x[i, 1] + y for i in 1:3 for j in 1:3)
end

function test_Model_all_constraints_scalar(::Any, ::Any)
    model = Model()
    @variable(model, x >= 0)
    @test 1 == @inferred num_constraints(
        model,
        VariableRef,
        MOI.GreaterThan{Float64},
    )
    ref =
        @inferred all_constraints(model, VariableRef, MOI.GreaterThan{Float64})
    @test ref == [LowerBoundRef(x)]
    @test 0 ==
          @inferred num_constraints(model, AffExpr, MOI.GreaterThan{Float64})
    aff_constraints = all_constraints(model, AffExpr, MOI.GreaterThan{Float64})
    @test isempty(aff_constraints)
    err = ErrorException(
        "`MathOptInterface.GreaterThan` is not a " *
        "concrete type. Did you miss a type parameter?",
    )
    @test_throws err num_constraints(model, AffExpr, MOI.GreaterThan)
    @test_throws err all_constraints(model, AffExpr, MOI.GreaterThan)
    err = ErrorException(
        "`GenericAffExpr` is not a concrete type. Did you miss a type parameter?",
    )
    @test_throws err try
        num_constraints(model, JuMP.GenericAffExpr, MOI.ZeroOne)
    catch e
        error(replace(e.msg, "JuMP." => ""))
    end
    @test_throws err try
        all_constraints(model, JuMP.GenericAffExpr, MOI.ZeroOne)
    catch e
        error(replace(e.msg, "JuMP." => ""))
    end
end

function test_Model_all_constraints_vector(::Any, ::Any)
    model = Model()
    @variable(model, x[1:2, 1:2], Symmetric)
    csdp = @constraint(model, x in PSDCone())
    csoc = @constraint(model, [x[1], 1] in SecondOrderCone())
    csos = @constraint(model, [x[2]^2, 1] in MOI.SOS1([1.0, 2.0]))
    @test_throws(
        DimensionMismatch,
        @constraint(model, [x[2]^2, 1] in MOI.SOS1([1.0, 2.0, 3.0]))
    )
    @test_throws(
        DimensionMismatch,
        @constraint(model, [x[2]^2, 1] in MOI.SOS2([1.0, 2.0, 3.0]))
    )
    @test 1 == @inferred num_constraints(
        model,
        Vector{VariableRef},
        MOI.PositiveSemidefiniteConeTriangle,
    )
    ref = all_constraints(
        model,
        Vector{VariableRef},
        MOI.PositiveSemidefiniteConeTriangle,
    )
    @test ref == [csdp]
    @test 1 ==
          @inferred num_constraints(model, Vector{AffExpr}, MOI.SecondOrderCone)
    ref = all_constraints(model, Vector{AffExpr}, MOI.SecondOrderCone)
    @test ref == [csoc]
    @test 1 ==
          @inferred num_constraints(model, Vector{QuadExpr}, MOI.SOS1{Float64})
    ref = all_constraints(model, Vector{QuadExpr}, MOI.SOS1{Float64})
    @test ref == [csos]
    @test 0 == @inferred num_constraints(
        model,
        Vector{AffExpr},
        MOI.PositiveSemidefiniteConeTriangle,
    )
    aff_constraints = all_constraints(
        model,
        Vector{AffExpr},
        MOI.PositiveSemidefiniteConeTriangle,
    )
    @test isempty(aff_constraints)
    err = ErrorException(
        "`$(GenericAffExpr{Float64})` is not a concrete type. Did you miss a " *
        "type parameter?",
    )
    @test_throws err try
        num_constraints(
            model,
            Vector{GenericAffExpr{Float64}},
            MOI.PositiveSemidefiniteConeTriangle,
        )
    catch e
        error(replace(e.msg, "JuMP." => ""))
    end
    @test_throws err try
        all_constraints(
            model,
            Vector{GenericAffExpr{Float64}},
            MOI.SecondOrderCone,
        )
    catch e
        error(replace(e.msg, "JuMP." => ""))
    end
    err = ErrorException(
        "`MathOptInterface.SOS1` is not a " *
        "concrete type. Did you miss a type parameter?",
    )
    @test_throws err all_constraints(
        model,
        Vector{GenericQuadExpr{Float64,VariableRef}},
        MOI.SOS1,
    )
end

function test_Model_list_of_constraint_types(::Any, ::Any)
    model = Model()
    @variable(model, x >= 0, Bin)
    @constraint(model, 2x <= 1)
    @constraint(model, [x, x] in SecondOrderCone())
    @constraint(model, [2x 1; 1 x] in PSDCone())
    @constraint(model, [x^2, x] in RotatedSecondOrderCone())
    constraint_types = @inferred list_of_constraint_types(model)
    @test Set(constraint_types) == Set([
        (VariableRef, MOI.ZeroOne),
        (VariableRef, MOI.GreaterThan{Float64}),
        (AffExpr, MOI.LessThan{Float64}),
        (Vector{VariableRef}, MOI.SecondOrderCone),
        (Vector{AffExpr}, MOI.PositiveSemidefiniteConeSquare),
        (Vector{QuadExpr}, MOI.RotatedSecondOrderCone),
    ])
end

function test_Model_dual_start(::Any, ::Any)
    model = Model()
    @variable(model, x)
    con = @constraint(model, 2x <= 1)
    @test dual_start_value(con) === nothing
    set_dual_start_value(con, 2)
    @test dual_start_value(con) == 2.0
    set_dual_start_value(con, nothing)
    @test dual_start_value(con) === nothing
end

function test_Model_dual_start_vector(::Any, ::Any)
    model = Model()
    @variable(model, x)
    con_vec = @constraint(model, [x, x] in SecondOrderCone())
    @test dual_start_value(con_vec) === nothing
    set_dual_start_value(con_vec, [1.0, 3.0])
    @test dual_start_value(con_vec) == [1.0, 3.0]
    set_dual_start_value(con_vec, nothing)
    @test dual_start_value(con_vec) === nothing
end

function test_Model_change_coefficient(::Any, ::Any)
    model = JuMP.Model()
    x = @variable(model)
    con_ref = @constraint(model, 2 * x == -1)
    @test JuMP.normalized_coefficient(con_ref, x) == 2.0
    JuMP.set_normalized_coefficient(con_ref, x, 1.0)
    @test JuMP.normalized_coefficient(con_ref, x) == 1.0
    JuMP.set_normalized_coefficient(con_ref, x, 3)  # Check type promotion.
    @test JuMP.normalized_coefficient(con_ref, x) == 3.0
    quad_con = @constraint(model, x^2 == 0)
    @test JuMP.normalized_coefficient(quad_con, x) == 0.0
    JuMP.set_normalized_coefficient(quad_con, x, 2)
    @test JuMP.normalized_coefficient(quad_con, x) == 2.0
    @test JuMP.isequal_canonical(
        JuMP.constraint_object(quad_con).func,
        x^2 + 2x,
    )
end

function test_Model_change_rhs(::Any, ::Any)
    model = JuMP.Model()
    x = @variable(model)
    con_ref = @constraint(model, 2 * x <= 1)
    @test JuMP.normalized_rhs(con_ref) == 1.0
    JuMP.set_normalized_rhs(con_ref, 2.0)
    @test JuMP.normalized_rhs(con_ref) == 2.0
    con_ref = @constraint(model, 2 * x - 1 == 1)
    @test JuMP.normalized_rhs(con_ref) == 2.0
    JuMP.set_normalized_rhs(con_ref, 3)
    @test JuMP.normalized_rhs(con_ref) == 3.0
    con_ref = @constraint(model, 0 <= 2 * x <= 1)
    @test_throws MethodError JuMP.set_normalized_rhs(con_ref, 3)
end

function test_Model_add_to_function_constant_scalar(::Any, ::Any)
    model = JuMP.Model()
    x = @variable(model)
    con_ref = @constraint(model, 2 <= 2 * x <= 3)
    con = constraint_object(con_ref)
    @test JuMP.isequal_canonical(JuMP.jump_function(con), 2x)
    @test JuMP.moi_set(con) == MOI.Interval(2.0, 3.0)
    JuMP.add_to_function_constant(con_ref, 1.0)
    con = constraint_object(con_ref)
    @test JuMP.isequal_canonical(JuMP.jump_function(con), 2x)
    @test JuMP.moi_set(con) == MOI.Interval(1.0, 2.0)
end

function test_Model_add_to_function_constant_vector(::Any, ::Any)
    model = JuMP.Model()
    x = @variable(model)
    con_ref = @constraint(model, [x + 1, x - 1] in MOI.Nonnegatives(2))
    con = constraint_object(con_ref)
    @test JuMP.isequal_canonical(JuMP.jump_function(con), [x + 1, x - 1])
    @test JuMP.moi_set(con) == MOI.Nonnegatives(2)
    JuMP.add_to_function_constant(con_ref, [2, 3])
    con = constraint_object(con_ref)
    @test JuMP.isequal_canonical(JuMP.jump_function(con), [x + 3, x + 2])
    @test JuMP.moi_set(con) == MOI.Nonnegatives(2)
end

function test_Model_value_constraint_var(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:2])
    @constraint(model, c1, x[1] + x[2] <= 3.0)
    @constraint(model, c2, x[1]^2 + x[2]^2 <= 3.0)
    @constraint(model, c3, [1.0, x[1], x[2]] in SecondOrderCone())

    vals = Dict(x[1] => 1.0, x[2] => 2.0)
    f = vidx -> vals[vidx]

    @test value(c1, f) === 3.0 # Affine expression
    @test value(c2, f) === 5.0 # Quadratic expression
    @test value(c3, f) == [1.0, 1.0, 2.0] # Vector expression
end

function _test_shadow_price_util(
    model_string,
    constraint_dual,
    constraint_shadow,
)
    model = Model()
    MOIU.loadfromstring!(backend(model), model_string)
    set_optimizer(
        model,
        () -> MOIU.MockOptimizer(
            MOIU.Model{Float64}();
            eval_objective_value = false,
            eval_variable_constraint_dual = false,
        ),
    )
    JuMP.optimize!(model)
    mock_optimizer = JuMP.backend(model).optimizer.model
    MOI.set(mock_optimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mock_optimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    JuMP.optimize!(model)
    for constraint_name in keys(constraint_dual)
        ci = MOI.get(backend(model), MOI.ConstraintIndex, constraint_name)
        constraint_ref = JuMP.ConstraintRef(model, ci, JuMP.ScalarShape())
        MOI.set(
            mock_optimizer,
            MOI.ConstraintDual(),
            JuMP.optimizer_index(constraint_ref),
            constraint_dual[constraint_name],
        )
        @test dual(constraint_ref) == constraint_dual[constraint_name]
        @test shadow_price(constraint_ref) == constraint_shadow[constraint_name]
    end
end

function test_Model_shadow_price(::Any, ::Any)
    _test_shadow_price_util(
        """
        variables: x, y
        minobjective: -1.0*x
        xub: x <= 2.0
        ylb: y >= 0.0
        c: x + y <= 1.0
        """,
        Dict("xub" => 0.0, "ylb" => 1.0, "c" => -1.0),
        Dict("xub" => 0.0, "ylb" => -1.0, "c" => -1.0),
    )

    _test_shadow_price_util(
        """
        variables: x, y
        maxobjective: 1.0*x
        xub: x <= 2.0
        ylb: y >= 0.0
        c: x + y <= 1.0
        """,
        Dict("xub" => 0.0, "ylb" => 1.0, "c" => -1.0),
        Dict("xub" => 0.0, "ylb" => 1.0, "c" => 1.0),
    )

    _test_shadow_price_util(
        """
        variables: x, y
        maxobjective: 1.0*x
        xub: x <= 2.0
        ylb: y >= 0.0
        """,
        Dict("xub" => -1.0, "ylb" => 0.0),
        Dict("xub" => 1.0, "ylb" => 0.0),
    )

    _test_shadow_price_util(
        """
        variables: x
        maxobjective: 1.0*x
        xeq: x == 2.0
        """,
        Dict("xeq" => -1.0),
        Dict("xeq" => 1.0),
    )

    return _test_shadow_price_util(
        """
        variables: x
        minobjective: 1.0*x
        xeq: x == 2.0
        """,
        Dict("xeq" => 1.0),
        Dict("xeq" => -1.0),
    )
end

function test_abstractarray_vector_constraint(ModelType, ::Any)
    model = ModelType()
    @variable(model, x[1:2, 1:2])
    c = @constraint(model, view(x, 1:4) in SOS1())
    obj = constraint_object(c)
    @test obj.func == x[1:4]
    @test obj.set == MOI.SOS1([1.0, 2.0, 3.0, 4.0])
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        f = getfield(@__MODULE__, name)
        @testset "$(name)" begin
            f(Model, VariableRef)
        end
        if !startswith("$(name)", "test_Model_")
            @testset "$(name)-JuMPExtension" begin
                f(JuMPExtension.MyModel, JuMPExtension.MyVariableRef)
            end
        end
    end
end

end

TestConstraint.runtests()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Julia example
const pi = 3.14

# create a parametric type
struct Point{T}
    x::T
    y::T
end
# create a point with integer coordinates
point1 = Point(1, 2)
println(point1)
# create a poit with real coordinates
point2 = Point{Float64}(1.0, 2.0)
println(point2)

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
function compare_to(obj::Enum, arg0::Enum)
    return jcall(obj, "compareTo", jint, (Enum,), arg0)
end

function describe_constable(obj::Enum)
    return jcall(obj, "describeConstable", Optional, ())
end

function equals(obj::Enum, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_angular_coordinates(obj::AEMAttitudeType, arg0::AbsoluteDate, arg1::Vector{jdouble}, arg2::jboolean, arg3::RotationOrder)
    return jcall(obj, "getAngularCoordinates", TimeStampedAngularCoordinates, (AbsoluteDate, Vector{jdouble}, jboolean, RotationOrder), arg0, arg1, arg2, arg3)
end

function get_attitude_data(obj::AEMAttitudeType, arg0::TimeStampedAngularCoordinates, arg1::jboolean, arg2::RotationOrder)
    return jcall(obj, "getAttitudeData", Vector{jdouble}, (TimeStampedAngularCoordinates, jboolean, RotationOrder), arg0, arg1, arg2)
end

function get_attitude_type(::Type{AEMAttitudeType}, arg0::JString)
    return jcall(AEMAttitudeType, "getAttitudeType", AEMAttitudeType, (JString,), arg0)
end

function get_declaring_class(obj::Enum)
    return jcall(obj, "getDeclaringClass", Class, ())
end

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

function hash_code(obj::Enum)
    return jcall(obj, "hashCode", jint, ())
end

function name(obj::Enum)
    return jcall(obj, "name", JString, ())
end

function ordinal(obj::Enum)
    return jcall(obj, "ordinal", jint, ())
end

function to_string(obj::Enum)
    return jcall(obj, "toString", JString, ())
end

function value_of(::Type{AEMAttitudeType}, arg0::JString)
    return jcall(AEMAttitudeType, "valueOf", AEMAttitudeType, (JString,), arg0)
end

function value_of(::Type{Enum}, arg0::Class, arg1::JString)
    return jcall(Enum, "valueOf", Enum, (Class, JString), arg0, arg1)
end

function values(::Type{AEMAttitudeType})
    return jcall(AEMAttitudeType, "values", Vector{AEMAttitudeType}, ())
end


text:
# mw_htest_i0_i1_pv - Construct H-test pvalue for the I(0) or Unit Root Null Model
# Syntax:  [output1,output2] = function_name(input1,input2,input3)
#
# Inputs:
#   x_t::Vector{Float64} - dct transformed data
#   param::Float64 - parameter for SD in H-test
#   itrend::Int64 - 1 for demean, 2 for detrend
#   ind_int::Int64 - 0 for I(0), 1 for I(1)
#   n_wsim::Int64 - number of simulated sd paths
#   rmatrix_string = name of matrix with standard random variables ... on disk ... w draws
#   nsim = number of simulations for p-value
# Outputs:
#    pv::FLoat64 - pvalue for h test
#
# Example:
#    Line 1 of example
#    Line 2 of example
#    Line 3 of example
#
# Other m-files required: none
# Subfunctions: none
# MAT-files required: none
#
# See also: OTHER_FUNCTION_NAME1,  OTHER_FUNCTION_NAME2

# Author: Pedro Chaim
# git repository: http://www.
# May 2004; Last revision: 12-May-2004
function mw_h_test_i0_i1_pv(x_t, param, itrend, ind_int, nwsim, rmatrix_string, nsim)
  q = size(x_t, 1)
  qq = q*(q+1)/2
  t = 100 # sample size for asymptotic approximations
  if itrend == 1
    psi = psi_compute(t, q)
  end
  # Construct Appropriate Covariance Matrices @
  cov_ha_i0, cov_ha_i1, cov_ha_i2 = mw_h_cov_ha_i0_i1_i2(psi, param, nwsim, rmatrix_string)
  if ind_int == 0
    cov_h0 = psi'psi
    cov_ha = cov_ha_i0
  elseif ind_int == 1
    cov_u = mw_rw_cov(t)
    cov_h0 = psi'cov_u*psi
    cov_ha = cov_ha_i1
  end
  # Construct Test Statistic @
  lr = mw_lrstat_vec(x_t, cov_h0, cov_ha_i0)
  # @ Compute Simulations for p-value @
  xsim = randn(q, nsim)
  a = chol(Hermitian(cov_h0, :L))
  xsim = a'xsim
  lr_sim = mw_lrstat_vec(xsim, cov_h0, cov_ha)
  pv = mean(lr_sim .> lr, 1)
return pv, lr#, lr, x_t, cov_h0, cov_ha_i0, cov_ha_i1, cov_ha_i2
end


text:
using Revise
using AdFem
using PyPlot

# forward computation
mmesh = Meh(joinpath(PDATA, "twoholes_large.stl"))
xy = gauss_nodes(mmesh)
κ = @. sin(xy[:,1]) * (1+xy[:,2]^2) + 1.0
f = 1e5 * @. xy[:,1] + xy[:,2]
K = compute_fem_laplace_matrix1(κ, mmesh)
F = compute_fem_source_term1(f, mmesh)
bdnode = bcnode(mmesh)
K, F = impose_Dirichlet_boundary_conditions(K, F, bdnode, zeros(length(bdnode)))
sol = K\F

# inverse modeling
nn_κ = squeeze(fc(xy, [20,20,20,1])) + 1
K = compute_fem_laplace_matrix1(nn_κ, mmesh)
F = compute_fem_source_term1(f, mmesh)
bdnode = bcnode(mmesh)
K, F = impose_Dirichlet_boundary_conditions(K, F, bdnode, zeros(length(bdnode)))
nn_sol = K\F
loss = sum((nn_sol - sol)^2)

sess = Session(); init(sess)
BFGS!(sess, loss)

nn_val = run(sess, nn_κ)
close("all")
figure(figsize=(18,5))
subplot(131)
visualize_scalar_on_gauss_points(κ, mmesh)
title("Reference")
subplot(132)
visualize_scalar_on_gauss_points(nn_val, mmesh)
title("DNN")
subplot(133)
visualize_scalar_on_gauss_points(abs.(κ-nn_val), mmesh)
title("Absolute Difference")
savefig("poisson_kappa.png")
close("all")
visualize_scalar_on_fem_points(sol, mmesh)
savefig("poisson_solution.png")


text:
module FourLEDs

using GtkReactive, Gtk.ShortNames, LibSerialPort, COBS

export main

baudrate = 115200
portfile = joinpath(DEPOT_PATH[1], "fourledsport")

function connect_slider(sp, id)
    sl = slider(0x00:0xff, value = 0x00)
    h = map(sl) do i
        encode(sp, [id, i])
    end
    return sl
end

function get_current_port_list(vb, msg, ok, win)
    hb = Box(:h)
    push!(vb, hb)
    msg = label(msg)
    ok = button(ok)
    push!(hb, msg)
    push!(hb, ok)
    ports = map(ok) do _
        LibSerialPort.get_port_list()
    end
    Gtk.showall(win)
    c = Condition()
    h = map(ok) do _
        notify(c)
    end
    wait(c)
    destroy(hb)
    return value(ports)
end

function get_port(vb, win)
    without = get_current_port_list(vb, "Is the LED strip disconnected?", "Yes", win)
    with = get_current_port_list(vb, "Now, connect the LED strip", "OK", win)
    port =  only(setdiff(with, without))
    open(portfile, "w") do o
        println(o, port)
    end
    return port
end

function get_sp(vb, win)
    if isfile(portfile)
        port = readline(portfile)
        try
            sp = LibSerialPort.open(port, baudrate)
            return sp
        catch ex
            rm(portfile)
        end
    end
    port = get_port(vb, win)
    sleep(1)
    return LibSerialPort.open(port, baudrate)
end

function main()
    win = Window("FourLEDs")
    vb = Box(:v)
    push!(win, vb)

    sp = get_sp(vb, win)

    for i in 0x00:0x03
        push!(vb, connect_slider(sp, i))
    end

    Gtk.showall(win)

    if !isinteractive()
        Gtk.waitforsignal(win, :destroy)
        close(sp)
    else
        signal_connect(win, :destroy) do widget
            close(sp)
        end
    end

    nothing

end

end


text:
using Test
using BenchmarkTools
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

@testset "Allocations" begin

  dir="./data/NAMD"
  atoms = readPDB("$dir/structure.pdb")

  options = Options(lastframe=1,nthreads=1,silent=true,seed=321,StableRNG=true)
  t_options = @benchmark Options(lastframe=1,seed=321,StableRNG=true,nthreads=1,silent=true) samples=1 evals=1
  @test t_options.allocs == 0

  protein = Selection(select(atoms,"protein"),nmols=1)
  t_selection1 = @benchmark Selection(select($atoms,"protein"),nmols=1) samples=1 evals=1
  @test t_selection1.allocs == 62082

  tmao = Selection(select(atoms,"resname TMAO"),natomspermol=14)
  t_selection2 = @benchmark Selection(select($atoms,"resname TMAO"),natomspermol=14) samples=1 evals=1
  @test t_selection2.allocs == 186144

  trajfile = "$dir/trajectory.dcd" # because of the interpolation of @benchmark
  traj = Trajectory(trajfile,protein,tmao)
  t_trajectory = @benchmark Trajectory($trajfile,$protein,$tmao) samples=1 evals=1
  @test t_trajectory.allocs == 839

  samples = CM.Samples(md=(traj.solvent.nmols-1)/2,random=options.n_random_samples)
  t_samples = @benchmark CM.Samples(md=($traj.solvent.nmols-1)/2,random=$options.n_random_samples) samples=1 evals=1
  @test t_samples.allocs == 0

  R = Result(traj,options)
  t_result = @benchmark Result($traj,$options) samples=1 evals=1
  @test t_result.allocs == 66

  framedata = CM.FrameData(traj,R)
  t_framedata = @benchmark CM.FrameData($traj,$R) samples=1 evals=1
  @test t_framedata.allocs == 206

  CM.nextframe!(traj)
  t_nextframe = @benchmark CM.nextframe!($traj) samples=1 evals=1
  @test t_nextframe.allocs == 34

  RNG = CM.init_random(options)
  t_RNG = @benchmark CM.init_random($options) samples=1 evals=1
  @test t_RNG.allocs == 2

  CM.mddf_frame!(1,framedata,options,RNG,R)
  t_mddf_frame = @benchmark CM.mddf_frame!(1,$framedata,$options,$RNG,$R) samples=1 evals=1
  @test t_mddf_frame.allocs == 0

end


text:
"""
abstractwca.jl

Common datatypes and functions for WCA and LWCA

References:
"""

abstract AbstractWCA <: DependentReferenceSystem
abstract AbstractOptimizedWCA <: AbstractWCA

ncomp(wca::AbstractWCA)::Int = ncomp(wca.trial)
numberdensity(wca::AbstractWCA)::Vector{Float64} = numberdensity(wca.trial)
totalnumberdensity(wca::AbstractWCA)::Float64 = totalnumberdensity(wca.trial)
composition(wca::AbstractWCA)::Vector{Float64} = composition(wca.trial)
temperature(wca::AbstractWCA)::Float64 = wca.temp
hsdiameter(wca::AbstractWCA)::Array{Float64,2} = hsdiameter(wca.trial)

function emptyradius(wca::AbstractOptimizedWCA)::Array{Float64,2}
  N::Int = ncomp(wca)
  σ::Matrix{Float64} = hsdiameter(wca)
  B::Matrix{Function} = blipfunction(wca)

  ret = Matrix{Float64}(N,N)
  for i in 1:N, j in 1:N
    i > j && continue

    opt = NLopt.Opt(:LN_BOBYQA, 1)
    NLopt.min_objective!(opt, (r, g) -> abs(B[i,j](r[1])))
    NLopt.lower_bounds!(opt, [0.])
    NLopt.upper_bounds!(opt, [σ[i,j]])
    NLopt.stopval!(opt, eps(Float64))

    (fmin, r₀, res) = NLopt.optimize(opt, [σ[i,j]])

    ret[i,j] = ret[j,i] = r₀[1]
  end
  ret
end

# Ref: N. E. Dubinin et al.: Thermochimica Acta, 518 (2011) 9-12.
function entropy(wca::AbstractOptimizedWCA)::Float64
  T::Float64 = temperature(wca)
  S_hs::Float64 = entropy(wca.trial)
  U₀_wca::Float64 = internal(wca)

  ΔS₀_wca = U₀_wca / T
  S = S_hs + ΔS₀_wca
end

# Perturbation-independent part of internal energy
# Ref: N. E. Dubinin et al.: Thermochimica Acta, 518 (2011) 9-12.
function internal(wca::AbstractOptimizedWCA)::Float64
  N::Int = ncomp(wca)
  T::Float64 = temperature(wca)
  ρ::Float64 = totalnumberdensity(wca)
  c::Vector{Float64} = composition(wca)
  hs::IndependentReferenceSystem = wca.trial
  σ::Array{Float64,2} = hsdiameter(hs)
  rmin::Array{Float64,2} = wca.rmin
  rcore::Array{Float64,2} = emptyradius(wca)
  u₀::Array{Function,2} = repulsivepotential(wca)
  g₀_wca::Array{Function,2} = paircorrelation(wca)

  U₀_wca::Float64 = 0

  for i in 1:N, j in 1:N
    i > j && continue

    n = i == j ? 1 : 2
    g₀s_wca::Function = spline(g₀_wca[i,j], rcore[i,j], rmin[i,j], 32)
    u₀s::Function = spline(u₀[i,j], rcore[i,j], rmin[i,j], 32)

    U₀_wca +=
      2π*ρ * n*c[i]*c[j] * ∫(r -> u₀s(r)*g₀s_wca(r)*r^2, rcore[i,j], rmin[i,j])
  end

  return U₀_wca
end


text:
#const ongoing_dir = "/Volumes/GoogleDrive/My Drive"
const ongoing_dir = "/home/beatriz/mainen.flipping.5ht@gmail.com/"
const Directory = joinpath(ongoing_dir,"Flipping/Datasets/Pharmacology/SwitchingData/rawdata")
const saving_dir = joinpath(ongoing_dir,"Flipping/Datasets/Pharmacology/SwitchingData/Results")

function files_dict(;dir = Directory)
    files = readdir(dir)
    filter!(t-> occursin("data_@switching2",t),files)
    return Dict("Dir" => dir,"Files" => files)
end

const Files_Dict = files_dict()

function get_rawfile(idx)
    dir = Files_Dict["Dir"]
    file = Files_Dict["Files"][idx]
    joinpath(dir,file)
end

function read_rawfile(file::String)
    MAT.matread(file)
end

function read_rawfile(idx::Int64)
    file_path = get_rawfile(idx)
    MAT.matread(file_path)
end


text:
using StatsModels, BlockDiagonals
"""
    olsranef(model::LinearMixedModel, method=:simultaneous)

Compute the group-level estimates using ordinary least squares.

This is somewhat similar to the conditional modes / BLUPs computed without
shrinkage.

There is no shrinkage of the group-level estimates with this approach, which
means singular estimates can be avoided. However, this also means the random
effects design matrix must not be singular.

Two methods are provided:
1. (default) OLS estimates computed for all strata (blocking variables)
   simultaneously with `method=simultaneous`. This pools the variance
   across estimates but does not shrink the estimates. Note that this method
   internal reparameterizes the random effects matrix `Z` to use effects coding
   and only use a single intercept shared across all grouping variables.
2. OLS estimates computed within each stratum with `method=stratum`. This
   method is equivalent for example to computing each subject-level and each
   item-level regression separately.

For fully balanced designs with a single blocking variable, these methods will
give the same results.

!!! warning
    If the design matrix for the random effects is rank deficient (e.g., through
    the use of `MixedModels.fulldummy` or missing cells in the data), then these
    methods will fail because no shrinkage/regularization is applied.
"""
function olsranef(model::LinearMixedModel{T}, method=:simultaneous) where {T}
    fixef_res = copy(response(model))
    # what's not explained by the fixed effects has to be explained by the RE
    X = model.X
    β = model.β # (X'X) \ (X'fixef_res)
    mul!(fixef_res, X, β, -one(T), one(T))

    return olsranef(model, fixef_res, Val(method))
end

function olsranef_org(model::LinearMixedModel{T}, fixef_res, ::Val{:stratum}) where {T}

    blups = Vector{Matrix{T}}()
    for trm in model.reterms
        re = []

        z = trm.z
        refs = unique(trm.refs)
        ngrps = length(refs)
        npreds = length(trm.cnames)

        j = 1
        for r in refs
            i = trm.refs .== r
            X = trm[i, j:(j+npreds-1)]
            b = (X'X) \ (X'fixef_res[i])
            push!(re, b)
            j += npreds
        end

        push!(blups, hcat(re...))
    end
    return blups, dummy_scalings(model.reterms)
end

function dummy_scalings(reterms)
     scalings = repeat([LinearAlgebra.I],length(reterms))
      scalings = vcat(scalings,1.) # add sigma scaling
      return scalings

end
function olsranef(model::LinearMixedModel{T}, fixef_res, ::Val{:simultaneous}) where {T}
    n_reterms = size(model.reterms)[1]

    code = Vector{Matrix{T}}(undef, n_reterms) # new contrast matrices
    mat = Vector{Matrix{T}}(undef, n_reterms) # new Z model matrices

    ### create new Z matrices with orthogonal contrasts
    for (i, trm) in enumerate(model.reterms)
        dim = size(trm.z)[1]
        cd = StatsModels.ContrastsMatrix(EffectsCoding(), trm.levels).matrix
        cd = kron(cd, I(dim))
        code[i] = cd
        mat[i] = trm * cd
    end

    # add in intercept
    # note that we pass two dims to ones to create a matrix
    pushfirst!(mat, ones(sum(first ∘ size, mat), 1))
    Z = hcat(mat...)
    # compute the BLUPs
    flatblups = Z'Z \ Z'fixef_res
    # get back to original coding
    flatblups = BlockDiagonal(code) * @view flatblups[2:end, :]

    blups = Vector{Matrix{T}}()
    offset = 1
    for trm in model.reterms
        chunksize = size(trm, 2)
        ngrps = length(trm.levels)
        npreds = length(trm.cnames)
        re = Matrix{T}(reshape(view(flatblups, offset:(offset+chunksize-1)),
                               npreds, ngrps))
        offset += chunksize
        push!(blups, re)
    end

    return blups, dummy_scalings(model.reterms)
end

"""
    residuals(model::LinearMixedModel, blups)

Compute the residuals of a mixed model using the specified group-level BLUPs/predictors.

This is useful for, e.g., comparing the residuals from a mixed-effects model with shrunken
group-level predictors against a non-shrunken classical OLS model fit within each group.
"""
function MixedModels.residuals(model::LinearMixedModel{T}, blups::Vector{<:AbstractMatrix{T}}) where T
    # XXX This is kinda type piracy, if it weren't developed by one of the MixedModels.jl devs....

    y = response(model) # we are now modifying the model

    ŷ = zeros(T, length(y))

    for (re, trm) in zip(blups, model.reterms)
        mul!(ŷ, trm, re, one(T), one(T))
    end

    mul!(ŷ, model.X, model.β, one(T), one(T))

    # TODO: do this inplace to avoid an allocation
    return y .- ŷ
end


text:
export ARCLDLt_HO_vs_Nwt

function ARCLDLt_HO_vs_Nwt(nlp 	   :: AbstractNLPModel,
              	   nlpstop :: NLPStopping;
				   corr_ho :: Symbol = :Shamanskii,
				   nwt_res_fact = 0.8,
				   λfact = 100.0,
				   kwargs...
               	   )

	T = eltype(nlp.meta.x0)

	return TRARC(nlp,
				  nlpstop;
				  TR = TrustRegion(T(10.0)),
				  c = Combi(hessian_dense, PDataLDLt{eltype(nlp.meta.x0)},  (x, y, z) -> solve_modelARCDiag_HO_vs_Nwt(x, y, z, λfact = λfact, nwt_res_fact = nwt_res_fact), preprocessLDLt, decreaseFact, Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
include("ulsp.jl")
Random.seed!(100003)

useCPLEX=0

periodos=10
ratioFijoVariable=250
mu=100
sigma=15
sinicial=0

periodos=parse(Int64,ARGS[1])
tipo=parse(Int64,ARGS[2])
precioFijo,precioStock,Demanda=generarInstancia(periodos,ratioFijoVariable,mu,sigma)
if tipo==1
  start=time()
  ulsp(periodos,precioFijo,precioStock,Demanda,sinicial,useCPLEX)
  println("tiempo modelo basico $periodos periodos: ",time()-start)
end
if tipo==2
  start=time()
  ulspAlt(periodos,precioFijo,precioStock,Demanda,sinicial,useCPLEX)
  println("tiempo modelo alternativo $periodos periodos: ",time()-start)
end
if tipo==3
  start=time()
  ulspAlt2(periodos,precioFijo,precioStock,Demanda,sinicial,useCPLEX)
  println("tiempo modelo alternativo2 $periodos periodos: ",time()-start)
end


text:
onethird(::Type{Float64}) = 0.3333333333333333
onethird(::Type{Double64}) = Double64((0.3333333333333333, 1.850371707708594e-17))

function Base.sqrt(x::Double64)
    x0 = sqrt(Hi(x))
    x1 = (x0 + x/x0) * 0.5
    x2 = (x1 + x/x1) * 0.5
    return x2
end

@inline function fastsqrt(x::Double64)
    x0 = sqrt(Hi(x))
    x1 = (x0 + x/x0) * 0.5
    return x1
end

function Base.cbrt(x::Double64)
    x0 = Double64(cbrt(Hi(x)))
    a =  fastsqrt( ((4*x/x0) - x0^2) / 3)
    x1 = (x0 + a)/2
    return x1
end

function fastcbrt(x::Double64)
    x0 = cbrt(Hi(x))
    x1 = (2*x0 + x/(x0*x0)) / 3.0
    x2 = (2*x1 + x/(x1*x1)) / 3.0
    return x2
end


text:
"""
    GridOperators

Provides local finite-difference
operator functions.
"""
module GridOperators

export ∇_z, Δ_z, adv_upwind
using ..Grids

"""
    ∇_z(f, grid::Grid)

Computes the local derivative of field ``f``:
``∇f``
"""
∇_z(f, grid::Grid) = (f[1] + f[3] - 2*f[2])*grid.dzi

"""
    Δ_z(f, grid::Grid)

Computes the local Laplacian of field ``f``:
``∇ • (∇f)``
"""
Δ_z(f, grid::Grid) = (f[3] - 2*f[2] + f[1])*grid.dzi2

"""
    Δ_z(f, grid::Grid, K)

Computes the local Laplacian of field ``f``
with a variable coefficient:
``∇ • (K ∇f)``
"""
function Δ_z(f, grid::Grid, K)
  K∇f_left  = (K[1]+K[2])/2*(f[2]-f[1])*grid.dzi
  K∇f_right = (K[2]+K[3])/2*(f[3]-f[2])*grid.dzi
  return (K∇f_right - K∇f_left)*grid.dzi
end

"""
    adv_upwind(ϕ, u, grid::Grid)

Local upwind advection operator ``u • ∇ϕ``. This
operator is stable but numerically diffusive.
"""
adv_upwind(ϕ, u, grid::Grid) = u[2] > 0 ? u[ 2 ]*(ϕ[2] - ϕ[1]) * grid.dzi :
                                          u[ 2 ]*(ϕ[3] - ϕ[2]) * grid.dzi

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
# using Revise
using FastGroupBy, StatsBase, DataFrames, Random
# using SortingAlgorithms
import DataFrames: DataFrame
using Test
using CategoricalArrays, PooledArrays

# include("fgroupreduce_fby.jl")

N=1_000_000; K=100;

# fastby for CategoricalArrays and PooledArrays

# TODO
pools = unique([randstring(rand(1:32)) for i = 1:N÷K]);
byvec = CategoricalArray{String, 1}(rand(UInt32(1):UInt32(length(pools)), N), CategoricalPool(pools, false));
valvec = rand(N);
@time cmres = countmap(byvec, weights(valvec)); #76
@time fnrs = fastby(sum, byvec, valvec); # 24

fnrsdict = Dict(zip(fnrs[1], fnrs[2]))

@test length(cmres) == length(fnrs[1])
@test all([cmres[k] ≈ fnrsdict[k] for k in keys(cmres)])

byvec = PooledArray(String.(byvec));
@time cmres = countmap(byvec, weights(valvec)); #38
@time fnrs = fastby(sum, byvec, valvec); # 21

fnrsdict = Dict(zip(fnrs[1], fnrs[2]))
@test length(cmres) == length(fnrsdict)
@test all([cmres[k] ≈ fnrsdict[k] for k in keys(cmres)])

# test sum
b = rand(1:1_000_000, 100_000_000)
@time res = fastby(sum, b, [1 for i in 1:100_000_000])
@time rescm = countmap(b)

resdict = Dict(zip(res...))

@test length(rescm) == length(resdict)
@test all([rescm[k] == resdict[k] for k in keys(resdict)])

# test mean
@time res = fastby(mean, b, [1 for i in 1:100_000_000])

resdict = Dict(zip(res...))

@test 1_000_000 == length(resdict)
@test all([1.0 == resdict[k] for k in keys(resdict)])

# Basic sumby and fastby!
a = [1, 1, 2, 3, 3];
aaa = fastby!(sum, a, copy(a))
b = Dict(1 => 2, 2 => 2, 3 => 6)
aaab = Dict(a=>b for (a,b) in zip(aaa...))
@test aaab == b

# Basic sumby and fastby
a = [1, 1, 2, 3, 3];
aaa = fastby(sum, a, copy(a))
b = Dict(1 => 2, 2 => 2, 3 => 6)
aaab = Dict(a=>b for (a,b) in zip(aaa...))
@test aaab == b

# fastby! sum
byvec  = [88, 888, 8, 88, 888, 88]
valvec = [1 , 2  , 3, 4 , 5  , 6]
grpsum = fastby!(sum, byvec, valvec)
expected_result = Dict(88 => 11, 8 => 3, 888 => 7)
@test Dict(a=>b for (a,b) in zip(grpsum...)) == expected_result

# fastby sum
byvec  = [88, 888, 8, 88, 888, 88]
valvec = [1 , 2  , 3, 4 , 5  , 6]
grpsum = fastby(sum, byvec, valvec)
expected_result = Dict(88 => 11, 8 => 3, 888 => 7)
@test Dict(a=>b for (a,b) in zip(grpsum...)) == expected_result

# fastby! string
byvec  = ["grpA", "grpC", "grpB", "grpA", "grpC", "grpA"]
valvec = [1     , 2     , 3     , 4     , 5     , 6     ]
grpsum = fastby!(sum, byvec, valvec)
expected_result = Dict("grpA" => 11, "grpB" => 3, "grpC" => 7)
@test Dict(a=>b for (a,b) in zip(grpsum...)) == expected_result

# fastby string
byvec  = ["grpA", "grpC", "grpB", "grpA", "grpC", "grpA"]
valvec = [1     , 2     , 3     , 4     , 5     , 6     ]
grpsum = fastby(sum, byvec, valvec)
expected_result = Dict("grpA" => 11, "grpB" => 3, "grpC" => 7)
@test Dict(a=>b for (a,b) in zip(grpsum...)) == expected_result

# fastby! Bool
byvec = rand(Bool, 1_000_000)
valvec = rand(1_000_000)
x = fastby!(sum, byvec, valvec)
y = countmap(byvec, weights(valvec))
@test length(Dict(zip(x...)))  == length(y) && [Dict(zip(x...))[k] ≈ y[k] for k = keys(y)] |> all

# fastby Bool
byvec = rand(Bool, 1_000_000)
valvec = rand(1_000_000)
x = fastby(sum, byvec, valvec)
y = countmap(byvec, weights(valvec))
@test length(Dict(zip(x...)))  == length(y) && [Dict(zip(x...))[k] ≈ y[k] for k = keys(y)] |> all

# sumby vs DataFrames.aggregate
Random.seed!(1);
id = rand(1:Int(round(N/K)), N);
val = rand(round.(rand(K)*100, digits = 4), N);
df = DataFrame(id = id, val = val);
@time x = DataFrames.aggregate(df, :id, mean); # 3.3 seconds
#@time x = by(df, :id, val = :val => sum); # 3.3 seconds
xdict = Dict(zip(x[!, :id], x[!, :val_mean]))
res = fastby(mean, df, :id, :val)
y = Dict(a=>b for (a,b) in zip(res[!,:id], res[!,:V1]))
length(xdict) == length(y) && [xdict[k] ≈ y[k] for k in keys(xdict)] |> all
#end


text:
"""
    struct immutable CartesianSpace <: EuclideanSpace{ℝ}

A CartesianSpace is a real euclidean space `ℝ^d` and therefore self-dual. It has no
additonal structure and is completely characterised by its dimension `d`. This is
the vector space that is implicitly assumed in most of matrix algebra.
"""
struct CartesianSpace <: EuclideanSpace{ℝ}
    d::Int
end

Base.getindex(::RealNumbers) = CartesianSpace
Base.getindex(::RealNumbers, d::Int) = CartesianSpace(d)
Base.:^(::RealNumbers, d::Int) = CartesianSpace(d)

# Corresponding methods:
#------------------------
dim(V::CartesianSpace) = V.d
Base.indices(V::CartesianSpace) = Base.OneTo(dim(V))

Base.show(io::IO, V::CartesianSpace) = print(io, "ℝ^$(V.d)")

# direct sum
⊕(V1::CartesianSpace, V2::CartesianSpace) = CartesianSpace(V1.d+V2.d)


text:
import BufferedStreams, Zlib, GZip, Libz
include("gzbufferedstream.jl")

# Source:
# ftp://ftp.ensembl.org/pub/release-81/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz
const filename = "Homo_sapiens.GRCh38.cds.all.fa"
const gzfilename = "Homo_sapiens.GRCh38.cds.all.fa.gz"

const methods = [
    ("IOStream",                      () -> open(filename)),
    ("BufferedInputStream/IOStream",  () -> BufferedStreams.BufferedInputStream(open(filename))),
    ("BufferedInputStream/Mmap",      () -> BufferedStreams.BufferedInputStream(
                                                Mmap.mmap(open(filename), Vector{UInt8}, (filesize(filename),)))),
    ("IOBuffer/Mmap",                 () -> IOBuffer(Mmap.mmap(open(filename), Vector{UInt8}, (filesize(filename),)))),
    ("GZip",                          () -> GZip.gzopen(gzfilename)),
    ("Zlib",                          () -> Zlib.Reader(open(gzfilename))),
    ("Libz",                          () -> Libz.ZlibInflateInputStream(open(gzfilename))),
    ("Libz/Mmap",                     () -> Libz.ZlibInflateInputStream(
                                                 Mmap.mmap(open(gzfilename), Vector{UInt8}, (filesize(gzfilename),)))),
    ("GZBufferedStream",              () -> GZBufferedStream(gzopen(gzfilename))),
    ("Pipe/gzip",                     () -> open(`gzip -cd $gzfilename`)[1]),
    #("BufferedInputStream/Pipe/gzip", () -> BufferedStreams.BufferedInputStream(open(`gzip -cd $gzfilename`)[1])),
]

function bench(f)
    for (method_name, open_func) in methods
        print("  ", method_name, ":  ")
        f(open_func())
        println(round(f(open_func()), 4))
    end
end

function bench_readbytes(input)
    tic()
    readbytes(input)
    return toq()
end

function bench_read(input)
    tic()
    while !eof(input)
        read(input, UInt8)
    end
    return toq()
end

function bench_readline(input)
    tic()
    for line in eachline(input)
    end
    return toq()
end

println("readbytes")
bench(bench_readbytes)

println("read")
bench(bench_read)

println("readline")
bench(bench_readline)


text:
# This file was generated, do not modify it.

using MLJ, PrettyPrinting
X, y = @load_iris
@load DecisionTreeClassifier

dtc = DecisionTreeClassifier()
r   = range(dtc, :max_depth, lower=1, upper=5)

tm = TunedModel(model=dtc, ranges=[r, ], measure=cross_entropy)

m = machine(tm, X, y)
fit!(m)

fitted_params(m).best_model.max_depth

tm = TunedModel(model=dtc, ranges=r, operation=predict_mode,
                measure=misclassification_rate)
m = machine(tm, X, y)
fit!(m)
fitted_params(m).best_model.max_depth

r = report(m)
r.best_measurement

using PyPlot
figure(figsize=(8,6))
plot(r.parameter_values, r.measurements)

xticks(1:5, fontsize=12)
yticks(fontsize=12)
xlabel("Maximum depth", fontsize=14)
ylabel("Misclassification rate", fontsize=14)
ylim([0, 1])

savefig("assets/literate/A-model-tuning-hpt.svg") # hide

X = (x1=rand(100), x2=rand(100), x3=rand(100))
y = 2X.x1 - X.x2 + 0.05 * randn(100);

dtr = @load DecisionTreeRegressor
forest = EnsembleModel(atom=dtr)

params(forest) |> pprint

r1 = range(forest, :(atom.n_subfeatures), lower=1, upper=3)
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0)
tm = TunedModel(model=forest, tuning=Grid(resolution=12),
                resampling=CV(nfolds=6), ranges=[r1, r2],
                measure=rms)
m = machine(tm, X, y)
fit!(m);

r = report(m)
r.best_measurement

figure(figsize=(8,6))

vals_sf = r.parameter_values[:, 1]
vals_bf = r.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, r.measurements)
xlabel("Number of sub-features", fontsize=14)
ylabel("Bagging fraction", fontsize=14)
xticks([1, 2, 3], fontsize=12)
yticks(fontsize=12)

savefig("assets/literate/A-model-tuning-hm.svg") # hide


text:
module JSONLines

import JSON3,
    Mmap,
    Tables,
    StructTypes
import CategoricalArrays

import Base.Threads.@spawn

import Base:
    length,
    size,
    iterate,
    IteratorSize,
    eltype,
    getindex,
    firstindex,
    lastindex,
    IndexStyle,
    checkbounds,
    summary,
    filter,
    findall,
    findnext,
    findfirst,
    findprev,
    findlast

import Tables:
    istable,
    rowaccess,
    rows

import OrderedCollections:
    OrderedDict

export writelines,
    @MStructType,
    readcols,
    columnnames,
    LineIndex,
    LineIterator,
    materialize,
    columnwise,
    gettypes,
    settypes!,
    settype!,
    columntypes,
    filter,
    findall,
    findnext,
    findfirst,
    findprev,
    findlast

include("LineIterator.jl")
include("LineIndex.jl")
include("utils.jl")
include("sugar.jl")
include("write.jl")
end


text:
val = UncertainValue([1, 2, 3], [0.3, 0.3, 0.3])
val2 = UncertainValue([1, 2, 3], [0.3, 0.3, 0.3])
vals = UncertainValue([val2, val], [1, 2])
val3 = UncertainValue(Normal, 0, 1)

###################################
# Single values
###################################
x = zeros(Float64, 10)

# There should be no zero entries after resampling in place
resample!(x, val)
@test any(x .== 0) == false

resample!(x, vals)
@test any(x .== 0) == false

resample!(x, val3)
@test any(x .== 0) == false

###################################
# Collections into vectors
###################################
v = fill(NaN, length(example_uvals))
resample!(v, example_uvals)
@test any(isnan.(v)) == false

###################################
# Collections into arrays
###################################
arr = fill(NaN, 100, length(example_uvals))
resample!(arr, example_uvals)
@test any(isnan.(arr)) == false

#################################################
# Uncertain index-value collections into vectors
#################################################
ids = UncertainIndexDataset(example_uidxs)
vals = UncertainValueDataset(example_uvals)
U = UncertainIndexValueDataset(ids, vals)

idxs = fill(NaN, length(U))
vals = fill(NaN, length(U))

resample!(idxs, vals, U)

@test any(isnan.(idxs)) == false
@test any(isnan.(vals)) == false

#################################################
# Uncertain index-value collections into arrays
#################################################
ids = UncertainIndexDataset(example_uidxs)
vals = UncertainValueDataset(example_uvals)
U = UncertainIndexValueDataset(ids, vals)

n_draws = 10
idxs = fill(NaN, n_draws, length(U))
vals = fill(NaN, n_draws, length(U))

resample!(idxs, vals, U)

@test any(isnan.(idxs)) == false
@test any(isnan.(vals)) == false

###################################################################################
# Draw N realisations of an uncertain value into vector-like containers of length N
###################################################################################

# A single uncertain value resampled multiple times into a N-element vector
N = 10
x = fill(NaN, 10)

resample!(x, UncertainValue(Normal(0, 1)))
@test any(isnan.(x)) == false

# A single uncertain value resampled multiple times into a N-element MVector
N = 10
x = MVector{N, Float64}(repeat([NaN], N))

resample!(x, UncertainValue(Normal(0, 1)))
@test any(isnan.(x)) == false

# A single uncertain value resampled multiple times into a 3-element FieldVector-type
mutable struct Vector3DType <: FieldVector{3, Float64}
    x::Float64
    y::Float64
    z::Float64
end

x = Vector3DType(NaN, NaN, NaN)

resample!(x, UncertainValue(Normal(0, 1)))
@test any(isnan.(x)) == false

#########################################################################################
# Draw single realisations of N uncertain values into vector-like containers of length N
#########################################################################################

# Three uncertain values resampled element-wise into a 3-element vector
x = repeat([NaN], 3)
resample!(x, (val, val, val))
@test any(isnan.(x)) == false

# Three uncertain values resampled element-wise into a 3-element MVector
N = 3
x = MVector{N, Float64}(repeat([NaN], N))
uvals = [UncertainValue([0], [1]), UncertainValue([1], [1]), UncertainValue([2], [1])]
resample!(x, uvals)
@test any(isnan.(x)) == false

# when the number of elements does not match the number of uncertain values
x = MVector{2, Float64}(repeat([NaN], 2))
uval = UncertainValue(Normal(0, 1))
resample!(x, uval)

# Two uncertain values resampled elementwise into a 2-element vector-like type
mutable struct VectorLikeType <: FieldVector{2, Float64}
    x::Float64
    y::Float64
end

x = VectorLikeType(NaN, NaN)
uvals = [UncertainValue([0], [1]), UncertainValue([1], [1])]
resample!(x, uvals)
@test any(isnan.(x)) == false


text:
using CSV
using DataFrames
using Gurobi
using Kombinator
using LightGraphs

n_repetitions = 10

run_um = true
run_st = true
run_ep = true

# Generic helpers.
function parse_algo_string(str::String)
    if str == "Exact"
        method = "Exact"
        algo = "Irrelevant"
        nls = ExactNonlinearSolver(Gurobi.Optimizer)
    elseif str == "Approx: DP"
        method = "NLCOP"
        algo = "Dynamic programming"
        nls = ApproximateNonlinearSolver(DynamicProgramming())
    elseif str == "Approx: LP"
        method = "NLCOP"
        algo = "Linear programming"
        nls = ApproximateNonlinearSolver(DefaultLinearFormulation(Gurobi.Optimizer))
    else
        error("Unrecognised algorithm: $(str).")
    end
    return method, algo, nls
end

# Uniform matroid.
if run_um
    um = DataFrame(method=String[], algo=String[], name=String[], nonlinearity=String[], size=Int[], repetition=Int[], objective=Float64[], time_ms=Float64[])
    um_sizes = [2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000]
    um_algos = ["Exact", "Approx: DP", "Approx: LP"]

    for repetition in 1:n_repetitions
        for size in um_sizes
            m = max(1, round(Int, size / 3))
            hi = 1.0
            lo = 0.5

            lw = [ifelse(i <= size / 2, hi, lo) for i in 1:size]
            nlw = [ifelse(i <= size / 2, lo, hi) for i in 1:size]

            for um_algo in um_algos
                method, algo, nls = parse_algo_string(um_algo)

                li = UniformMatroidInstance(lw, m, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    0.5,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )

                t0 = time_ns()
                s = solve(nli, nls)
                t1 = time_ns()

                time_ms = (t1 - t0) / 1_000_000
                obj_lin = sum(lw[e] for e in s.variables)
                obj_nl = sum(nlw[e] for e in s.variables)
                obj = obj_lin + sqrt(obj_nl)

                push!(um, Dict(:method => method, :algo => algo, :name => um_algo, :nonlinearity => "Square root", :size => size, :repetition => repetition, :objective => obj, :time_ms => time_ms))
            end
        end
    end

    CSV.write(@__DIR__() * "/um.csv", um)
end

# Spanning tree.
if run_st
    st = DataFrame(method=String[], algo=String[], name=String[], nonlinearity=String[], size=Int[], repetition=Int[], objective=Float64[], time_ms=Float64[])
    st_sizes = [2, 5, 10, 20, 50, 100, 200, 500]
    st_algos = ["Exact", "Approx: DP", "Approx: LP"]

    for repetition in 1:n_repetitions
        for size in st_sizes
            g = complete_graph(size)
            hi = 1.0
            lo = 0.5

            lw = Dict{Edge{Int}, Float64}(e => ifelse(i <= size / 2, hi, lo) for (i, e) in enumerate(edges(g)))
            nlw = Dict{Edge{Int}, Float64}(e => ifelse(i <= size / 2, lo, hi) for (i, e) in enumerate(edges(g)))

            for st_algo in st_algos
                method, algo, nls = parse_algo_string(st_algo)

                li = SpanningTreeInstance(g, lw, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    0.5,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )

                t0 = time_ns()
                s = solve(nli, nls)
                t1 = time_ns()

                time_ms = (t1 - t0) / 1_000_000
                obj_lin = sum(e in keys(lw) ? lw[e] : lw[reverse(e)] for e in s.variables)
                obj_nl = sum(e in keys(nlw) ? nlw[e] : nlw[reverse(e)] for e in s.variables)
                obj = obj_lin + sqrt(obj_nl)

                push!(st, Dict(:method => method, :algo => algo, :name => st_algo, :nonlinearity => "Square root", :size => size, :repetition => repetition, :objective => obj, :time_ms => time_ms))
            end
        end
    end

    CSV.write(@__DIR__() * "/st.csv", st)
end

# Elementary path.
if run_ep
    ep = DataFrame(method=String[], algo=String[], name=String[], nonlinearity=String[], size=Int[], repetition=Int[], objective=Float64[], time_ms=Float64[])
    ep_sizes = [2, 5, 10, 20, 50, 100, 200, 500]
    ep_algos = ["Exact", "Approx: DP", "Approx: LP"]

    for repetition in 1:n_repetitions
        for size in ep_sizes
            g0 = complete_graph(size) # Only for determining the edges
            g = complete_digraph(size)
            hi = 1.0
            lo = 0.5

            lw = Dict{Edge{Int}, Float64}(e => ifelse(i <= size / 2, hi, lo) for (i, e) in enumerate(edges(g0)))
            nlw = Dict{Edge{Int}, Float64}(e => ifelse(i <= size / 2, lo, hi) for (i, e) in enumerate(edges(g0)))

            for ep_algo in ep_algos
                method, algo, nls = parse_algo_string(ep_algo)

                li = ElementaryPathInstance(g, lw, 1, size, Maximise())
                nli = NonlinearCombinatorialInstance(
                    li,
                    lw,
                    nlw,
                    SquareRoot,
                    0.5,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                    true,
                    DefaultLinearFormulation(Gurobi.Optimizer),
                )

                t0 = time_ns()
                s = solve(nli, nls)
                t1 = time_ns()

                time_ms = (t1 - t0) / 1_000_000
                obj_lin = sum(e in keys(lw) ? lw[e] : lw[reverse(e)] for e in s.variables)
                obj_nl = sum(e in keys(nlw) ? nlw[e] : nlw[reverse(e)] for e in s.variables)
                obj = obj_lin + sqrt(obj_nl)

                push!(ep, Dict(:method => method, :algo => algo, :name => ep_algo, :nonlinearity => "Square root", :size => size, :repetition => repetition, :objective => obj, :time_ms => time_ms))
            end
        end
    end

    CSV.write(@__DIR__() * "/ep.csv", ep)
end


text:
# LSP indexing solution
# https://github.com/julia-vscode/julia-vscode/issues/800#issuecomment-650085983
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/SpinDoctor.jl")
    using .SpinDoctor
end

using SpinDoctor
using LinearAlgebra
using GLMakie

## Chose a plotting theme
set_theme!(theme_light())
set_theme!(theme_dark())
set_theme!(theme_black())

## Create model from setup recipe
# include("setups/axon.jl")
# include("setups/sphere.jl")
include("setups/plates.jl")
# include("setups/cylinders.jl")
# include("setups/spheres.jl")
# include("setups/neuron.jl")

mesh, surfaces, cells = @time create_geometry(setup; recreate = true);
model = Model(; mesh, coeffs...);
volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)
@info "Number of nodes per compartment:" length.(model.mesh.points)

## Plot mesh
plot_surfaces(surfaces, 1:3)
plot_mesh(model.mesh, 1:1)

## Assemble finite element matrices
matrices = @time assemble_matrices(model);

## Magnetic field gradient
dir = [1.0, 0.0, 0.0]
profile = PGSE(2000.0, 6000.0)
# profile = CosOGSE(5000.0, 5000.0, 2)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)

## Solve BTPDE

# Callbacks for time stepping (plot solution, save time series)
printer = Printer(; nupdate = 1, verbosity = 2)
writer = VTKWriter(; nupdate = 5)
plotter = Plotter{T}(; nupdate = 5)
# callbacks = [printer, plotter]
callbacks = [printer, plotter, writer]

# Choose BTDPE solver (specialized solver only for PGSE)
solver = IntervalConstantSolver{T}(; θ = 0.5, timestep = 5.0)
solver = QNDF(autodiff = false)

# Solve BTPDE
btpde = BTPDE(; model, matrices)
ξ = @time solve(btpde, gradient, solver; callbacks)

## Plot magnetization
plot_field(model.mesh, ξ)

## Compute signal
compute_signal(matrices.M, ξ)
compute_signal.(matrices.M_cmpts, split_field(model.mesh, ξ))

## Save magnetization
savefield(model.mesh, ξ, "output/magnetization")

## Matrix Formalism

# Perform Laplace eigendecomposition
laplace = Laplace{T}(; model, matrices, neig_max = 400)
lap_eig = @time solve(laplace)
length_scales = eig2length.(lap_eig.values, D_avg)

# Truncate basis at minimum length scale
length_scale = 3
λ_max = length2eig(length_scale, D_avg)
lap_eig = limit_lengthscale(lap_eig, λ_max)

# Compute magnetization using the matrix formalism reduced order model
mf = MatrixFormalism(; model, matrices, lap_eig)
ξ = @time solve(mf, gradient; ninterval = 500)

## Solve HADC
hadc = HADC(; model, matrices)
adc_cmpts = @time solve(hadc, gradient)

## Solve Karger model

# Compute HADC and fit difftensors
directions = unitsphere(50)
gradients = [
    ScalarGradient(collect(d), gradient.profile, gradient.amplitude) for
    d ∈ eachcol(directions)
]
hadc = HADC(; model, matrices)
adcs, = @time solve_multigrad(hadc, gradients)
difftensors = fit_tensors(directions, adcs)

# Solve Karger
karger = Karger(; model, difftensors)
signal = @time solve(karger, gradient; timestep = 5.0)

## Solve analytical model
# Compute analytical Laplace eigenfunctions
length_scale = 0.3
eigstep = 1e-8
eiglim = length2eig(length_scale, D_avg)
analytical_coeffs = analytical_coefficients(setup, coeffs)
analytical_laplace = AnalyticalLaplace(; analytical_coeffs..., eiglim, eigstep)
lap_mat = @time solve(analytical_laplace)

# Compute analytical matrix formalism signal truncation
analytical_mf = AnalyticalMatrixFormalism(; analytical_laplace, lap_mat, volumes)
signal = solve(analytical_mf, gradient)