text:
export PDATA, get_edge_dof,
    impose_Dirichlet_boundary_conditions, dof_to_gauss_points, get_boundary_edge_orientation,
    compute_pml_term, solve_slip_law

"""
    PDATA

The folder where built-in meshes are stored.
"""
PDATA = abspath(joinpath(@__DIR__, "..", "..",  "deps", "MFEM", "MeshData"))

function get_meshio()
    try
        pyimport("meshio")
    catch
        PIP = joinpath(ADCME.BINDIR, "pip")
        run_with_env(`$PIP install meshio==4.2`)
        pyimport("meshio")
    end
end

"""
    Mesh(filename::String; file_format::Union{String, Missing} = missing)

Reads a mesh file and extracts element, coordinates and boundaries.

# Example
```julia
mesh = Mesh(joinpath(PDATA, "twoholes.stl"))
```
"""
function Mesh(filename::String; file_format::Union{String, Missing} = missing,
                order::Int64 = 2, degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = 2)
    if splitext(filename)[2] == ".mat"
        d = matread(filename)
        return Mesh(Float64.(d["nodes"]), Int64.(d["elems"]), order, degree, lorder)
    end
    meshio = AdFem.get_meshio()
    if !ismissing(file_format)
        mesh = meshio.read(filename, file_format = file_format)
    else
        mesh = meshio.read(filename)
    end
    elem = py"list($mesh.cells[0])"[2]
    if length(elem)==0
        error("No triangles found in the mesh file.")
    end
    Mesh(Float64.(mesh.points[:,1:2]), Int64.(elem) .+ 1, order, degree, lorder)
end

"""
    get_edge_normal(mmesh::Mesh)
"""
function get_edge_normal(mmesh::Mesh)
    _edge_to_elem_map = Set{Tuple{Int64, Int64}}()
    dic = (x,y)->(minimum([x,y]), maximum([x,y]))
    add_dic = (e1, e2)->begin
        if dic(e1, e2) in _edge_to_elem_map
            delete!(_edge_to_elem_map, dic(e1, e2))
        else
            push!(_edge_to_elem_map, dic(e1, e2))
        end
    end
    edict = Dict()
    for i = 1:mmesh.nelem
        e1, e2, e3 = mmesh.elems[i,:]
        add_dic(e1, e2); edict[dic(e1,e2)] = e3
        add_dic(e3, e2); edict[dic(e3,e2)] = e1
        add_dic(e1, e3); edict[dic(e1,e3)] = e2
    end
    out = zeros(Float64, length(_edge_to_elem_map), 2)
    for (k,s) in enumerate(_edge_to_elem_map)
        p1 = s[1]
        p2 = s[2]
        p3 = edict[s]
        x1 = mmesh.nodes[p1,:] - mmesh.nodes[p3,:]
        x2 = mmesh.nodes[p1,:] - mmesh.nodes[p2,:]
        n = [x2[2]; -x2[1]]
        n = n/norm(n)
        if dot(n, x1)<0
            n = -n
        end
        out[k,:] = n
    end
    out
end

"""
    get_edge_normal(bdedge::Array{Int64, 2}, mmesh::Mesh)

Returns the outer normal for all edges in `bdedge`.
"""
function get_edge_normal(bdedge::Array{Int64, 2}, mmesh::Mesh)
    _edge_to_elem_map = Set{Tuple{Int64, Int64}}()
    dic = (x,y)->(minimum([x,y]), maximum([x,y]))
    add_dic = (e1, e2)->begin
        if dic(e1, e2) in _edge_to_elem_map
            delete!(_edge_to_elem_map, dic(e1, e2))
        else
            push!(_edge_to_elem_map, dic(e1, e2))
        end
    end
    edict = Dict()
    for i = 1:mmesh.nelem
        e1, e2, e3 = mmesh.elems[i,:]
        add_dic(e1, e2); edict[dic(e1,e2)] = e3
        add_dic(e3, e2); edict[dic(e3,e2)] = e1
        add_dic(e1, e3); edict[dic(e1,e3)] = e2
    end
    out = zeros(Float64, size(bdedge, 1), 2)
    for k = 1:size(bdedge, 1)
        s = bdedge[k,:]
        p1 = s[1]
        p2 = s[2]
        p3 = edict[dic(p1, p2)]
        x1 = mmesh.nodes[p1,:] - mmesh.nodes[p3,:]
        x2 = mmesh.nodes[p1,:] - mmesh.nodes[p2,:]
        n = [x2[2]; -x2[1]]
        n = n/norm(n)
        if dot(n, x1)<0
            n = -n
        end
        out[k,:] = n
    end
    out
end

"""
    save(filename::String, mesh::Mesh)

Saves the mesh to the file `filename`.
"""
function save(filename::String, mesh::Mesh)
    matwrite(filename, Dict(
        "nodes"=>mesh.nodes, "elems"=>mesh.elems
    ))
end

"""
    get_edge_dof(edges::Array{Int64, 2}, mesh::Mesh)
    get_edge_dof(edges::Array{Int64, 1}, mesh::Mesh)

Returns the DOFs for `edges`, which is a `K × 2` array containing vertex indices.
The DOFs are not offset by `nnode`, i.e., the smallest edge DOF could be 1.

When the input is a length 2 vector, it returns a single index for the corresponding edge DOF.
"""
function get_edge_dof(edges::Array{Int64, 2}, mesh::Mesh)
    d = Dict{Tuple{Int64, Int64}, Int64}()
    for i = 1:mesh.nedge
        m = minimum(mesh.edges[i,:])
        M = maximum(mesh.edges[i,:])
        d[(m, M)] = i
    end
    idx = Int64[]
    for i = 1:size(edges, 1)
        m = minimum(edges[i,:])
        M = maximum(edges[i,:])
        push!(idx, d[(m, M)])
    end
    idx
end

function get_edge_dof(edges::Array{Int64, 1}, mesh::Mesh)
    idx = get_edge_dof(reshape(edges, 1, 2), mesh)
    return idx[1]
end

@doc raw"""
    impose_Dirichlet_boundary_conditions(A::Union{SparseArrays, Array{Float64, 2}}, rhs::Array{Float64,1}, bdnode::Array{Int64, 1},
        bdval::Array{Float64,1})
    impose_Dirichlet_boundary_conditions(A::SparseTensor, rhs::Union{Array{Float64,1}, PyObject}, bdnode::Array{Int64, 1},
        bdval::Union{Array{Float64,1}, PyObject})

Algebraically impose the Dirichlet boundary conditions. We want the solutions at indices `bdnode` to be `bdval`. Given the matrix and the right hand side

$$\begin{bmatrix} A_{II} & A_{IB} \\ A_{BI} & A_{BB} \end{bmatrix}, \begin{bmatrix}r_I \\ r_B \end{bmatrix}$$

The function returns

$$\begin{bmatrix} A_{II} & 0 \\ 0 & I \end{bmatrix}, \begin{bmatrix}r_I - A_{IB} u_B \\ r_B \end{bmatrix}$$
"""
function impose_Dirichlet_boundary_conditions(A::Union{SparseMatrixCSC, Array{Float64, 2}},
    rhs::Array{Float64,1}, bdnode::Array{Int64, 1},
    bdval::Array{Float64,1})
    N = length(rhs)
    r = copy(rhs)
    idx = ones(Bool, N)
    idx[bdnode] .= false
    A11 = A[idx, idx]
    A12 = A[idx, bdnode]
    r[idx] = r[idx] - A12 * bdval
    r[bdnode] = bdval
    B = spzeros(N, N)
    B[idx, idx] = A11
    B[bdnode, bdnode] = spdiagm(0=>ones(length(bdnode)))
    B, r
end

function impose_Dirichlet_boundary_conditions(A::SparseTensor, rhs::Union{Array{Float64,1}, PyObject},
    bdnode::Array{Int64, 1},
    bdval::Union{Array{Float64,1}, PyObject})
    indices = A.o.indices
    vv = A.o.values
    @assert size(A, 1) == size(A, 2) == length(rhs)
    @assert length(bdnode)==length(bdval)
    @assert length(bdnode)<=length(rhs)
    impose_dirichlet_ = load_op_and_grad(AdFem.libmfem,"impose_dirichlet", multiple=true)
    indices,vv,bd,rhs,bdval = convert_to_tensor(Any[indices,vv,bdnode,rhs,bdval], [Int64,Float64,Int64,Float64,Float64])
    indices, vv, rhs = impose_dirichlet_(indices,vv,bd,rhs,bdval)
    RawSparseTensor(indices, vv, size(A)...), set_shape(rhs, (size(A,2),))
end

"""
    impose_Dirichlet_boundary_conditions(A::SparseTensor, bdnode::Array{Int64, 1})

A helper function to impose homogeneous Dirichlet boundary condition.
"""
function impose_Dirichlet_boundary_conditions(A::SparseTensor, bdnode::Array{Int64, 1})
    rhs = zeros(size(A, 1))
    bdval = zeros(length(bdnode))
    B, _ = impose_Dirichlet_boundary_conditions(A, rhs, bdnode, bdval)
    B
end

"""
    fem_to_gauss_points(u::PyObject, mesh::Mesh)
"""
function fem_to_gauss_points(u::PyObject, mesh::Mesh)
    fem_to_gauss_points_mfem_ = load_op_and_grad(AdFem.libmfem,"fem_to_gauss_points_mfem")
    u = convert_to_tensor(Any[u], [Float64]); u = u[1]
    out = fem_to_gauss_points_mfem_(u)
    set_shape(out, (get_ngauss(mesh),))
end

"""
    fem_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)
"""
function fem_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)
    out = zeros(get_ngauss(mesh))
    @eval ccall((:FemToGaussPointsMfem_Julia, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}),
                $out, $u)
    return out
end

"""
    dof_to_gauss_points(u::PyObject, mesh::Mesh)
    dof_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)

Similar to [`fem_to_gauss_points`](@ref). The only difference is that the function uses all DOFs---which means,
for quadratic elements, the nodal values on the edges are also used.
"""
function dof_to_gauss_points(u::PyObject, mesh::Mesh)
    @assert length(u)==mesh.ndof
    dof_to_gauss_points_mfem_ = load_op_and_grad(AdFem.libmfem,"dof_to_gauss_points_mfem")
    u = convert_to_tensor(Any[u], [Float64]); u = u[1]
    out = dof_to_gauss_points_mfem_(u)
    set_shape(out, (get_ngauss(mesh),))
end

function dof_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)
    @assert length(u)==mesh.ndof
    out = zeros(get_ngauss(mesh))
    @eval ccall((:DofToGaussPointsMfem_forward_Julia, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}),
                $out, $u)
    return out
end

"""
    get_boundary_edge_orientation(bdedge::Array{Int64, 2}, mmesh::Mesh)

Returns the orientation of the edges in `bdedge`. For example, if for a boundary element `[1,2,3]`, assume `[1,2]` is the boundary edge,
then

```
get_boundary_edge_orientation([1 2;2 1], mmesh) = [1.0;-1.0]
```

The return values for non-boundary edges in `bdedge` is undefined.
"""
function get_boundary_edge_orientation(bdedge::Array{Int64, 2}, mmesh::Mesh)
    edge_orientation = Dict{Tuple{Int64, Int64}, Float64}()
    elems = mmesh.elems
    bd = bcedge(mmesh)
    bd = Set([(bd[i,1], bd[i,2]) for i = 1:size(bd, 1)])
    add_to_dict = (k, i, j)->begin
        if (elems[k, i], elems[k, j]) in bd || (elems[k, j], elems[k, i]) in bd
            edge_orientation[(elems[k, i], elems[k, j])] = 1.0
            edge_orientation[(elems[k, j], elems[k, i])] = -1.0
        end
    end
    for i = 1:mmesh.nelem
        add_to_dict(i, 1, 2)
        add_to_dict(i, 2, 3)
        add_to_dict(i, 3, 1)
    end
    sng = zeros(size(bdedge, 1))
    for i = 1:size(bdedge, 1)
        if !haskey(edge_orientation, (bdedge[i,1], bdedge[i,2]))
            error("($(bdedge[i,1]), $(bdedge[i,2])) is not a boundary edge")
        end
        sng[i] = edge_orientation[(bdedge[i,1], bdedge[i,2])]
    end
    sng
end

@doc raw"""
    compute_pml_term(u::Union{Array{Float64,1}, PyObject},βprime::Union{Array{Float64,1}, PyObject},
    c::Union{Array{Float64,1}, PyObject},nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh)

- `u`: a vector of length `mmesh.ndof`
- `βprime`: a vector of length $n_{\text{gauss}}$
- `c`: a tensor of size $n_{\text{gauss}}$
- `nv`: a matrix of size $n_{\text{gauss}}\times 2$

This function returns four outputs

$$\begin{aligned}k_1&=(c^2 n\partial_n u, n\partial_n\delta u)\\k_2&=(\beta'n\cdot(c^2 n\partial_n u), \delta u)\\ k_3&=(c^2\nabla^\parallel u, n\partial_n \delta u) + (c^2 n\partial_ n u, \nabla^\parallel \delta u)\\ k_4&= (c^2 \nabla^\parallel u, \nabla^\parallel \delta u)\end{aligned}$$
"""
function compute_pml_term(u::Union{Array{Float64,1}, PyObject},βprime::Union{Array{Float64,1}, PyObject},
        c::Union{Array{Float64,1}, PyObject},nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh)
    @assert length(u)==mmesh.ndof
    @assert length(βprime)==length(c)==size(nv,1)==get_ngauss(mmesh)
    compute_pml_term_ = load_op_and_grad(AdFem.libmfem,"compute_pml_term", multiple=true)
    u,betap,c,nv = convert_to_tensor(Any[u,βprime,c,nv], [Float64,Float64,Float64,Float64])
    out = compute_pml_term_(u,betap,c,nv)
    set_shape(out[1], (mmesh.ndof, )), set_shape(out[2], (mmesh.ndof, )), set_shape(out[3], (mmesh.ndof, )), set_shape(out[4], (mmesh.ndof, ))
end

@doc raw"""
    compute_pml_term(u::Union{Array{Float64,1}, PyObject},βprime::Union{Array{Float64,1}, PyObject},
    λ::Union{Array{Float64,1}, PyObject},μ::Union{Array{Float64,1}, PyObject}, nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh)

- `u`: a vector of length `mmesh.ndof`
- `βprime`: a vector of length $n_{\text{gauss}}$
- `λ`, `μ`: Lam\'e parameters
- `nv`: a matrix of size $n_{\text{gauss}}\times 2$

This is a 2D version of [`compute_pml_term`](@ref).
"""
function compute_pml_term(u::Union{Array{Float64,1}, PyObject},βprime::Union{Array{Float64,1}, PyObject},
    λ::Union{Array{Float64,1}, PyObject},μ::Union{Array{Float64,1}, PyObject}, nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh)
    @assert length(u)==2mmesh.ndof
    @assert length(βprime)==length(λ)==length(μ)==size(nv,1)==get_ngauss(mmesh)
    compute_pml_elastic_term_ = load_op_and_grad(AdFem.libmfem,"compute_pml_elastic_term", multiple=true)
    u,betap,λ,μ,nv = convert_to_tensor(Any[u,βprime,λ,μ,nv], [Float64,Float64,Float64,Float64,Float64])
    out = compute_pml_elastic_term_(u,betap,λ,μ,nv)
    set_shape(out[1], (2mmesh.ndof, )), set_shape(out[2], (2mmesh.ndof, )), set_shape(out[3], (2mmesh.ndof, )), set_shape(out[4], (2mmesh.ndof, ))
end

@doc raw"""
    solve_slip_law(v, ψ, dc, v0, a, b, f0, Δt::Float64)

Solves one step of the slip law equation

$$\dot \psi = - \frac{|V|}{d_C}\left( a \sinh^{-1}\left( \frac{|V|}{2V_0}e^{\frac{\psi}{a}} \right) - f_0 + (b-a)*\log \frac{|V|}{V_0} \right)$$

We discretize the equation with a central difference scheme

$$\frac{\psi^{n+1} - \psi^{n-1}}{2\Delta t} =  - \frac{|V|}{d_c}\left( a \sinh^{-1} \left( \frac{|V|}{2V_0} e^{\frac{\psi^{n+1} + \psi^{n-1}}{2a}{}} \right) - f_0 + (b-a) \log \frac{|V|}{V_0} \right)$$

- `dc`, `v0`, `a`, `b`, and `f0` are scalars
"""
function solve_slip_law(v, ψ, dc, v0, a, b, f0, Δt::Float64)
    @assert size(v0)==size(a)==size(b)==size(f0)==size(dc)==()
    @assert size(v)==size(ψ)
    @assert length(size(v))==length(size(ψ))==1
    v, ψ, dc, v0, a, b, f0 = convert_to_tensor(Any[v, ψ, dc, v0, a, b, f0], [Float64, Float64, Float64, Float64, Float64, Float64, Float64])
    v = abs(v)
    A = (-Δt*v)/dc
    B = v/2v0*exp(ψ/2a)
    C = (ψ + 2Δt*v/dc*(f0 - (b-a)*log(v/v0)))/2a
    solve_slip_law(A, B, C, ψ)
end

@doc raw"""
    solve_slip_law(
        A::Union{Array{Float64,1}, PyObject},
        B::Union{Array{Float64,1}, PyObject},
        C::Union{Array{Float64,1}, PyObject},
        X0::Union{Array{Float64,1}, PyObject})

A helper function for [`solve_slip_law`](@ref) the nonlinear equation

$$x - A\sinh^{-1}(Bx) - C = 0$$

`A`, `B`, and `C` are vectors of the same length. `X0` is the initial guess.
"""
function solve_slip_law(
        A::Union{Array{Float64,1}, PyObject},
        B::Union{Array{Float64,1}, PyObject},
        C::Union{Array{Float64,1}, PyObject},
        X0::Union{Array{Float64,1}, PyObject})
    @assert length(X0) == length(A) == length(B) == length(C)
    solve_slip_law_ = load_op_and_grad(AdFem.libmfem,"solve_slip_law")
    a,b,c,xinit = convert_to_tensor(Any[A, B, C, X0], [Float64,Float64,Float64,Float64])
    out = solve_slip_law_(a,b,c,xinit)
    set_shape(out, (length(a),))
end


text:
struct LessEqual{L, T <: Integer} <: Integer
    val::T

    LessEqual{L}(x::T) where {L, T} = new{L, T}(x)
end

import Base: +
addone(x::LessEqual{L}) where L = x.val == L ? LessEqual{L}(1) : LessEqual{L}(x + 1)
+(lhs::LessEqual{L}, rhs::Integer) where L = rhs == 1 ? addone(lhs) : LessEqual{L}(lhs.val + rhs)

"""
    fastmod1(::Val{X}, ::Val{Y})

`mod1` in compile time, 40x faster if the value can be determined
"""
@generated fastmod1(::Val{X}, ::Val{Y}) where X where Y =
    X == 0   ? Y :
    X <= Y   ? X :
    X == Y+1 ? 1 :
    X == Y+2 ? 2 :
    mod1(X, Y)

ext: jl


text:
module Slurm

include("Directives.jl")
include("Commands.jl")
include("SystemModules.jl")
include("Shells.jl")
include("Scriptify.jl")

end # module


text:
mutable struct GPSSim
    # Ephemeris data initial and final epochs
    startGPSWeek::Int
    startWeekDay::Int
    endGPSWeek::Int
    endWeekDay::Int

    # Raw gps data with ephemeris in IGS reference frame (currently IGS14)
    rawSp3Data::DataFrame

    # GPS Processed Data
    data::Matrix{Float64}

    # Interpolators
    interps::DataFrame

    # ERPs
    erps::Matrix{Float64}

    # Measurement statistics
    σρ::Float64
    σr::Float64

    function GPSSim(startGPSWeek::Integer, startWeekDay::Integer, endGPSWeek::Integer, endWeekDay::Integer;
                    σρ::Float64 = 1.0e-3, σr::Float64 = 5e-3)
        # Read in sp3 unprocessed data
        rawSp3Data = readSP3s(startGPSWeek, startWeekDay, endGPSWeek, endWeekDay)

        # Read in Earth Orientation Parameters
        erps = readERPs(startGPSWeek, startWeekDay, endGPSWeek, endWeekDay)

        # Rotate data from ITRF14 to GCRF and data required for gps simulation in fast to read matrix
        data = rotateData(rawSp3Data, erps)

        # Create interpolators
        interps = genInterpolants(data)

        new(startGPSWeek, startWeekDay, endGPSWeek, endWeekDay, rawSp3Data, data, interps, erps, σρ, σr)
    end
end

function gpsMeasurement(gps::GPSSim, r::AbstractVector, t; type = :true)
    # Grab vectors of interpolators
    Xs = gps.interps.X
    Ys = gps.interps.Y
    Zs = gps.interps.Z

    # Compute time of transmition from each satelite
    tts = zeros(32)
    for i in 1:32
        try
            func(tt) = transmissionTimeFunction(tt, t, r, Xs[i], Ys[i], Zs[i])
            tts[i] = find_zero(func, t)
            #tts[i] = transmissionTimeFunction(t, r, Xs[i], Ys[i], Zs[i])
        catch
            tts[i] = 0.0
        end
    end

    # Determine visable GPS satelites
    visable = Vector{Bool}(undef, 32); visable .= false
    for i in 1:32
        if tts[i] != 0
            # Get vector to GPS Sat if
            rsat = @SVector [Xs[i](tts[i]), Ys[i](tts[i]), Zs[i](tts[i])]

            # Compute vector from r to rsat
            rdiff = @SVector [rsat[1] - r[1], rsat[2] - r[2], rsat[3] - r[3]]

            # Check for intersection with sphere by solving quadratic equation
            nrDiff = norm(rdiff)
            u = @SVector [rdiff[1] / nrDiff, rdiff[2] / nrDiff, rdiff[3] / nrDiff]
            b = 2.0*transpose(u)*r
            c = transpose(r)*r - 6378.137^2
            d = b^2 - 4*c
            # If d is negative, rdiff does not intersect the sphere of the Earth
            if d < 0.0
                visable[i] = true
            end
        end
    end

    # Generate measurements from visable satelites
    meas    = zeros(2*sum(visable) + 1)
    meas[1] = t
    count   = 0
    for i in 1:32
        if visable[i] == true
            count += 1

            # Set sat num
            meas[2*(count - 1) + 2] = i

            # Get vector to GPS Sat
            if type == :true # Use true position of satelite
                rsat = @SVector [Xs[i](tts[i]), Ys[i](tts[i]), Zs[i](tts[i])]
            else # If generating expected measurements, corrupt position to simulate ArgumentError
                 # in GPS broadcast ephemeris
                rsat = @SVector [Xs[i](tts[i]) + gps.σr*randn(),
                                 Ys[i](tts[i]) + gps.σr*randn(),
                                 Zs[i](tts[i]) + gps.σr*randn()]
            end

            # Compute vector from r to rsat
            rdiff = @SVector [rsat[1] - r[1], rsat[2] - r[2], rsat[3] - r[3]]

            # Get GPS clock bias and corrupt with noise
            Δtsat = 0.0
            Δtm = 0.0
            n, m  = size(gps.data)
            for j in 32 + i:31:n
                if gps.data[j - 31,1] < t && gps.data[j, 1] > t
                    #t1      = gps.data[j - 31,1]
                    #t2      = gps.data[j,1]
                    Δt1     = gps.data[j - 31,6]
                    #Δt2     = gps.data[j,6]
                    Δtσ1    = gps.data[j - 31,10]
                    #Δtσ2    = gps.data[j,10]
                    #Δttrue  = (Δt1 + (t - t1)*(Δt2 - Δt1)/(t2 - t1)) * 1e-6
                    #Δtσ     = Δtσ1 + (t - t1)*(Δtσ2 - Δtσ1)/(t2 - t1) * 1e-12
                    Δtm     = Δt1*1e-6
                    Δtσ     = Δtσ1*1e-12
                    Δtsat   = Δtm + Δtσ*randn()
                    if isnan(Δtsat) # Currently data in GPSSim.data is getting set to NaN so this is a workaround for now.
                        Δtsat = 0.0
                    end
                end
            end

            # Compute pseudorange
            if type == :true
                meas[2*(count - 1) + 3] = norm(rdiff) - 299792.458*Δtsat + gps.σρ*randn()
            else
                meas[2*(count - 1) + 3] = norm(rdiff) - 299792.458*Δtm
            end
        end
    end

    return meas
end

function gpsMeasurement(gps::GPSSim, r::AbstractVector, gpsSats, t; type = :true)
    # Grab vectors of interpolators
    Xs = gps.interps.X
    Ys = gps.interps.Y
    Zs = gps.interps.Z

    # Compute time of transmition from each satelite
    tts = zeros(32)
    for ifloat in gpsSats
        i = Int(ifloat)
        try
            #func(tt) = transmissionTimeFunction(tt, t, r, Xs[i], Ys[i], Zs[i])
            #tts[i] = find_zero(func, t)
            tts[i] = transmissionTimeFunction(t,r,Xs[i],Ys[i],Zs[i])
        catch
            tts[i] = 0.0
        end
    end

    # Generate measurements from visable satelites
    meas    = zeros(2*length(gpsSats) + 1)
    ps      = zeros(length(gpsSats), 3)
    meas[1] = t
    count   = 0
    for ifloat in gpsSats
        i = Int(ifloat)
        count += 1

        # Set sat num
        meas[2*(count - 1) + 2] = i

        # Get vector to GPS Sat
        if type == :true # Use true position of satelite
            rsat = @SVector [Xs[i](tts[i]), Ys[i](tts[i]), Zs[i](tts[i])]
        else # If generating expected measurements, corrupt position to simulate ArgumentError
                # in GPS broadcast ephemeris
            rsat = @SVector [Xs[i](tts[i]) + gps.σr*randn(),
                             Ys[i](tts[i]) + gps.σr*randn(),
                             Zs[i](tts[i]) + gps.σr*randn()]
        end

        # Compute vector from r to rsat
        rdiff = @SVector [rsat[1] - r[1], rsat[2] - r[2], rsat[3] - r[3]]

        # Get GPS clock bias and corrupt with noise
        Δtsat = 0.0
        Δtm   = 0.0
        n, m  = size(gps.data)
        for j in 32 + i:31:n
            if gps.data[j - 31,1] < t && gps.data[j, 1] > t
                #t1      = gps.data[j - 31,1]
                #t2      = gps.data[j,1]
                Δt1     = gps.data[j - 31,6]
                #Δt2     = gps.data[j,6]
                Δtσ1    = gps.data[j - 31,10]
                #Δtσ2    = gps.data[j,10]
                #Δttrue  = (Δt1 + (t - t1)*(Δt2 - Δt1)/(t2 - t1)) * 1e-3
                #Δtσ     = Δtσ1 + (t - t1)*(Δtσ2 - Δtσ1)/(t2 - t1) * 1e-12
                Δtm     = Δt1*1e-6
                Δtσ     = Δtσ1*1e-12
                Δtsat   = Δtm + Δtσ*randn()
            end
        end

        # Compute pseudorange
        if type == :true
            meas[2*(count - 1) + 3] = norm(rdiff) - 299792.458*Δtsat + gps.σρ*randn()
        else
            meas[2*(count - 1) + 3] = norm(rdiff) - 299792.458*Δtm
        end

        # Set GPS Satelite position
        ps[count, :] .= rdiff
    end

    return (meas, ps)
end

function plotGPS(gps::GPSSim, sats::AbstractVector)
    mat"""
    data = $(gps.data);
    sats = $sats;
    satData = nan(floor(size(data,1)/30),3);

    figure()
    hold on
    for i = 1:length(sats)
        satData = satData .* NaN;
        sat     = sats(i);
        count   = 0;
        for j = 1:size(data,1)
            if data(j,2) == sat
                count = count + 1;
                satData(count,1:3)  = data(j,3:5);
            end
        end
        plot3(satData(:,1),satData(:,2),satData(:,3))
    end

    % Plot Earth
    [A,B,C] = ellipsoid(0.0, 0.0, 0.0, 6378, 6378, 6378);
    surf(A,B,-C, imread("./figures/WM.jpg"), "FaceColor", "texturemap", "EdgeColor", "none")

    grid on
    axis equal
    view(3)

    xlabel("X, km", "Interpreter", "latex")
    ylabel("Y, km", "Interpreter", "latex")
    zlabel("Z, km", "Interpreter", "latex")

    set(gca, "fontname", "Times New Roman", "fontsize", 10)
    set(gcf, "PaperUnits", "inches", "PaperPosition", [0.25, 0.25, 5.0, 4.0])
    set(gcf, "PaperPositionMode", "Manual")
    print("./figures/GPSConst.eps", "-depsc", "-r300")
    """
end

plotGPS(gps::GPSSim) = plotGPS(gps, 1:32)

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
using DimensionalData, Test
using DimensionalData: _set, AutoSampling

a = [1 2; 3 4]
dimz = (X((143.0, 145.0); mode=Sampled(order=Ordered()), metadata=Metadata(Dict(:meta => "X"))),
        Y((-38.0, -36.0); mode=Sampled(order=Ordered()), metadata=Metadata(Dict(:meta => "Y"))))
da = DimArray(a, dimz; name=:test)

a2 = [1 2 3 4
      3 4 5 6
      4 5 6 7]
dimz2 = (Dim{:row}(10.0:10.0:30.0), Dim{:column}(-2:1.0:1.0))
da2 = DimArray(a2, dimz2; name=:test2)

s = DimStack(da2, DimArray(2a2, dimz2; name=:test3))

@testset " Array fields" begin
    @test name(set(da2, :newname)) == :newname
    @test name(set(da2, Name(:newname))) == Name{:newname}()
    @test name(set(da2, NoName())) == NoName()
    @test metadata(set(da2, Metadata(Dict(:testa => "test")))).val == Dict(:testa => "test")
    @test parent(set(da2, fill(9, 3, 4))) == fill(9, 3, 4)
    # A differently sized array can't be set
    @test_throws ArgumentError parent(set(da2, [9 9; 9 9])) == [9 9; 9 9]
end

@testset "DimStack fields" begin
    @test_throws ArgumentError set(s, (x=a2, y=3a2))
    @test_throws ArgumentError set(s, (test2=a2, test3=hcat(a2, a2)))
    s2 = set(s, (test2=a2, test3=3a2))
    @test keys(s2) == (:test2, :test3)
    @test values(s2) == (a2, 3a2)
    @test metadata(set(s, Metadata(Dict(:testa => "test")))).val == Dict(:testa => "test")
end

@testset "DimStack Dimension" begin
    @test typeof(dims(set(s, row=X, column=Z))) <: Tuple{<:X,<:Z}
    @test typeof(dims(set(s, row=X(), column=Z()))) <: Tuple{<:X,<:Z}
    @test typeof(dims(set(s, row=:row2, column=:column2))) <: Tuple{<:Dim{:row2},<:Dim{:column2}}
    @test typeof(dims(set(s, :column => Ti(), :row => Z))) <: Tuple{<:Z,<:Ti}
    @test typeof(dims(set(s, Dim{:row}(Y()), Dim{:column}(X())))) <: Tuple{<:Y,<:X}
    @test typeof(dims(set(s, (Dim{:row}(Y), Dim{:column}(X))))) <: Tuple{<:Y,<:X}
    @test index(set(s, Dim{:row}([:x, :y, :z])), :row) == [:x, :y, :z]
end

@testset "DimArray dim Dimension" begin
    @test typeof(dims(set(da, X=:a, Y=:b))) <: Tuple{<:Dim{:a},<:Dim{:b}}
    @test typeof(dims(set(da2, Dim{:row}(Y()), Dim{:column}(X())))) <: Tuple{<:Y,<:X}
    @test typeof(dims(set(da, X => Ti(), Y => Z()))) <: Tuple{<:Ti,<:Z}
    @test typeof(dims(set(da2, :column => Ti(), :row => Z()))) <: Tuple{<:Z,<:Ti}
    @test typeof(dims(set(da2, (Dim{:row}(Y), Dim{:column}(X))))) <: Tuple{<:Y,<:X}
    @test typeof(dims(set(da2, row=X, column=Z))) <: Tuple{<:X,<:Z}
    @test typeof(dims(set(da2, row=X(), column=Z()))) <: Tuple{<:X,<:Z}
    @test typeof(dims(set(da2, row=:row2, column=:column2))) <: Tuple{<:Dim{:row2},<:Dim{:column2}}
    @test index(set(da2, Dim{:row}([:x, :y, :z])), :row) == [:x, :y, :z]
end

@testset "Dimension index" begin
    @test index(set(da2, :column => [:a, :b, :c, :d], :row => 4:6)) ==
        (4:6, [:a, :b, :c, :d])
    @test index(set(da2, column=Val((:a, :b, :c, :d)), row=Val((4:6...,)))) ==
        ((4:6...,), (:a, :b, :c, :d))
    @test index(set(s, :column => 10:5:20, :row => 4:6)) == (4:6, 10:5:20)
    @test step.(span(set(da2, :column => 10:5:20, :row => 4:6))) == (1, 5)
end

@testset "dim mode" begin
    interval_da = set(da, X=Intervals(), Y=Intervals())
    @test mode(set(da2, :column => NoIndex(), :row => Sampled(sampling=Intervals(Center())))) ==
        (Sampled(Ordered(), Regular(10.0), Intervals(Center())), NoIndex())
    @test mode(set(da2, column=NoIndex())) ==
        (Sampled(Ordered(), Regular(10.0), Points()), NoIndex())
    @test mode(set(da2, :column => NoIndex(), :row => Sampled())) ==
        (Sampled(Ordered(), Regular(10.0), Points()), NoIndex())
    @test mode(set(da, X=NoIndex(), Y=Categorical())) == (NoIndex(), Categorical())

    @testset "span" begin
        # TODO: should this error? the span step doesn't match the index step
        @test span(set(da2, row=Irregular(10, 12), column=Regular(9.9))) ==
            (Irregular(10, 12), Regular(9.9))
        @test _set(Sampled(), AutoSpan()) == Sampled()
        @test _set(Sampled(), Irregular()) == Sampled(AutoOrder(), Irregular(), AutoSampling())
    end

    @testset "locus" begin
        @test_throws ArgumentError set(da2, (End(), Center()))
        @test locus(set(interval_da, X(End()), Y(Center()))) == (End(), Center())
        @test locus(set(interval_da, X=>End(), Y=>Center())) == (End(), Center())
        @test locus(set(interval_da, X=End, Y=Center)) == (End(), Center())
        @test locus(set(da, Y=Center())) == (Center(), Center())
        @test _set(Intervals(Center()), Start()) == Intervals(Start())
        @test _set(Intervals(Center()), AutoLocus()) == Intervals(Center())
        @test _set(Points(), Center()) == Points()
        @test_throws ArgumentError _set(Points(), Start())
        @test_throws ArgumentError _set(Points(), End())
    end

    @testset "sampling" begin
        @test sampling(interval_da) == (Intervals(), Intervals())
        @test sampling(set(da, (X(Intervals(End())), Y(Intervals(Start()))))) ==
            (Intervals(End()), Intervals(Start()))
        @test _set(Sampled(), AutoSampling()) == Sampled()
        @test _set(Sampled(), Intervals()) == Sampled(AutoOrder(), AutoSpan(), Intervals())
        @test _set(Points(), AutoSampling()) == Points()
        @test _set(AutoSampling(), Intervals()) == Intervals()
        @test _set(AutoSampling(), AutoSampling()) == AutoSampling()
    end

    @testset "order" begin
        uda = set(da, Y(Unordered()))
        @test order(uda) == (Ordered(), Unordered(ForwardRelation()))
        @test order(set(uda, Y=ReverseRelation())) == (Ordered(), Unordered(ReverseRelation()))
        @test order(ArrayOrder, set(da, X=ReverseArray)) == (ReverseArray(), ForwardArray())
        @test relation(set(da, Y(ReverseRelation())), Y) == ReverseRelation()
        @test order(set(X(), AutoOrder())) == AutoOrder()
        @test indexorder(set(X(; mode=Sampled(order=Ordered())), ReverseIndex())) == ReverseIndex()
        @test arrayorder(set(X(; mode=Sampled(order=Ordered())), ReverseArray())) == ReverseArray()
        @test arrayorder(set(X(; mode=Sampled(order=Ordered())), ReverseArray())) == ReverseArray()
        @test indexorder(_set(Ordered(), ReverseIndex())) == ReverseIndex()
        @test arrayorder(_set(Ordered(), ReverseArray())) == ReverseArray()
        @test relation(_set(Ordered(), ReverseRelation())) == ReverseRelation()
        @test relation(_set(Unordered(), ReverseRelation())) == ReverseRelation()
        @test_throws ArgumentError _set(Unordered(), ReverseIndex())
        @test_throws ArgumentError _set(Unordered(), ReverseArray())
        @test _set(Sampled(), Ordered()) == Sampled(Ordered(), AutoSpan(), AutoSampling())
        @test _set(Sampled(; order=Unordered()), AutoOrder()) == Sampled(Unordered(), AutoSpan(), AutoSampling())
    end
end

@testset "metadata" begin
    @test metadata(set(X(), Metadata(Dict(:a=>1, :b=>2)))).val == Dict(:a=>1, :b=>2)
    dax = set(da, X(Metadata(Dict(:a=>1, :b=>2))))
    @test metadata(dims(dax), X).val == Dict(:a=>1, :b=>2)
    @test metadata(dims(dax), Y).val == Dict(:meta => "Y")
    dax = set(da, X(; metadata=Metadata(Dict(:a=>1, :b=>2))))
    @test metadata(dims(dax, X)).val == Dict(:a=>1, :b=>2)
    dax = set(da2, row=Metadata(Dict(:a=>1, :b=>2)))
    @test metadata(dims(dax, :row)).val == Dict(:a=>1, :b=>2)
    dax = set(da2, column=Metadata(Dict(:a=>1, :b=>2)))
    @test metadata(dims(dax, :column)).val == Dict(:a=>1, :b=>2)
    @test metadata(set(da, NoMetadata())) == NoMetadata()
    @test metadata(set(da, NoMetadata)) == NoMetadata()
end

@testset "all dim fields" begin
    dax = set(da, X(20:-10:10; mode=Sampled(), metadata=Metadata(Dict(:a=>1, :b=>2))))
    x = dims(dax, X)
    @test val(x) == 20:-10:10
    @test order(x) == Ordered(ReverseIndex(), ForwardArray(), ForwardRelation())
    @test span(x) == Regular(-10)
    @test mode(x) == Sampled(Ordered(ReverseIndex(), ForwardArray(), ForwardRelation()), Regular(-10), Points())
    @test metadata(x).val == Dict(:a=>1, :b=>2)
end

@testset "errors with set" begin
    @test_throws ArgumentError set(da, Sampled())
    @test_throws ArgumentError set(da, X=7)
    @test_throws ArgumentError _set(dims(da, X), X(7))
    @test_throws ArgumentError set(da, notadimname=Sampled())
end

@testset "_set nothing" begin
    @test _set(nothing, nothing) == nothing
    @test _set(1, nothing) == 1
    @test _set(nothing, 2) == 2
end


text:
function JacobiDG3(DG,F,z,Topo,Param)
ksi=DG.xw;
eta=DG.xw;
zeta=DG.xwZ;
z1=z[1];
z2=z[2];
n=DG.OrdPoly+1;
nz=DG.OrdPolyZ+1;
X=zeros(n,n,nz,3);
dXdx=zeros(n,n,nz,3,3);
dXdxI=zeros(n,n,nz,3,3);
J=zeros(n,n,nz);
theta=zeros(n,n);
for j=1:n
  for i=1:n
    for k=1:nz
    X[i,j,k,1:2]=0.25*((1-ksi[i])*(1-eta[j])*F.P[1:2,1]+
                    (1+ksi[i])*(1-eta[j])*F.P[1:2,2]+
                    (1+ksi[i])*(1+eta[j])*F.P[1:2,3]+
                    (1-ksi[i])*(1+eta[j])*F.P[1:2,4]);
    z=0.5*((1-zeta[k])*z1+(1+zeta[k])*z2);
    (X[i,j,k,3],dXdx[i,j,k,3,3])=Topo(X[i,j,k,1],X[i,j,k,2],z,Param);
    dXdx[i,j,k,3,3]=0.5*dXdx[i,j,k,3,3]*(z2-z1);
    end
  end
end

for k=1:nz
  dXdx[:,:,k,1,1]=DG.DS*X[:,:,k,1];
  dXdx[:,:,k,2,1]=DG.DS*X[:,:,k,2];
  dXdx[:,:,k,3,1]=DG.DS*X[:,:,k,3];
  dXdx[:,:,k,1,2]=reshape(X[:,:,k,1],n,n)*DG.DS';
  dXdx[:,:,k,2,2]=reshape(X[:,:,k,2],n,n)*DG.DS';
  dXdx[:,:,k,3,2]=reshape(X[:,:,k,3],n,n)*DG.DS';
end
for j=1:n
  for i=1:n
    for k=1:nz
      J[i,j,k]=det(reshape(dXdx[i,j,k,:,:],3,3));
      dXdxI[i,j,k,:,:]=inv(reshape(dXdx[i,j,k,:,:],3,3))*J[i,j,k];
    end
  end
end
return X,J,dXdx,dXdxI,theta
end


text:
## Import imbalance functions
include("Imbalance.jl")

srand(123456)

## Calculate the imbalance and probability of imbalance for given group/class vectors
datafile = ARGS[1]
numsamps = parse(Int, ARGS[2])
cc = ARGS[3]
outdir = ARGS[4]

## ensure trailing slash
if outdir[end] != '/'
    outdir = outdir * "/"
end

## read in the class/group vectors
data = readdlm(datafile, Any)
classvect = data[:,1]
groupvect = data[:,2]
## get matrix size
numspp = length(groupvect)
## Sanity check
if numspp != length(classvect)
    error("Group and class vectors differ in length")
end
## for comparison to earlier analysis, reduce to two classes (cc and not-cc)
if cc != "all"
    classvect[find(xx -> xx != cc, classvect)] = "Not-" * cc
    numclasses = 2
else
    numclasses = length(unique(classvect))
end
## count the number of species in each class
sppcounts = zeros(Int16, numclasses)
for ii in classvect
    sppcounts[findfirst(xx -> xx == ii, unique(classvect))] += 1
end
## get group dimensions/distribution
## (note that the groupings are always read as an integer vector starting at 0)
groupsizes = zeros(Int16, maximum(groupvect) + 1)
for ii in groupvect
    groupsizes[ii+1] += 1
end
## assemble empirical distribution
empirical = [zeros(Int16, maximum(groupvect) + 1) for xx in 1:numclasses]
for ii in unique(classvect)
    for jj in unique(groupvect)
        empirical[findfirst(xx -> xx == ii, unique(classvect))][jj+1] =
            sum(groupvect[classvect .== ii] .== jj)
    end
end
## Sanity checks
if sum(groupsizes) != sum(sppcounts) != numspp
    error("Group sizes do not sum to total number of species")
end
if size(hcat(empirical...)) != (length(groupsizes), length(sppcounts))
    error("Empirical observation does not conform with number of groups, classes")
end
if sum(empirical) != groupsizes
    error("Impossible empirical distribution given group sizes")
end
if [sum(xx) for xx in empirical] != sppcounts
    error("Impossible empirical distribution given species counts per class")
end
## Find pvalue of empirical observation
if numsamps > 0
    psi = Imbalance(empirical...)
    ## Repeatedly, uniformly sample potential distributions
    outfile = open(outdir * "$(rsplit(basename(datafile), ".", limit=2)[1])-$cc", "w")
    SampleDists(cc, psi, groupvect, classvect, numsamps, outfile)
    println(outdir * "$(rsplit(basename(datafile), ".", limit=2)[1])-$cc")
    close(outfile)
else
    presencevects = [GetPotentialPresences(groupsizes, pp) for pp in sppcounts]
    fullpresences = GetCompatiblePresences(groupsizes, presencevects...)
    imbatable = GetAllImbalances(fullpresences)
    psi = Imbalance(empirical...)
    prob = Emppvalue(empirical, imbatable)
    println("ψ = $psi, p(ψ) = $prob")
    outfile = open(outdir * "$(rsplit(basename(datafile), ".", limit=2)[1])", "w")
    write(outfile, "$cc,$psi,$prob\n")
    close(outfile)
end


text:
export frame3dd,
       frame3DD_circular_tube_truss_bar_family

Base.@kwdef struct Frame3DDBackend{K,T} <: LazyBackend{K,T}
  realized::Parameter{Bool}=Parameter(false)
  truss_nodes::Vector{<:TrussNode}=TrussNode[]
  truss_bars::Vector{<:TrussBar}=TrussBar[]
  shapes::Vector{<:Shape}=Shape[] # This contains all the rest that is not treated yet
  truss_node_data::Vector{TrussNodeData}=TrussNodeData[]
  truss_bar_data::Vector{TrussBarData}=TrussBarData[]
end

abstract type FR3DDKey end
const FR3DDId = Any
const FR3DDIds = Vector{FR3DDId}
const FR3DDRef = GenericRef{FR3DDKey, FR3DDId}
const FR3DDRefs = Vector{FR3DDRef}
const FR3DDEmptyRef = EmptyRef{FR3DDKey, FR3DDId}
const FR3DDUniversalRef = UniversalRef{FR3DDKey, FR3DDId}
const FR3DDNativeRef = NativeRef{FR3DDKey, FR3DDId}
const FR3DDUnionRef = UnionRef{FR3DDKey, FR3DDId}
const FR3DDSubtractionRef = SubtractionRef{FR3DDKey, FR3DDId}
const FR3DD = Frame3DDBackend{FR3DDKey, FR3DDId}

void_ref(b::FR3DD) = FR3DDNativeRef(-1)

const frame3dd = FR3DD()

# Frame3DD needs to merge nodes and bars
save_shape!(b::FR3DD, s::TrussNode) = maybe_merged_node(b, s)
save_shape!(b::FR3DD, s::TrussBar) = maybe_merged_bar(b, s)

# Frame3DD does not need layers
with_family_in_layer(f::Function, backend::FR3DD, family::Family) = f()

# Frame3DD Families
abstract type Frame3DDFamily <: BackendFamily{Any} end

Base.@kwdef struct Frame3DDTrussBarFamily <: Frame3DDFamily
  Ax  # cross section area
  Asy # shear area y-direction
  Asz # shear area z-direction
  Jxx # torsional moment of inertia - x axis
  Iyy # bending moment of inertia - y axis
  Izz # bending moment of inertia - z axis
  E   # elastic modulus
  G   # shear   modulus
  p   # roll angle
  d   # mass density
end

truss_bar_family_cross_section_area(f::Frame3DDTrussBarFamily) = f.Ax

#=
Circular Tube (outer radius= Ro, inner radius = Ri):

Ax = π ( Ro2 - Ri2 )
Asy = Asz = Ax / ( 0.54414 + 2.97294(Ri/Ro) - 1.51899(Ri/Ro)2 ) ± 0.05%
Jxx = (1/2) π ( Ro4 - Ri4 )
Ixx = Iyy = (1/4) π ( Ro4 - Ri4 )
=#

frame3DD_circular_tube_truss_bar_family(rₒ, rᵢ; E, G, p, d) =
  let Ax = annulus_area(rₒ, rᵢ),
      Asyz = Ax/(0.54414 + 2.97294*(rᵢ/rₒ) - 1.51899*(rᵢ/rₒ)^2),
      Jxx = π/2*(rₒ^4 - rᵢ^4),
      Ixxyy = Jxx/2
    Frame3DDTrussBarFamily(
      Ax  =Ax,
      Asy =Asyz,
      Asz =Asyz,
      Jxx =Jxx,
      Iyy =Ixxyy,
      Izz =Ixxyy,
      E   =E,
      G   =G,
      p   =p,
      d   =d)
 end

#=
 @rad1 = 25.0                   # Section diameter (mm)
 @thk1 = 1.6                    # Section thickness (mm)

 @section = "CHS_25x1.6"        # Tube name
 @material = "Steel C450"       # Duragal DualGrade RHS Grade C450Lo (AS 1163)
 @ftype = "Tube"                # Tube type (Tube, Beam, Channel, Angle)
 @emod = "210.0"                # Young's Modulus (GN per square metre) for steel
 @gmod = "79.0"                 # Modulus of Rigidity (GN per square metre) for steel
 @pang = "0.0"                  # section axis rotation angle (deg)
 @ax = "568.3"                  # cross-sectional area (mm^2)
 @asy = "3.2"                   # 2 x thickness (mm) shear area
 @asz = "3.2"                   # 2 x thickness (mm) shear area
 @jxx = "0.0585E+06"            # J (mm^4)
 @iyy = "0.0702E+06"            # Iyy (mm^4)
 @izz = "0.0237E+06"            # Izz (mm^4)
=#

set_backend_family(
   default_truss_bar_family(),
   frame3dd,
   frame3DD_circular_tube_truss_bar_family(0.0213/2, 0.0213/2-0.0026,
     E=210000000000.0, # (Young's modulus)
     G=81000000000.0, # (Kirchoff's or Shear modulus)
     p=0.0, # Roll angle
     d=77010.0)) # Density

#
delete_all_shapes(b::FR3DD) =
  begin
    empty!(b.truss_nodes)
    empty!(b.truss_bars)
    b.realized(false)
    b
  end

realize(b::FR3DD, s::TrussNode) =
  error("BUM")

realize(b::FR3DD, s::TrussBar) =
  error("BUM")

realize(b::FR3DD, s::Panel) =
  error("BUM")

#new_truss_analysis(v=nothing; self_weight=false, backend=frame3dd) =
displacements_from_frame3dd(b::FR3DD, filename, load) =
  let nodes = process_nodes(b.truss_nodes, load),
      bars = process_bars(b.truss_bars, nodes),
      supports = unique(filter(s -> s != false, map(n -> n.family.support, nodes))),
      loaded_nodes = filter(!truss_node_is_supported, nodes),
      supported_nodes = filter(truss_node_is_supported, nodes),
      shear = 0,		     # 1: include shear deformation effects, 0: don't
      geom  = 0,		     # 1: include geometric stiffness effects, 0: don't
      exagg = 1,		     # exaggeration factor for Gnuplot output
      scale = 1.0,		   # zoom scale factor for Gnuplot plotting
      dx = 1.0		       # x-axis increment for internal force calc's
    empty!(b.truss_node_data)
    append!(b.truss_node_data, nodes)
    empty!(b.truss_bar_data)
    append!(b.truss_bar_data, bars)
    open(filename, "w") do io
      println(io, "Frame analysis for Khepri");
      println(io, "%% node data ...");
      println(io, "$(length(nodes))\t\t%% number of nodes");
      println(io, "%% J\t\tX\t\tY\t\tZ\t\tr\t\tnode coordinates");
      for n in nodes
        i = n.id
        p = n.loc
        println(io, "$(i)\t$(p.x)\t$(p.y)\t$(p.z)\t0")
      end
      println(io, "");
      println(io, "%% reaction data ...");
      #nR = count(has_support, nodes);
      println(io, "$(length(supported_nodes))    %% number of nodes with reaction forces");
      println(io, "%% j\tRx\tRy\tRz\tRxx\tRyy\tRzz");
      for n in supported_nodes
        s = n.family.support
        print(io, "$(n.id),\t$(Int(s.ux)),\t$(Int(s.uy)),\t$(Int(s.uz))\t")
        println(io, "$(Int(s.rx)),\t$(Int(s.ry)),\t$(Int(s.rz))")
      end
      println(io, "")
      println(io, "%% frame element section property data ...")
      println(io, "$(length(bars))\t\t%% number of frame elements")
      println(io, "%% m\tn1\tn2\t\tAx\t\tAsy\t\tAsz\t\tJxx\t\tIyy\t\tIzz\t\tE\t\tG\t\tp\tdensity");
      for bar in bars
        f = family_ref(b, bar.family)
        print(io, "$(bar.id),\t$(bar.node1.id),\t$(bar.node2.id),")
        print(io, "\t$(f.Ax),\t$(f.Asy),\t$(f.Asz),")
        print(io, "\t$(f.Jxx),\t$(f.Iyy),\t$(f.Izz),")
        println(io, "\t$(f.E),\t$(f.G),\t$(f.p),\t$(f.d)")
      end
      println(io, "\n");
      println(io, "$(shear)\t\t%% 1: include shear deformation, 0: do not");
      println(io, "$(geom)\t\t%% 1: include geometric stiffness, 0: do not");
      println(io, "$(exagg)\t%% exagerate deformations in plotting");
      println(io, "$(scale)\t%% zoom scale factor for 3D plotting");
      println(io, "$(dx)\t%% x-axis increment for internal forces calc");
      println(io, "\n");
      println(io, "%% static load data ...");
      println(io, "1\t\t%% number of static load cases");
      println(io, "\t\t%% begin static load case 1 of 1");
      println(io, "%% gravitational acceleration for self-weight loading");
      println(io, "%% gX         gY         gZ");
      println(io, "  0.0        0.0        0.0\n");
      println(io, "$(length(loaded_nodes))\t\t%% number of loaded nodes");
      println(io, "%% j\t\tFx\t\tFy\t\tFz\t\tMxx\t\tMyy\t\tMzz");
      for n in loaded_nodes
        println(io, "$(n.id),\t$(load.x),\t$(load.y),\t$(load.z),\t0,\t0,\t0")
      end
      println(io, "\n");
      println(io, "0\t\t%% number of members with uniform distributed loads")
      println(io, "0\t\t%% number of members with trapezoidal loads");
      println(io, "0\t\t%% number of members with internal point loads");
      println(io, "0\t\t%% number of members with temperature loads");
      println(io, "0\t\t%% number of nodes with prescribed displacements");
      println(io, "\t\t%% end   static load case 1 of 1\n");
      println(io, "%% inertial load data ...");
      println(io, "0\t\t%% number of dynamic modes to analyze");
    end
  end

frame3dd_plugin = joinpath(dirname(dirname(abspath(@__FILE__))), "Plugins", "Frame3DD", "frame3dd.exe")

frame3dd_simulation_path() =
  mktempdir(tempdir(), prefix="Frame3DD_")

##########################################
backend_truss_analysis(b::FR3DD, load::Vec) =
  # Ensure extension is FMM to force Matlab mode
  let simulation_folder = frame3dd_simulation_path(),
#      input_path = joinpath(simulation_folder, "IOdata.FMM"),
#      output_path = joinpath(simulation_folder, "IOdata.OUT") ,
#      output_m_path = joinpath(simulation_folder, "IOdata_out.m")
       input_path = joinpath(simulation_folder, "IOdata.IN"),
       output_path = joinpath(simulation_folder, "IOdata.OUT")
    @info input_path
    displacements_from_frame3dd(b, input_path, load)
    try
      #run(`$frame3dd_plugin -w -i $input_path -o $output_path`)
      #run(`$frame3dd_plugin -i $input_path -o $output_path`)
      withenv("FRAME3DD_OUTDIR"=>simulation_folder) do
        run(`$frame3dd_plugin -q -i $input_path -o $output_path`)
      end
    catch e
    end
    lines = readlines(output_path)
    idx1 = findfirst(l->startswith(l, "N O D E   D I S P L A C E M E N T S"), lines) + 2
    idx2 = findnext(l->startswith(l, "F R A M E   E L E M E N T   E N D   F O R C E S"), lines, idx1)
    Dict([let strs = split(line)
            parse(Int, strs[1])=>vxyz(parse.(Float64, strs[2:4])...)
          end
          for line in lines[idx1:idx2-1]])
  end

node_displacement_function(b::FR3DD, results) =
  n -> get(results, n.id, vx(0))


text:
module SSLTestBasics
using SimpleSDMLayers
using Test

M = rand(Bool, (3,5))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)

@test longitudes(S) == range(0.1, 0.9; length=5)
@test latitudes(S) == range(1/6, 5/6; length=3)

M = rand(Bool, (4,3))
S = SimpleSDMPredictor(M, 0.2, 1.8, -1.0, 2.0)

@test longitudes(S) == range(S.left+stride(S,1), S.right-stride(S,1); length=size(S,2))
@test latitudes(S) == range(S.bottom+stride(S,2), S.top-stride(S,2); length=size(S,1))

bbox = boundingbox(S)
@test bbox.left == S.left
@test bbox.right == S.right
@test bbox.bottom == S.bottom
@test bbox.top == S.top

A = grid(S)
@test M == A

end


text:
module Bioceananigans

export LightGrowth!
export MixedLayerDepth!

include("LightGrowth.jl")
include("MixedLayerDepth.jl")

end # end module

max_stars_repo_name: iuryt/biocenanigans

max_issues_repo_name: iuryt/biocenanigans


text:
@testset "Parametric dimension" begin
  @test paramdim(Segment) == 1
  @test paramdim(Triangle) == 2
  @test paramdim(Quadrangle) == 2
  @test paramdim(PolyArea) == 2
  @test paramdim(Hexahedron) == 3
  @test paramdim(Pyramid) == 3
  @test paramdim(Tetrahedron) == 3
end


text:
"""
    abstract type GridapType end
"""
abstract type GridapType end

function Base.show(io::IO,object::GridapType)
  print(io,"$(nameof(typeof(object)))()")
end

"""
"""
function operate(op,a)
  s = "Unary operation $op is not defined for objects of type $(typeof(a))\n"
  s *= "Possible fix, define:\n"
  s *= "Gridap.Helpers.operate(::$(typeof(op)),::$(typeof(a)))"
  @unreachable s
end

"""
"""
function operate(op,a,b)
  s = "Binary operation $op is not defined for objects of type $(typeof(a)) and $(typeof(b))\n"
  s *= "Possible fix, define:\n"
  s *= "Gridap.Helpers.operate(::$(typeof(op)),::$(typeof(a)),::$(typeof(b)))"
  @unreachable s
end

for op in (:+,:-,:tr, :transpose, :adjoint)
  @eval begin
    function ($op)(a::GridapType)
      operate($op,a)
    end
  end
end

for op in (:+,:-,:*,:cross,:dot,:/)
  @eval begin

    function ($op)(a::GridapType,b::GridapType)
      operate($op,a,b)
    end

    function ($op)(a::GridapType,b)
      operate($op,a,b)
    end

    function ($op)(a,b::GridapType)
      operate($op,a,b)
    end

  end
end


text:
abstract type AbstractAxis{T, BL, BR, I} <: AbstractVector{T} end
abstract type AbstractDiscreteAxis{T, BL, BR, I} <: AbstractAxis{T, BL, BR, I} end

"""
    DiscreteAxis{T, BL, BR} <: AbstractAxis{T, BL, BR}

* T: Type of ticks
* BL, BR ∈ {:periodic, :reflecting, :infinite, :r0, :fixed}
* BL: left boundary condition
* BR: right boundary condition
* I: IntervalSets.Interval (closed or open boundaries)
"""
struct DiscreteAxis{T, BL, BR, I} <: AbstractAxis{T, BL, BR, I}
    interval::I
    ticks::Vector{T}
end
@inline size(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR} = size(dx.ticks)
@inline IndexStyle(::Type{<:DiscreteAxis}) = IndexLinear()
@inline length(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR} = length(dx.ticks)
@inline getindex(dx::DiscreteAxis{T, BL, BR}, i::Int) where {T, BL, BR} = dx.ticks[i]
@inline setindex!(dx::DiscreteAxis{T, BL, BR}, v::T, i::Int) where {T, BL, BR} = setindex!(dx.ticks, v, i)
@inline axes(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR} = axes(dx.ticks)

function DiscreteAxis{T, BL, BR}(int::I, ticks::Vector{T})::DiscreteAxis{T, BL, BR, typeof(int)} where {T, BL, BR, I}
    return DiscreteAxis{T, BL, BR, typeof(int)}( int, ticks )
end

"""
    DiscreteAxis(left_endpoint::T, right_endpoint::T, BL::Symbol, BR::Symbol, L::Symbol, R::Symbol, ticks::AbstractVector{T}) where {T}

* T: Type of ticks
* BL, BR ∈ {:periodic, :reflecting, :infinite, :r0, :fixed}
* L, R {:closed, :open}
* ticks: Ticks of the axis
"""
function DiscreteAxis(left_endpoint::T, right_endpoint::T, BL::Symbol, BR::Symbol, L::Symbol, R::Symbol, ticks::AbstractVector{T}) where {T}
    int::Interval{L, R, T} = Interval{L, R, T}( left_endpoint, right_endpoint )
    return DiscreteAxis{T, BL, BR, typeof(int)}( int, ticks )
end

function sizeof(dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    return sizeof(dx.interval) + sizeof(dx.ticks)
end

function print(io::IO, dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    print(io, dx.interval, " - length = ", length(dx))
end
function println(io::IO, dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    println(io, dx.interval)
    println(io, "length = ", length(dx))
end
function show(io::IO, dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    println(io, dx)
end
function show(io::IO, ::MIME"text/plain", dx::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    show(io, dx)
end

function get_boundary_types(int::Interval{L, R})::Tuple{Symbol, Symbol} where {L, R}
    return L, R
end

function get_boundary_types(ax::DiscreteAxis{T,LB,RB})::NTuple{4, Symbol} where {T, LB, RB}
    return LB, RB, get_boundary_types(ax.interval)...
end

function uniq(v::Vector{T})::Vector{T} where {T <: Real}
    v1::Vector{T} = Vector{T}()
    if length(v) > 0
        laste::T = v[1]
        push!(v1, laste)
        for e in v
            if e != laste
                laste = e
                push!(v1, laste)
            end
        end
    end
    return v1
end

function merge_axis_ticks_with_important_ticks(ax::DiscreteAxis{T}, impticks::Vector{T}; atol::Real = 0.0001 )::Vector{T} where {T}
    v::Vector{T} = T[]
    for r in impticks if in(r, ax.interval) push!(v, r) end end
    for r in ax push!(v, r) end
    sort!(v)
    v = uniq(v)
    delete_idcs::Vector{Int} = Int[]
    for i in 1:(length(v) - 1)
        if (v[i + 1] - v[i]) < atol
            if !in(v[i], impticks) push!(delete_idcs, i) end
            if !in(v[i + 1], impticks) push!(delete_idcs, i + 1) end
        end
    end
    delete_idcs = sort(uniq(delete_idcs))
    deleteat!(v, delete_idcs)
    for impv in impticks
        if !in(impv, v) && in(impv, ax.interval)
            error("Important ticks were removed.")
        end
    end
    return v
end

function range(interval::Interval{:closed, :closed, T}; step::Union{Missing, T} = missing, length::Union{Missing, Int} = missing) where {T}
    stop::T = interval.right
    if ismissing(step) && ismissing(length)
        range(interval.left, stop = stop)
    elseif ismissing(step)
        range(interval.left, stop = stop, length=length)
    elseif ismissing(length)
        range(interval.left, stop = stop, step=step)
    else
        error(KeyError, ": Both keywords `step` and `length` were given. But only one is allowed.")
    end
end

function range(interval::Interval{:closed, :open, T}; step::Union{Missing, T} = missing, length::Union{Missing, Int} = missing) where {T}
    if ismissing(step) && ismissing(length)
        length::Int = 2
        stop::T = (interval.right + interval.left) / 2
        range(interval.left, stop = stop, length=2)
    elseif ismissing(step)
        stop = interval.right - ( interval.right - interval.left ) / length
        range(interval.left, stop = stop, length=length)
    elseif ismissing(length)
        # stop = interval.right - interval.right % step
        stop = geom_round(interval.right - step)
        range(interval.left, stop = stop, step=step)
    else
        error(KeyError, ": Both keywords `step` and `length` were given. But only one is allowed.")
    end
end

function range(interval::Interval{:open, :closed, T}; step::Union{Missing, T} = missing, length::Union{Missing, Int} = missing) where {T}
    stop::T = interval.right
    if ismissing(step) && ismissing(length)
        step::T = (stop - interval.left) / 2
        range(interval.left + step, stop = stop, length=2)
    elseif ismissing(step)
        step = (stop - interval.left) / length
        range(interval.left + step, stop = stop, length=length)
    elseif ismissing(length)
        range(interval.left + step, stop = stop, step=step)
    else
        error(KeyError, ": Both keywords `step` and `length` were given. But only one is allowed.")
    end
end

function range(interval::Interval{:open, :open, T}; step::Union{Missing, T} = missing, length::Union{Missing, Int} = missing) where {T}
    if ismissing(step) && ismissing(length)
        step::T = ( interval.right - interval.left ) / 3
        range(interval.left + step, stop = interval.right - step, length=2)
    elseif ismissing(step)
        step = ( interval.right - interval.left ) / (length + 1)
        range(interval.left + step, stop = interval.right - step, length=length)
    elseif ismissing(length)
        tmp::T = interval.right % step
        if tmp == 0 tmp = step end
        stop = interval.right - tmp
        range(interval.left + step, stop = stop, step=step)
    else
        error(KeyError, ": Both keywords `step` and `length` were given. But only one is allowed.")
    end
end

function DiscreteAxis{BL, BR}(interval::Interval{L, R, T}; step::Union{Missing, T} = missing, length::Union{Missing, Int} = missing)::DiscreteAxis{T, BL, BR} where {L, R, T, BL, BR}
    ticks::Vector{T} = collect(range(interval, step=step, length=length))
    if T == Float32 || T == Float64
        ticks = round.(ticks, sigdigits = geom_sigdigits(T))
        for iv in eachindex(ticks)
            if isapprox(ticks[iv], 0, atol = geom_atol_zero(T))
                ticks[iv] = zero(T)
            end
        end
    end
    DiscreteAxis{T, BL, BR}(interval, ticks)
end

function midpoints(a::Vector{T})::Vector{T} where {T}
    @inbounds r::Vector{T} = a[1:end-1]
    @simd for i in eachindex(r)
        @inbounds r[i] += 0.5 * (a[i + 1] - a[i])
    end
    return r
end

function get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :fixed, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :fixed} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :periodic, :periodic} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :infinite, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[1] = ticks_ext[2] - Δ
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :r0, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    ticks_ext[1] = ticks_ext[2] - (ticks_ext[3] - ticks_ext[2])
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :r0, :fixed} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    ticks_ext[1] = ticks_ext[2] - (ticks_ext[3] - ticks_ext[2])
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :r0, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    ticks_ext[1] = ticks_ext[2] - (ticks_ext[3] - ticks_ext[2])
    Δ::T = ticks_ext[end-1] - ticks_ext[end - 2]
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :fixed, :fixed} )::Vector{T} where {T}
    # same as get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :infinite, :fixed} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[1] = ticks_ext[2] - Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :infinite, :reflecting} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[1] = ticks_ext[2] - Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :fixed, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end
function get_extended_ticks( ax::DiscreteAxis{T, :reflecting, :infinite} )::Vector{T} where {T}
    ticks_ext::Vector{T} = Array{T}(undef, length(ax.ticks) + 2)
    ticks_ext[2:end-1] = ax.ticks
    set_periodic_bondary_ticks!(ticks_ext, ax.interval)
    Δ::T = 1 * (ticks_ext[end-1] - ticks_ext[2])
    ticks_ext[end] = ticks_ext[end - 1] + Δ
    return ticks_ext
end

function set_periodic_bondary_ticks!( ticks::Vector{T}, interval::Interval{:closed, :open, T})::Nothing where {T}
    ticks[1] = ticks[2] - (interval.right - ticks[end - 1])
    ticks[end] = interval.right
    nothing
end
function set_periodic_bondary_ticks!( ticks::Vector{T}, interval::Interval{:open, :closed, T})::Nothing where {T}
    ticks[1] = interval.left
    ticks[end] = ticks[end - 1] + (ticks[2] - interval.left)
    nothing
end
function set_periodic_bondary_ticks!( ticks::Vector{T}, interval::Interval{:open, :open, T})::Nothing where {T}
    ticks[1] = interval.left
    ticks[end] = interval.right
    nothing
end

function set_periodic_bondary_ticks!( ticks::Vector{T}, interval::Interval{:closed, :closed, T})::Nothing where {T, ispolaraxis}
    if length(ticks) == 3
        ticks[1] = ticks[2] - 2π
        ticks[end] = ticks[2] + 2π # -> Δmidpoint_φ = 2π -> area of circle is 2π * 0.5*r^2
    else
        ticks[1] = ticks[2] - (ticks[3] - ticks[2])
        ticks[end] = ticks[end - 1] + (ticks[end - 1] - ticks[end - 2])
    end
    nothing
end

function searchsortednearest(a::AbstractVector{T}, x::T)::Int where {T <: Real}
    idx::Int = searchsortedfirst(a, x)
    if (idx == 1) return idx end
    if (idx > length(a)) return length(a) end
    if (a[idx] == x) return idx end
    if (abs(a[idx] - x) < abs(a[idx - 1] - x))
        return idx
    else
        return idx - 1
    end
end
@inline function searchsortednearest(ax::DiscreteAxis{T}, x::T)::Int where {T <: Real}
    return searchsortednearest(ax.ticks, x)
end
function searchsortednearest(ax::DiscreteAxis{T, :periodic, :periodic}, x::T)::Int where {T <: Real}
    if x in ax.interval
        return searchsortednearest(ax.ticks, x)
    else
        period::T = ax.interval.right - ax.interval.left
        v::T = x
        while v >= ax.interval.right
            v -= period
        end
        while v < ax.interval.left
            v += period
        end
        return searchsortednearest(ax.ticks, v)
    end

end

function DiscreteAxis(nt::NamedTuple; unit = u"m/m")
    T = typeof(ustrip(nt.knots[1]))
    knots::Vector{T} = convert(Vector{T}, ustrip.(uconvert.(unit, nt.knots)))
    lep::T = ustrip(uconvert.(unit, nt.interval.left_boundary.endpoint ))
    rep::T = ustrip(uconvert.(unit, nt.interval.right_boundary.endpoint))
    int = Interval{nt.interval.left_boundary.closedopen, nt.interval.right_boundary.closedopen}( lep, rep )
    return DiscreteAxis{T, nt.interval.left_boundary.boundaryhandling, nt.interval.right_boundary.boundaryhandling, typeof(int)}(
        int, knots
    )
end

Base.convert(T::Type{DiscreteAxis}, x::NamedTuple; unit = u"m/m") = T(x)

function NamedTuple(ax::DiscreteAxis{T, BL, BR}; unit = u"m/m") where {T, BL, BR}
    int::Interval = ax.interval
    int_types::Tuple{Symbol, Symbol} = get_boundary_types(int)
    return (
        knots = ax.ticks * unit,
        interval = (
            left_boundary = (
                endpoint = int.left * unit,
                closedopen = int_types[1],
                boundaryhandling = BL,
            ),
            right_boundary = (
                endpoint = int.right * unit,
                closedopen = int_types[2],
                boundaryhandling = BR,
            ),
        )
    )
end

Base.convert(T::Type{NamedTuple}, x::DiscreteAxis; unit = u"m/m") = T(x)


text:
export shuffle_indices!, add_remove_kink_chain!

"""
    shuffle_indices!(m::Model, e::Ensemble, c::Configuration)
Performs an Update that takes a Kink and randomly shuffles its annihlators with each other and its creators.
The Sets of creators and annihilators of the kinks stay the same the acceptance propability is therefore always 1.
"""
function shuffle_indices!(m::Model, e::Ensemble, c::Configuration)
    if isempty(c.kinks)
        return 1.0, Step()
    end
    kink = rand(c.kinks)
    if rand() > 0.5
        # shuffle creators
        Δ = Step(nothing, (kink,), nothing, (first(kink) => T4(last(kink).j,last(kink).i,last(kink).k,last(kink).l),))
    else
        # shuffle annihilators
        Δ = Step(nothing, (kink,), nothing, (first(kink) => T4(last(kink).i,last(kink).j,last(kink).l,last(kink).k),))
    end
    apply_step!(c,Δ)
    return 1.0, Δ
end

""" perform a number of subsequent updates
    first perform a number of updates which add kinks
    second perform the same number of updates which remove kinks.
    For now it is not sure wether it is recommandable or not to use this update"""
function add_remove_kink_chain!(m::Model, e::Ensemble, c::Configuration)
    if isempty(c.kinks)
        return 1.0, Step()
    end
    add_single_kinks = [add_type_C!, add_type_D!, add_type_E!]#
    remove_single_kinks = [remove_type_C!, remove_type_D!, remove_type_E!]#

    chain_length = rand(1:1)
    acc_prob = 1.0
    step_list = Array{Step,1}()# TODO: it may be more efficient to pre-allocate this output and only set it here
    for i in  1:chain_length
        dv, Δ = rand(add_single_kinks)(m,e,c)
        acc_prob *= dv
        push!(step_list, Δ)
        if (iszero(acc_prob) || (Δ == Step()))
            apply_back_step!(c, step_list)
            return 0.0, Step()
        end
    end
    for i in  1:chain_length
        dv, Δ = rand(remove_single_kinks)(m,e,c)
        acc_prob *= dv
        push!(step_list, Δ)
        if (iszero(acc_prob) || (Δ == Step()))
            apply_back_step!(c, step_list)
            return 0.0, Step()
        end
    end
    return acc_prob, step_list
end

"""
    isuseful(c::Configuration, up::Function)
Function with the name isuseful checkt wether it is senseble to propose a spezific Update-Type.
If no explizit function is defined for the type spezific type of Update its called on, this function will be called, which just returns true.
"""
function isuseful(c::Configuration, up::Function)
    return true
end

function isuseful(c::Configuration, up::typeof(shuffle_indices!))
    if isempty(c.kinks)
        return false
    else
        return true
    end
end


text:
@with_kw struct Fiber{T, S}
   foo::Int = 6
   diameter::T#fiber diameter
   length::T# fiber length
   Tk::T# temperture
   p::T#pressure
   loss::Bool# attenuation Need modification to accept array
   plasma::Bool
   pltype::Symbol
   dispersion::Bool
   shock::Bool
   fr::Float64
   thg::Bool
   nsaves::Int64=1000
end

function systemInfo(pulse::Pulse, fiber::Fiber, gases::Array{Gas{Float64, Symbol},1})
    gasMix = Vector()
    for i in 1:length(gases)
        t = gases[i].type
        append!(gasMix, Dict("$t" => gases[i].pp))
    end
    println(gasMix)
    dict = Dict("FWHM" => pulse.FWHM,
    "wavelength" => pulse.λ0,
    "P0" => pulse.P0,
    "natural_pulse_width" => pulse.τ0,
    "energy" => pulse.energy,
    "gas" => gasMix,
    "pressure" => gases[1].p,
    "radius" => fiber.diameter/2,
    "plasma" => fiber.plasma)
    if fiber.plasma == true
        dict["plasma"] = fiber.pltype
    else
        dict["plasma"] = false
    end
    return JSON.json(dict)
end

function systemInfo(pulse::Pulse,gases::Array{Gas{Float64, Symbol},1})
    gasMix = ""
    for i in 1:length(gases)
        t = string(gases[i].type)*"_"*string(gases[i].pp)*"_"
        gasMix *= t
    end

    return gasMix*string(Dates.minute(now())+Dates.hour(now())*100+Dates.second(now())*0.01)
end

function effectiveArea(diameter)
    return 1.5*(diameter/2)^2
end

function rNGas(p, Tk)
    T0 = 273.15
    p0 = 1.01325
    return p*T0/(p0*Tk)
end

function χ30(gas, p, Tk)
    #=
    Chi3 of a 'gas' [m^2/V^2], mostly from Lehmeier 1985
    pressure: P [bar]
    temperature: T [K]
    =#
    Chi3He = 3.43e-28
    Chi3N2 = 21.1*Chi3He
    Chi3Ne = 1.8*Chi3He
    Chi3Ar = 23.5*Chi3He
    Chi3Kr = 64.0*Chi3He
    Chi3Xe = 188.2*Chi3He
    Chi3O2 = 20.93*Chi3He #estimated from n2 from https://www.osapublishing.org/view_article.cfm?gotourl=https%3A%2F%2Fwww%2Eosapublishing%2Eorg%2FDirectPDFAccess%2FE77D9C07-DA73-08A6-1DBCC7FF6E9E255C_34926%2Fjosab-14-3-650%2Epdf%3Fda%3D1%26id%3D34926%26seq%3D0%26mobile%3Dno&org=Heriot-Watt%20University%20Library
    Chi3Air = Chi3O2*0.20946 + Chi3N2*(1-0.20946-0.009340) + Chi3Ar*0.009340 #estimated from N2 and O2

    if gas == :N2
        return rNGas(p, Tk)*4*Chi3N2
    elseif gas == :He
        return rNGas(p, Tk)*4*Chi3He
    elseif gas == :Ne
        return rNGas(p, Tk)*4*Chi3Ne
    elseif gas == :Ar
        return rNGas(p, Tk)*4*Chi3Ar
    elseif gas == :Kr
        return rNGas(p, Tk)*4*Chi3Kr
    elseif gas == :Xe
        return rNGas(p, Tk)*4*Chi3Xe
    elseif gas == :O2
        return rNGas(p, Tk)*4*Chi3O2
    elseif gas == :O3
        return rNGas(p, Tk)*4*Chi3O2*1.5
    elseif gas == :Air
        return rNGas(p, Tk)*4*Chi3Air
    end
end

function χ3(gas, p, Tk)
    sus = 0
    for i in keys(gas)
        sus += χ30(i, gas[i]*p, Tk)
    end
    return sus
end

function n_2(gas::Symbol,  p, Tk)
    res = 3.0*χ3(gas, p, Tk)/(4.0*ϵ0*c)
    return res
end

function n_2(χ3)
    res = 3.0*χ3/(4.0*ϵ0*c)
    return res
end

function getGamma(gas::Symbol,  p, Tk, ω0, diameter)
    return n_2(gas, p, Tk)*ω0/(c*effectiveArea(diameter))
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
function to_arrow_symbol(button_set)
    for b in button_set
        GLFW.KEY_RIGHT == b && return :right
        GLFW.KEY_LEFT  == b && return :left
        GLFW.KEY_DOWN  == b && return :down
        GLFW.KEY_UP    == b && return :up
    end
    return :nothing
end

function mousedragg_objectid(mouse_dragg, mouse_hover)
    map(mouse_dragg) do dragg
        value(mouse_hover), dragg
    end
end

function add_complex_signals!(screen)
    @materialize keyboard_buttons, mouse_buttons = screen.inputs
    no_scancode = map(remove_scancode, keyboard_buttons)
    button_s = merge(
        button_signals(no_scancode, :button),
        button_signals(mouse_buttons, :mouse_button)
    )
    #mousedragdiff_id = mousedragg_objectid(screen.inputs, mouse2id(screen))
    #selection        = foldp(drag2selectionrange, 0:0, mousedragdiff_id)
    arrow_navigation = const_lift(to_arrow_symbol, button_s[:buttons_pressed])
    merge!(
        screen.inputs,
        Dict{Symbol, Any}(
            #:mousedragg_objectid => mousedragdiff_id,
           # :selection           => selection,
            :arrow_navigation    => arrow_navigation
        ),
        button_s
    )
    screen
end

"""
Builds a Set of keys, accounting for released and pressed keys
"""
function currently_pressed_keys(v0::Set{Int}, button_action_mods)
    button, action, mods = button_action_mods
    if button != GLFW.KEY_UNKNOWN
        if action == GLFW.PRESS
            push!(v0, button)
        elseif action == GLFW.RELEASE
            delete!(v0, button)
        elseif action == GLFW.REPEAT
            # nothing needs to be done, besides returning the same set of keys
        else
            error("Unrecognized enum value for GLFW button press action: $action")
        end
    end
    return v0
end

function remove_scancode(button_scancode_action_mods)
    button, scancode, action, mods = button_scancode_action_mods
    button, action, mods
end
isreleased(button) = button[2] == GLFW.RELEASE
isdown(button) = button[2] == GLFW.PRESS

"""
Creates high level signals from the raw GLFW button signals.
Returns a dictionary with button released and down signals.
It also creates a signal, which is the set of all currently pressed buttons.
`name` is used to name the dictionary keys.
`buttons` is a tuple of (button, action, mods)::NTuple{3, Int}
"""
function button_signals(buttons::Signal{NTuple{3, Int}}, name::Symbol)
    keyset = Set{Int}()
    sizehint!(keyset, 10) # make it less suspicable to growing/shrinking
    released = filter(isreleased, buttons.value, buttons)
    down     = filter(isdown, buttons.value, buttons)
    Dict{Symbol, Any}(
        Symbol("$(name)_released") => map(first, released),
        Symbol("$(name)_down")     => map(first, down),
        Symbol("$(name)s_pressed") => foldp(
            currently_pressed_keys, keyset, buttons
        )
    )
end

"""
Selection of random objects on the screen is realized by rendering an
object id + plus an arbitrary index into the framebuffer.
The index can be used for e.g. instanced geometries.
"""
immutable SelectionID{T <: Integer} <: FixedVectorNoTuple{2, T}
    id::T
    index::T
    function SelectionID(args::NTuple{2, T})
        new{T}(args[1], args[2])
    end
end

begin
global push_selectionqueries!

const selection_data = Array(SelectionID{UInt16}, 1, 1)
const old_mouse_position = Array(Vec{2, Float64}, 1)

function push_selectionqueries!(screen)
    mouse_position   = value(mouseposition(screen))
    selection_signal = mouse2id(screen)
    window_size      = widths(screen)
    buff  = framebuffer(screen).objectid
    if old_mouse_position[] != mouse_position
        glReadBuffer(GL_COLOR_ATTACHMENT1)
        x,y = Vec{2, Int}(map(floor, mouse_position))
        w,h = window_size
        if x > 0 && y > 0 && x <= w && y <= h
            glReadPixels(x, y, 1, 1, buff.format, buff.pixeltype, selection_data)
            val = convert(Matrix{SelectionID{Int}}, selection_data)[1,1]
            push!(selection_signal, val)
        end
        old_mouse_position[] = mouse_position
    end
end

end

export screenshot, screenbuffer


text:
ν(G, L, Es, p) = p[1]*(G*L-Es/p[2])/(L + p[3]*G)

function rate(C)
    G, L, Es = C
    p = (9.03f-2, 5.0f-2, 1.9f1)
    ν(G, L, Es, p)
end

function reaction(du, u, p, t)
    G, L, Es = u
    dG  = du[1] = - ν(G, L, Es, p)
    dL  = du[2] = - ν(G, L, Es, p)
    dEs = du[3] =   ν(G, L, Es, p)
    return du
end

function reaction!(du, u, p, t)
    du = reaction(du, u, p, t)
end

function generate_true_solution(u₀, tₑ; nPoints=20)
    p = (9.03f-2, 5.0f-2, 1.9f1)
    tspan = (0.0, tₑ)

    prob = ODEProblem(reaction!, u₀, tspan, p, saveat=0:tₑ/nPoints:tₑ)
    sol = solve(prob, Tsit5())
    G = [c[1] for c in sol.u]
    L = [c[2] for c in sol.u]
    Es = [c[3] for c in sol.u]
    time = sol.t

  return permutedims([G L Es]), time
end

function plot_reaction(data)
  C, time = data

  pl = plot(xlabel="Time (min)", ylabel="Concentration (mM)",
      ylims=(-2, 1.2*maximum(C)),
      xlims=(-2, time[end]*1.01))
  scatter!(pl, time, C[1,:], label="Sucrose")
  scatter!(pl, time, C[2,:], label="Stearic acid")
  scatter!(pl, time, C[3,:], label="Sucr-Stea-Esther")
  return pl
end

function prediction_vs_data(data, pred)
    C, time = data

    pl = plot_reaction(data)
    lcs = pl.series_list
    plot!(pl, time, pred[1,:], label = "", lc=lcs[1][:linecolor])
    plot!(pl, time, pred[2,:], label = "",  lc=lcs[2][:linecolor])
    plot!(pl, time, pred[3,:], label = "",  lc=lcs[3][:linecolor])
    return pl
end


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.MVLinearGaussian
import SMCExamples.Particles.MVFloat64Particle
using StaticArrays
using Plots
Plots.gr()

setRNGs(0)

model, theta, ys, ko = MVLinearGaussian.defaultMVLGModel(2, 10)

nsamples = 2^14

smcio = SMCIO{model.particle, model.pScratch}(16, model.maxn, 1, true, 2.0)

v = Vector{MVFloat64Particle{2}}(undef, 10)
for p = 1:10
  v[p] = MVFloat64Particle{2}()
  v[p].x .= zeros(MVector{2, Float64})
end

function smoothingMeans(model, smcio, m, v)
  tmp = Vector{MVector{2,Float64}}(undef, 10)
  for p = 1:10
    tmp[p] = zeros(MVector{2, Float64})
  end
  for i = 1:m
    csmc!(model, smcio, v, v)
    for p = 1:10
      tmp[p] .+= v[p].x
    end
  end
  return tmp ./ m
end

@time out = smoothingMeans(model, smcio, nsamples, v)
plot((x -> x[1]).(out))
plot!((x -> x[1]).(ko.smoothingMeans), color=:red)

plot((x -> x[2]).(out))
plot!((x -> x[2]).(ko.smoothingMeans), color=:red)


text:
using Nextion
using Nextion: checkcommandcomplete, _write_end_of_command
using Test
using LibSerialPort: sp_blocking_read

struct NexResponse
    code::Return.Code.ReturnCode
end

@testset "raw_cmd_invalid" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    function wait_response(ser)
        delim = Char.([0xff, 0xff, 0xff])
        timeout_ms = 1000
        s = readuntil(ser, delim, timeout_ms)
        r = unsafe_wrap(Vector{UInt8}, s)
        NexResponse(Return.code(r[1]))
    end

    @testset "valid" begin
        pid = 0
        cmd = "page $pid"
        send(nexSerial, cmd)
        cmd = "t0.txt=\"Hello\""
        send(nexSerial, cmd)
    end

    @testset "invalid" begin
        pid = 0
        cmd = "page $pid"
        send(nexSerial, cmd)
        cmd = "t0.txt=\"Valid cmd\""
        send(nexSerial, cmd)
        cmd = "t0.txt=\"Invalid cmd"
        send(nexSerial, cmd)
        r = wait_response(nexSerial._serial)
        println(r)
#    cmd = "pageEEE 0"
    #    _code = send(nexSerial, cmd)
    #    checkcommandcomplete(nexSerial)
    #    #@test _code == Return.Code.INVALID_CMD
    end

    close(nexSerial)
    @test true

end


text:
using MendelKinship
using Test
using DataFrames
using CSV

datadir = joinpath(@__DIR__, "../data/documentation_data/")
cd(datadir)

@testset "Example 2 in Tutorial" begin
    Kinship("control_compare_29a.txt")
    @test isfile("kinship_file_output.txt")

    df = CSV.read("kinship_file_output.txt", DataFrame) # read result
    @test size(df) == (22578, 7)
    @test df[1, 1] == 14
    @test df[1, 2] == 14
    @test df[1, 3] == 26732
    @test df[1, 4] == 264
    @test df[1, 5] == 0.0
    @test isapprox(df[1, 6], 0.10955177290095937, atol=1e-7)
    @test isapprox(df[1, 7], 5.319415124407314, atol=1e-7)
end


text:
################################################################################
#
#  FakeFmpqMat.jl : Model fmpq_mat's as fmpz_mat's with denominator
#
################################################################################

iszero(x::FakeFmpqMat) = iszero(x.num)

issquare(x::FakeFmpqMat) = issquare(numerator(x))

function numerator(x::FakeFmpqMat; copy::Bool = true)
  if copy
    return deepcopy(x.num)
  else
    return x.num
  end
end

function denominator(x::FakeFmpqMat; copy::Bool = true)
  if copy
    return deepcopy(x.den)
  else
    return x.den
  end
end

nrows(x::FakeFmpqMat) = nrows(x.num)

ncols(x::FakeFmpqMat) = ncols(x.num)

function simplify_content!(x::FakeFmpqMat)
  c = content(x.num)
  gcd!(c, c, x.den)
  if !isone(c)
    divexact!(x.num, x.num, c)
    divexact!(x.den, x.den, c)
  end
  y = canonical_unit(x.den)
  if !isone(y)
    mul!(x.den, x.den, y)
    mul!(x.num, x.num, y)
  end
end

################################################################################
#
#  Getindex
#
################################################################################

Base.getindex(a::FakeFmpqMat, i::Int, j::Int) = fmpq(a.num[i, j], a.den)

################################################################################
#
#  Hashing
#
################################################################################

function hash(a::FakeFmpqMat, b::UInt)
  h = xor(Base.hash(a.num, b), Base.hash(a.den, b))
  h = xor(h, Base.hash(b))
  h = (h << 1) | (h >> (sizeof(Int)*8 - 1))
  return h
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FakeFmpqMat)
  print(io, "FakeFmpqMat with numerator\n", x.num, "\nand denominator ", x.den)
end

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::FakeFmpqMat)
  return FakeFmpqMat(-x.num, x.den, true)
end

#TODO: may be possible to simplify more efficiently.
#The content of the numerator of the inverse may be non trivial!
function inv(x::FakeFmpqMat)
  i, d_i = pseudo_inv(x.num)
  g = gcd(d_i, x.den)
  if isone(g)
    return FakeFmpqMat(i * x.den, d_i)
  elseif g == d_i
    return FakeFmpqMat(i * divexact(x.den, d_i))
  elseif g == x.den
    return FakeFmpqMat(i, divexact(d_i, x.den))
  else
    return FakeFmpqMat(i * divexact(x.den, g), divexact(d_i, g))
  end
end

################################################################################
#
#  Binary operations
#
################################################################################

function +(x::FakeFmpqMat, y::FakeFmpqMat)
  t = y.den * x.num + x.den * y.num
  d = x.den * y.den
  z = FakeFmpqMat(t, d)
  return z
end

function *(x::FakeFmpqMat, y::FakeFmpqMat)
  t = x.num * y.num
  d = x.den * y.den
  z = FakeFmpqMat(t, d)
  return z
end

function mul!(z::FakeFmpqMat, x::FakeFmpqMat, y::FakeFmpqMat)
  @req ncols(x) == nrows(y) "Wrong dimensions"
  @req nrows(z) == nrows(x) "Wrong dimensions"
  @req ncols(z) == ncols(y) "Wrong dimensions"
  mul!(z.num, x.num, y.num)
  mul!(z.den, x.den, y.den)
  simplify_content!(z)
  return z
end

################################################################################
#
#  Adhoc binary operations
#
################################################################################

function *(x::FakeFmpqMat, y::fmpz_mat)
  t = x.num * y
  z = FakeFmpqMat(t, denominator(x))
  return z
end

function *(x::fmpz_mat, y::FakeFmpqMat)
  t = x * y.num
  z = FakeFmpqMat(t, denominator(y))
  return z
end

function *(x::Integer, y::FakeFmpqMat)
  g = gcd(x, y.den)
  if isone(g)
    return FakeFmpqMat(y.num * x, y.den, true)
  else
    return FakeFmpqMat(y.num * divexact(x, g), divexact(y.den, g), true)
  end
end

function *(x::fmpz, y::FakeFmpqMat)
  g = gcd(x, y.den)
  if isone(g)
    return FakeFmpqMat(y.num * x, y.den, true)
  else
    return FakeFmpqMat(y.num * divexact(x, g), divexact(y.den, g), true)
  end
end

*(x::FakeFmpqMat, y::Integer) = y * x

*(x::FakeFmpqMat, y::fmpz) = y * x

function *(x::fmpq, y::FakeFmpqMat)
  n = numerator(x)
  d = denominator(x)
  g = gcd(n, y.den)
  if isone(g)
    return FakeFmpqMat(y.num * n, y.den * d, true)
  else
    return FakeFmpqMat(y.num * divexact(n, g), d * divexact(y.den, g), true)
  end
end

*(x::FakeFmpqMat, y::fmpq) = y * x

*(x::FakeFmpqMat, y::Rational{<:Integer}) = x * fmpq(y)

*(x::Rational{<:Integer}, y::FakeFmpqMat) = y * x

################################################################################
#
#  Comparison
#
################################################################################

==(x::FakeFmpqMat, y::FakeFmpqMat) = (x.num == y.num) && (x.den == y.den)

isequal(x::FakeFmpqMat, y::FakeFmpqMat) = (x.num == y.num) && (x.den == y.den)

################################################################################
#
#  Conversion
#
################################################################################

to_array(x::FakeFmpqMat) = (x.num, x.den)

function to_fmpz_mat(x::FakeFmpqMat)
  !isone(x.den) && error("Denominator has to be 1")
  return numerator(x)
end

function FakeFmpqMat(x::Vector{fmpq})
  dens = fmpz[denominator(x[i]) for i=1:length(x)]
  den = lcm(dens)
  M = zero_matrix(FlintZZ, 1, length(x))
  for i in 1:length(x)
    M[1,i] = numerator(x[i])*divexact(den, dens[i])
  end
  return FakeFmpqMat(M,den)
end

function fmpq_mat(x::FakeFmpqMat)
  z = fmpq(1, x.den) * fmpq_mat(x.num)
  return z
end

function fmpq_mat(x::fmpz_mat)
  z = zero_matrix(FlintQQ, nrows(x), ncols(x))
  ccall((:fmpq_mat_set_fmpz_mat, libflint), Nothing, (Ref{fmpq_mat}, Ref{fmpz_mat}), z, x)
  return z
end

function _fmpq_mat_to_fmpz_mat_den(x::fmpq_mat)
  z = zero_matrix(FlintZZ, nrows(x), ncols(x))
  d = fmpz()
  ccall((:fmpq_mat_get_fmpz_mat_matwise, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz}, Ref{fmpq_mat}), z, d, x)
  return z, d
end

################################################################################
#
#  Hermite normal form for numerator
#
################################################################################

function hnf!(x::FakeFmpqMat, shape = :lowerleft)
  x.num = _hnf(x.num, shape)
  return x
end

function hnf(x::FakeFmpqMat, shape = :lowerleft; triangular_top::Bool = false, compute_det::Bool = false)
  if triangular_top
    @assert ncols(x) <= nrows(x)
    z = one(FlintZZ)
    for i in 1:(ncols(x) - 1)
      for j in (i + 1):ncols(x)
        @assert iszero(x.num[i, j])
      end
    end
    for i in 1:ncols(x)
      z = mul!(z, z, x.num[i, i])
    end
    h = _hnf_modular_eldiv(x.num, z, shape)
  elseif compute_det
    d = det(numerator(x))
    h = _hnf_modular_eldiv(x.num, d, shape)
  else
    h = _hnf(x.num, shape)
  end
  return FakeFmpqMat(h, denominator(x))
end

function hnf_modular_eldiv(x::FakeFmpqMat, g::fmpz, shape = :lowerleft)
  h = _hnf_modular_eldiv(x.num, g, shape)
  return FakeFmpqMat(h, denominator(x))
end

function hnf!!(x::FakeFmpqMat, shape = :lowerleft)
  _hnf!(x.num, shape)
end

################################################################################
#
#  Sub
#
################################################################################

function sub(x::FakeFmpqMat, r::UnitRange{Int}, c::UnitRange{Int})
  z = FakeFmpqMat(sub(x.num, r, c), x.den)
  return z
end

function Base.deepcopy_internal(x::FakeFmpqMat, dict::IdDict)
  z = FakeFmpqMat()
  z.num = Base.deepcopy_internal(x.num, dict)
  z.den = Base.deepcopy_internal(x.den, dict)
  z.rows = nrows(x)
  z.cols = ncols(x)
  return z
end

################################################################################
#
#  Zero row
#
################################################################################

function iszero_row(M::FakeFmpqMat, i::Int)
  return iszero_row(M.num, i)
end

################################################################################
#
#  Zero row
#
################################################################################

function iszero_column(M::FakeFmpqMat, i::Int)
  return iszero_column(M.num, i)
end

################################################################################
#
#  Determinant
#
################################################################################

function det(x::FakeFmpqMat)
  nrows(x) != ncols(x) && error("Matrix must be square")

  return det(x.num)//(x.den)^nrows(x)
end

################################################################################
#
#  Rank
#
################################################################################

rank(x::FakeFmpqMat) = rank(x.num)

################################################################################
#
#  vcat / hcat
#
################################################################################

function vcat(M::FakeFmpqMat, N::FakeFmpqMat)
  g = gcd(denominator(M, copy = false), denominator(N, copy = false))
  d1 = divexact(denominator(M, copy = false), g)
  d2 = divexact(denominator(N, copy = false), g)
  M1 = numerator(M, copy = false)*d2
  N1 = numerator(N, copy = false)*d1
  return FakeFmpqMat(vcat(M1, N1), g*d1*d2)
end

function reduce(::typeof(vcat), A::Vector{FakeFmpqMat})
  if length(A) == 1
    return A[1]
  end

  d = reduce(lcm, (denominator(a) for a in A), init = fmpz(1))
  res = zero_matrix(FlintZZ, sum(nrows, A), ncols(A[1]))
  k = 1
  for i in 1:length(A)
    _copy_matrix_into_matrix(res, k, 1, divexact(d, denominator(A[i])) * numerator(A[i], copy = false))
    k += nrows(A[i])
  end
  return FakeFmpqMat(res, d)
end

function hcat(M::FakeFmpqMat, N::FakeFmpqMat)
  g = gcd(denominator(M, copy = false), denominator(N, copy = false))
  d1 = divexact(denominator(M, copy = false), g)
  d2 = divexact(denominator(N, copy = false), g)
  M1 = numerator(M, copy = false)*d2
  N1 = numerator(N, copy = false)*d1
  return FakeFmpqMat(hcat(M1, N1), g*d1*d2)
end

################################################################################
#
#  Transpose
#
################################################################################

function transpose(M::FakeFmpqMat)
  return FakeFmpqMat(transpose(numerator(M, copy = false)), denominator(M))
end


text:
####################################################################
# In this example we demonstrate using the Proximal Gradient solver
# with a prox operator that enforces that the real part of the
# eigenvalues alpha is nonpositive.
#
# We test the l2 penalty here (a more or less arbitrary choice)

using RobustDMD

T = Float32

# using PyPlot
#--------------------------------------------------------------------
# Generate DMD Synthetic Data
#--------------------------------------------------------------------
println("generate synthetic data...")
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
h = n; # h is not used by this solver (only Trim)
k = 3;      # number of modes
# generate data
sigma = T(1e-4); # size of background noise
mu = T(0.0); # size of spikes (set to zero here for l2 fitting example)
p = T(0.1); # frequency of spikes (percentage of corrupted entries)
X, Xclean, t, alphat, betat = genDMD(m, n, k, sigma, mu; seed=123456, p=p);

#--------------------------------------------------------------------
# Initial Guess
#--------------------------------------------------------------------
println("initialize alpha and B...")
# create initial guess using trap. rule with exact dmd
alpha0, B0 = dmdexactestimate(m,n,k,X,t);

#--------------------------------------------------------------------
# Proximal Gradient L2
#--------------------------------------------------------------------
println("Proximal Gradient l2 loss experiment...")
# loss functions
lossf = (z) -> l2_func(z);
lossg = (z) -> l2_grad!(z);
l2Params = DMDParams(k, X, t, lossf, lossg);
l2DMD = DMDVariables(alpha0, B0, l2Params);
# wrapper: sets up inner solve specialized for l2
l2Svars = DMDVPSolverVariablesLSQ(l2Params)

# prox function, projects to real part <= 0.0
function prox_mr(αr)
    k = length(αr) >> 1;
    for i = 1:k
        αr[i<<1-1] = min(T(0.0), αr[i<<1-1]);
    end
end
options = PG_options(1000,T(1e-6),true,true,10,prox_mr);

# apply solver
stats = PG_solve_DMD!(l2DMD, l2Params, l2Svars, options);

# check performance
err0 = besterrperm(alpha0, alphat);
err1 = besterrperm(l2DMD.alpha, alphat);

println("err exact DMD ..... ", err0)
println("err l2 fit ........ ", err1)