text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# LQ Factorizations
"""
    LQ <: Factorization

Matrix factorization type of the `LQ` factorization of a matrix `A`. The `LQ`
decomposition is the [`QR`](@ref) decomposition of `transpose(A)`. This is the return
type of [`lq`](@ref), the corresponding matrix factorization function.

If `S::LQ` is the factorization object, the lower triangular component can be
obtained via `S.L`, and the orthogonal/unitary component via `S.Q`, such that
`A ≈ S.L*S.Q`.

Iterating the decomposition produces the components `S.L` and `S.Q`.

# Examples
```jldoctest
julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> S = lq(A)
LQ{Float64,Array{Float64,2}} with factors L and Q:
[-8.60233 0.0; 4.41741 -0.697486]
[-0.581238 -0.813733; -0.813733 0.581238]

julia> S.L * S.Q
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> l, q = S; # destructuring via iteration

julia> l == S.L &&  q == S.Q
true
```
"""
struct LQ{T,S<:AbstractMatrix{T}} <: Factorization{T}
    factors::S
    τ::Vector{T}

    function LQ{T,S}(factors, τ) where {T,S<:AbstractMatrix{T}}
        require_one_based_indexing(factors)
        new{T,S}(factors, τ)
    end
end
LQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = LQ{T,typeof(factors)}(factors, τ)
function LQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T}
    LQ(convert(AbstractMatrix{T}, factors), convert(Vector{T}, τ))
end

# iteration for destructuring into components
Base.iterate(S::LQ) = (S.L, Val(:Q))
Base.iterate(S::LQ, ::Val{:Q}) = (S.Q, Val(:done))
Base.iterate(S::LQ, ::Val{:done}) = nothing

struct LQPackedQ{T,S<:AbstractMatrix} <: AbstractMatrix{T}
    factors::Matrix{T}
    τ::Vector{T}
    LQPackedQ{T,S}(factors::AbstractMatrix{T}, τ::Vector{T}) where {T,S<:AbstractMatrix} = new(factors, τ)
end
LQPackedQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = LQPackedQ{T,typeof(factors)}(factors, τ)

"""
    lq!(A) -> LQ

Compute the [`LQ`](@ref) factorization of `A`, using the input
matrix as a workspace. See also [`lq`](@ref).
"""
lq!(A::StridedMatrix{<:BlasFloat}) = LQ(LAPACK.gelqf!(A)...)
"""
    lq(A) -> S::LQ

Compute the LQ decomposition of `A`. The decomposition's lower triangular
component can be obtained from the [`LQ`](@ref) object `S` via `S.L`, and the
orthogonal/unitary component via `S.Q`, such that `A ≈ S.L*S.Q`.

Iterating the decomposition produces the components `S.L` and `S.Q`.

The LQ decomposition is the QR decomposition of `transpose(A)`, and it is useful
in order to compute the minimum-norm solution `lq(A) \\ b` to an underdetermined
system of equations (`A` has more columns than rows, but has full row rank).

# Examples
```jldoctest
julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> S = lq(A)
LQ{Float64,Array{Float64,2}} with factors L and Q:
[-8.60233 0.0; 4.41741 -0.697486]
[-0.581238 -0.813733; -0.813733 0.581238]

julia> S.L * S.Q
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> l, q = S; # destructuring via iteration

julia> l == S.L &&  q == S.Q
true
```
"""
lq(A::StridedMatrix{<:BlasFloat})  = lq!(copy(A))
lq(x::Number) = lq(fill(x,1,1))

copy(A::LQ) = LQ(copy(A.factors), copy(A.τ))

LQ{T}(A::LQ) where {T} = LQ(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ))
Factorization{T}(A::LQ{T}) where {T} = A
Factorization{T}(A::LQ) where {T} = LQ{T}(A)
AbstractMatrix(A::LQ) = A.L*A.Q
AbstractArray(A::LQ) = AbstractMatrix(A)
Matrix(A::LQ) = Array(AbstractArray(A))
Array(A::LQ) = Matrix(A)

adjoint(A::LQ) = Adjoint(A)
Base.copy(F::Adjoint{T,<:LQ{T}}) where {T} =
    QR{T,typeof(F.parent.factors)}(copy(adjoint(F.parent.factors)), copy(F.parent.τ))

function getproperty(F::LQ, d::Symbol)
    m, n = size(F)
    if d === :L
        return tril!(getfield(F, :factors)[1:m, 1:min(m,n)])
    elseif d === :Q
        return LQPackedQ(getfield(F, :factors), getfield(F, :τ))
    else
        return getfield(F, d)
    end
end

Base.propertynames(F::LQ, private::Bool=false) =
    (:L, :Q, (private ? fieldnames(typeof(F)) : ())...)

getindex(A::LQPackedQ, i::Integer, j::Integer) =
    lmul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i]

function show(io::IO, ::MIME"text/plain", C::LQ)
    println(io, typeof(C), " with factors L and Q:")
    io = IOContext(io, :compact => true)
    show(io, C.L)
    println(io)
    show(io, C.Q)
end

LQPackedQ{T}(Q::LQPackedQ) where {T} = LQPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ))
AbstractMatrix{T}(Q::LQPackedQ) where {T} = LQPackedQ{T}(Q)
Matrix{T}(A::LQPackedQ) where {T} = convert(Matrix{T}, LAPACK.orglq!(copy(A.factors),A.τ))
Matrix(A::LQPackedQ{T}) where {T} = Matrix{T}(A)
Array{T}(A::LQPackedQ{T}) where {T} = Matrix{T}(A)
Array(A::LQPackedQ) = Matrix(A)

size(F::LQ, dim::Integer) = size(getfield(F, :factors), dim)
size(F::LQ)               = size(getfield(F, :factors))

# size(Q::LQPackedQ) yields the shape of Q's square form
function size(Q::LQPackedQ)
    n = size(Q.factors, 2)
    return n, n
end
function size(Q::LQPackedQ, dim::Integer)
    if dim < 1
        throw(BoundsError())
    elseif dim <= 2 # && 1 <= dim
        return size(Q.factors, 2)
    else # 2 < dim
        return 1
    end
end

## Multiplication by LQ
function lmul!(A::LQ, B::StridedVecOrMat)
    lmul!(LowerTriangular(A.L), view(lmul!(A.Q, B), 1:size(A,1), axes(B,2)))
    return B
end
function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
    TAB = promote_type(TA, TB)
    _cut_B(lmul!(Factorization{TAB}(A), copy_oftype(B, TAB)), 1:size(A,1))
end

## Multiplication by Q
### QB
lmul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
function (*)(A::LQPackedQ, B::StridedVecOrMat)
    TAB = promote_type(eltype(A), eltype(B))
    lmul!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB))
end

### QcB
lmul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} =
    (A = adjA.parent; LAPACK.ormlq!('L','T',A.factors,A.τ,B))
lmul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
    (A = adjA.parent; LAPACK.ormlq!('L','C',A.factors,A.τ,B))

function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat)
    A = adjA.parent
    TAB = promote_type(eltype(A), eltype(B))
    if size(B,1) == size(A.factors,2)
        lmul!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB))
    elseif size(B,1) == size(A.factors,1)
        lmul!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))])
    else
        throw(DimensionMismatch("first dimension of B, $(size(B,1)), must equal one of the dimensions of A, $(size(A))"))
    end
end

### QBc/QcBc
function *(A::LQPackedQ, adjB::Adjoint{<:Any,<:StridedVecOrMat})
    B = adjB.parent
    TAB = promote_type(eltype(A), eltype(B))
    BB = similar(B, TAB, (size(B, 2), size(B, 1)))
    adjoint!(BB, B)
    return lmul!(A, BB)
end
function *(adjA::Adjoint{<:Any,<:LQPackedQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat})
    A, B = adjA.parent, adjB.parent
    TAB = promote_type(eltype(A), eltype(B))
    BB = similar(B, TAB, (size(B, 2), size(B, 1)))
    adjoint!(BB, B)
    return lmul!(adjoint(A), BB)
end

# in-place right-application of LQPackedQs
# these methods require that the applied-to matrix's (A's) number of columns
# match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place
# operation, and the underlying LAPACK routine (ormlq) treats the implicit Q
# as its (nQ-by-nQ) square form)
rmul!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} =
    LAPACK.ormlq!('R', 'N', B.factors, B.τ, A)
rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} =
    (B = adjB.parent; LAPACK.ormlq!('R', 'T', B.factors, B.τ, A))
rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} =
    (B = adjB.parent; LAPACK.ormlq!('R', 'C', B.factors, B.τ, A))

# out-of-place right application of LQPackedQs
#
# LQPackedQ's out-of-place multiplication behavior is context dependent. specifically,
# if the inner dimension in the multiplication is the LQPackedQ's second dimension,
# the LQPackedQ behaves like its square form. if the inner dimension in the
# multiplication is the LQPackedQ's first dimension, the LQPackedQ behaves like either
# its square form or its truncated form depending on the shape of the other object
# involved in the multiplication. we treat these cases separately.
#
# (1) the inner dimension in the multiplication is the LQPackedQ's second dimension.
# in this case, the LQPackedQ behaves like its square form.
#
function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ})
    Q = adjQ.parent
    TR = promote_type(eltype(A), eltype(Q))
    return rmul!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
end
function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ})
    A, Q = adjA.parent, adjQ.parent
    TR = promote_type(eltype(A), eltype(Q))
    C = adjoint!(similar(A, TR, reverse(size(A))), A)
    return rmul!(C, adjoint(AbstractMatrix{TR}(Q)))
end
#
# (2) the inner dimension in the multiplication is the LQPackedQ's first dimension.
# in this case, the LQPackedQ behaves like either its square form or its
# truncated form depending on the shape of the other object in the multiplication.
#
# these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension
# (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q),
# and if so effectively apply Q's square form to A without additional shenanigans; and
# (2) if the preceding dimensions do not match, check whether the appropriate dimension of
# A instead matches the number of rows of the matrix of which Q is a factor (i.e.
# size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending
# A as necessary for check (1) to pass (if possible) and then applying Q's square form
#
function *(A::StridedVecOrMat, Q::LQPackedQ)
    TR = promote_type(eltype(A), eltype(Q))
    if size(A, 2) == size(Q.factors, 2)
        C = copy_oftype(A, TR)
    elseif size(A, 2) == size(Q.factors, 1)
        C = zeros(TR, size(A, 1), size(Q.factors, 2))
        copyto!(C, 1, A, 1, length(A))
    else
        _rightappdimmismatch("columns")
    end
    return rmul!(C, AbstractMatrix{TR}(Q))
end
function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ)
    A = adjA.parent
    TR = promote_type(eltype(A), eltype(Q))
    if size(A, 1) == size(Q.factors, 2)
        C = adjoint!(similar(A, TR, reverse(size(A))), A)
    elseif size(A, 1) == size(Q.factors, 1)
        C = zeros(TR, size(A, 2), size(Q.factors, 2))
        adjoint!(view(C, :, 1:size(A, 1)), A)
    else
        _rightappdimmismatch("rows")
    end
    return rmul!(C, AbstractMatrix{TR}(Q))
end
_rightappdimmismatch(rowsorcols) =
    throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ",
        "must match either (1) the number of columns of the (LQPackedQ) matrix on the right ",
        "or (2) the number of rows of that (LQPackedQ) matrix's internal representation ",
        "(the factorization's originating matrix's number of rows)")))

function (\)(A::LQ{TA},B::StridedVecOrMat{TB}) where {TA,TB}
    S = promote_type(TA,TB)
    m, n = size(A)
    m ≤ n || throw(DimensionMismatch("LQ solver does not support overdetermined systems (more rows than columns)"))
    m == size(B,1) || throw(DimensionMismatch("Both inputs should have the same number of rows"))
    AA = Factorization{S}(A)
    X = _zeros(S, B, n)
    X[1:size(B, 1), :] = B
    return ldiv!(AA, X)
end
# With a real lhs and complex rhs with the same precision, we can reinterpret
# the complex rhs as a real rhs with twice the number of columns
function (\)(F::LQ{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
    require_one_based_indexing(B)
    X = zeros(T, size(F,2), 2*size(B,2))
    X[1:size(B,1), 1:size(B,2)] .= real.(B)
    X[1:size(B,1), size(B,2)+1:size(X,2)] .= imag.(B)
    ldiv!(F, X)
    return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(X, div(length(X), 2), 2))))),
                           isa(B, AbstractVector) ? (size(F,2),) : (size(F,2), size(B,2)))
end

function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T
    require_one_based_indexing(B)
    ldiv!(LowerTriangular(A.L), view(B, 1:size(A,1), axes(B,2)))
    return lmul!(adjoint(A.Q), B)
end

# In LQ factorization, `Q` is expressed as the product of the adjoint of the
# reflectors.  Thus, `det` has to be conjugated.
det(Q::LQPackedQ) = conj(_det_tau(Q.τ))

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using MPI
using CLIMA
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.Vtk

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray, )
else
  const ArrayTypes = (Array, )
end

const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E)
const statenames = ("ρ", "U", "V", "W", "E")

const _nviscstates = 6
const _τ11, _τ22, _τ33, _τ12, _τ13, _τ23 = 1:_nviscstates

const _ngradstates = 3
const _states_for_gradient_transform = (_ρ, _U, _V, _W)

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
  using Random
end

include("mms_solution_generated.jl")

# preflux computation
@inline function preflux(Q, _...)
  γ::eltype(Q) = γ_exact
  @inbounds ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]
  ρinv = 1 / ρ
  u, v, w = ρinv * U, ρinv * V, ρinv * W
  ((γ-1)*(E - ρinv * (U^2 + V^2 + W^2) / 2), u, v, w, ρinv)
end

# max eigenvalue
@inline function wavespeed(n, Q, aux, t, P, u, v, w, ρinv)
  γ::eltype(Q) = γ_exact
  @inbounds abs(n[1] * u + n[2] * v + n[3] * w) + sqrt(ρinv * γ * P)
end

# flux function
cns_flux!(F, Q, VF, aux, t) = cns_flux!(F, Q, VF, aux, t, preflux(Q)...)

@inline function cns_flux!(F, Q, VF, aux, t, P, u, v, w, ρinv)
  @inbounds begin
    ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]

    τ11, τ22, τ33 = VF[_τ11], VF[_τ22], VF[_τ33]
    τ12 = τ21 = VF[_τ12]
    τ13 = τ31 = VF[_τ13]
    τ23 = τ32 = VF[_τ23]

    # inviscid terms
    F[1, _ρ], F[2, _ρ], F[3, _ρ] = U          , V          , W
    F[1, _U], F[2, _U], F[3, _U] = u * U  + P , v * U      , w * U
    F[1, _V], F[2, _V], F[3, _V] = u * V      , v * V + P  , w * V
    F[1, _W], F[2, _W], F[3, _W] = u * W      , v * W      , w * W + P
    F[1, _E], F[2, _E], F[3, _E] = u * (E + P), v * (E + P), w * (E + P)

    # viscous terms
    F[1, _U] -= τ11; F[2, _U] -= τ12; F[3, _U] -= τ13
    F[1, _V] -= τ21; F[2, _V] -= τ22; F[3, _V] -= τ23
    F[1, _W] -= τ31; F[2, _W] -= τ32; F[3, _W] -= τ33

    F[1, _E] -= u * τ11 + v * τ12 + w * τ13
    F[2, _E] -= u * τ21 + v * τ22 + w * τ23
    F[3, _E] -= u * τ31 + v * τ32 + w * τ33
  end
end

# Compute the velocity from the state
@inline function velocities!(vel, Q, _...)
  @inbounds begin
    # ordering should match states_for_gradient_transform
    ρ, U, V, W = Q[1], Q[2], Q[3], Q[4]
    ρinv = 1 / ρ
    vel[1], vel[2], vel[3] = ρinv * U, ρinv * V, ρinv * W
  end
end

# Visous flux
@inline function compute_stresses!(VF, grad_vel, _...)
  μ::eltype(VF) = μ_exact
  @inbounds begin
    dudx, dudy, dudz = grad_vel[1, 1], grad_vel[2, 1], grad_vel[3, 1]
    dvdx, dvdy, dvdz = grad_vel[1, 2], grad_vel[2, 2], grad_vel[3, 2]
    dwdx, dwdy, dwdz = grad_vel[1, 3], grad_vel[2, 3], grad_vel[3, 3]

    # strains
    ϵ11 = dudx
    ϵ22 = dvdy
    ϵ33 = dwdz
    ϵ12 = (dudy + dvdx) / 2
    ϵ13 = (dudz + dwdx) / 2
    ϵ23 = (dvdz + dwdy) / 2

    # deviatoric stresses
    VF[_τ11] = 2μ * (ϵ11 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ22] = 2μ * (ϵ22 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ33] = 2μ * (ϵ33 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ12] = 2μ * ϵ12
    VF[_τ13] = 2μ * ϵ13
    VF[_τ23] = 2μ * ϵ23
  end
end

@inline function stresses_penalty!(VF, nM, velM, QM, aM, velP, QP, aP, t)
  @inbounds begin
    n_Δvel = similar(VF, Size(3, 3))
    for j = 1:3, i = 1:3
      n_Δvel[i, j] = nM[i] * (velP[j] - velM[j]) / 2
    end
    compute_stresses!(VF, n_Δvel)
  end
end

@inline stresses_boundary_penalty!(VF, _...) = VF.=0

# initial condition
function initialcondition!(dim, Q, t, x, y, z, _...)
  DFloat = eltype(Q)
  ρ::DFloat = ρ_g(t, x, y, z, dim)
  U::DFloat = U_g(t, x, y, z, dim)
  V::DFloat = V_g(t, x, y, z, dim)
  W::DFloat = W_g(t, x, y, z, dim)
  E::DFloat = E_g(t, x, y, z, dim)

  if integration_testing
    @inbounds Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E] = ρ, U, V, W, E
  else
    @inbounds Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E] =
    10+rand(), rand(), rand(), rand(), 10+rand()
  end
end

const _nauxstate = 3
const _a_x, _a_y, _a_z = 1:_nauxstate
@inline function auxiliary_state_initialization!(aux, x, y, z)
  @inbounds begin
    aux[_a_x] = x
    aux[_a_y] = y
    aux[_a_z] = z
  end
end

@inline function source3D!(S, Q, aux, t)
  @inbounds begin
    x,y,z = aux[_a_x], aux[_a_y], aux[_a_z]
    S[_ρ] = Sρ_g(t, x, y, z, Val(3))
    S[_U] = SU_g(t, x, y, z, Val(3))
    S[_V] = SV_g(t, x, y, z, Val(3))
    S[_W] = SW_g(t, x, y, z, Val(3))
    S[_E] = SE_g(t, x, y, z, Val(3))
  end
end

@inline function source2D!(S, Q, aux, t)
  @inbounds begin
    x,y,z = aux[_a_x], aux[_a_y], aux[_a_z]
    S[_ρ] = Sρ_g(t, x, y, z, Val(2))
    S[_U] = SU_g(t, x, y, z, Val(2))
    S[_V] = SV_g(t, x, y, z, Val(2))
    S[_W] = SW_g(t, x, y, z, Val(2))
    S[_E] = SE_g(t, x, y, z, Val(2))
  end
end

@inline function bcstate2D!(QP, QVP, auxP, nM, QM, QVM, auxM, bctype, t, _...)
  @inbounds begin
    x, y, z = auxM[_a_x], auxM[_a_y], auxM[_a_z]
    if integration_testing
      initialcondition!(Val(2), QP, t, x, y, z)
    else
      for s = 1:length(QP)
        QP[s] = QM[length(QP)+1-s]
      end
      for s = 1:_nviscstates
        QVP[s] = QVM[s]
      end
    end
  end
  nothing
end

@inline function bcstate3D!(QP, QVP, auxP, nM, QM, QVM, auxM, bctype, t, _...)
  @inbounds begin
    x, y, z = auxM[_a_x], auxM[_a_y], auxM[_a_z]
    if integration_testing
      initialcondition!(Val(3), QP, t, x, y, z)
    else
      for s = 1:length(QP)
        QP[s] = QM[length(QP)+1-s]
      end
      for s = 1:_nviscstates
        QVP[s] = QVM[s]
      end
    end
  end
  nothing
end

function run(mpicomm, ArrayType, dim, topl, warpfun, N, timeend, DFloat, dt)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DFloat,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = warpfun,
                                         )

  # spacedisc = data needed for evaluating the right-hand side function
  numflux!(x...) = NumericalFluxes.rusanov!(x..., cns_flux!, wavespeed,
                                            preflux)
  bcstate! = dim == 2 ?  bcstate2D! : bcstate3D!
  numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(x..., cns_flux!,
                                                            bcstate!,
                                                            wavespeed, preflux)
  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = _nstate,
                           flux! = cns_flux!,
                           numerical_flux! = numflux!,
                           numerical_boundary_flux! = numbcflux!,
                           number_gradient_states = _ngradstates,
                           states_for_gradient_transform =
                             _states_for_gradient_transform,
                           number_viscous_states = _nviscstates,
                           gradient_transform! = velocities!,
                           viscous_transform! = compute_stresses!,
                           viscous_penalty! = stresses_penalty!,
                           viscous_boundary_penalty! =
                           stresses_boundary_penalty!,
                           auxiliary_state_length = _nauxstate,
                           auxiliary_state_initialization! =
                           auxiliary_state_initialization!,
                           source! = dim == 2 ? source2D! : source3D!)

  # This is a actual state/function that lives on the grid
  initialcondition(Q, x...) = initialcondition!(Val(dim), Q, DFloat(0), x...)
  Q = MPIStateArray(spacedisc, initialcondition)

  lsrk = LowStorageRungeKutta(spacedisc, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  npoststates = 5
  _P, _u, _v, _w, _ρinv = 1:npoststates
  postnames = ("P", "u", "v", "w", "ρinv")
  postprocessarray = MPIStateArray(spacedisc; nstate=npoststates)

  step = [0]
  mkpath("vtk")
  cbvtk = GenericCallbacks.EveryXSimulationSteps(100) do (init=false)
    DGBalanceLawDiscretizations.dof_iteration!(postprocessarray, spacedisc,
                                               Q) do R, Q, QV, aux
      @inbounds let
        (R[_P], R[_u], R[_v], R[_w], R[_ρinv]) = preflux(Q)
      end
    end

    outprefix = @sprintf("vtk/cns_%dD_mpirank%04d_step%04d", dim,
                         MPI.Comm_rank(mpicomm), step[1])
    @debug "doing VTK output" outprefix
    writevtk(outprefix, Q, spacedisc, statenames,
             postprocessarray, postnames)
    step[1] += 1
    nothing
  end

  # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, ))
  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  # Print some end of the simulation information
  engf = norm(Q)
  if integration_testing
    Qe = MPIStateArray(spacedisc,
                       (Q, x...) -> initialcondition!(Val(dim), Q,
                                                      DFloat(timeend), x...))
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf/eng0 engf-eng0 errf errf / engfe
  else
    @info @sprintf """Finished
    norm(Q)            = %.16e
    norm(Q) / norm(Q₀) = %.16e
    norm(Q) - norm(Q₀) = %.16e""" engf engf/eng0 engf-eng0
  end
  integration_testing ? errf : (engf / eng0)
end

using Test
let
  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
    device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end

  polynomialorder = 4
  base_num_elem = 4
  if integration_testing
    expected_result = Array{Float64}(undef, 2, 3) # dim-1, lvl
    expected_result[1,1] = 1.6687745307357629e-01
    expected_result[1,2] = 5.4179126727473799e-03
    expected_result[1,3] = 2.3066157635992409e-04
    expected_result[2,1] = 3.3669188610024728e-02
    expected_result[2,2] = 1.7603468555920912e-03
    expected_result[2,3] = 9.1108572847298699e-05
    lvls = size(expected_result, 2)
  else
    expected_result = Dict{Tuple{Int64, Int64, DataType}, AbstractFloat}()
    expected_result[2, 1, Float64] = 9.9075897196717488e-01
    expected_result[3, 1, Float64] = 1.0099522817205739e+00
    expected_result[2, 3, Float64] = 9.9072475063319887e-01
    expected_result[3, 3, Float64] = 1.0099521111150005e+00
    lvls = 1
  end

  for ArrayType in ArrayTypes
    for DFloat in (Float64,) #Float32)
      result = zeros(DFloat, lvls)
      for dim = 2:3
        for l = 1:lvls
          integration_testing || Random.seed!(0)
          if dim == 2
            Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
            brickrange = (range(DFloat(0); length=Ne[1]+1, stop=1),
                          range(DFloat(0); length=Ne[2]+1, stop=1))
            topl = BrickTopology(mpicomm, brickrange,
                                 periodicity = (false, false))
            dt = 1e-2 / Ne[1]
            warpfun = (x, y, _) -> begin
              (x + sin(x*y), y + sin(2*x*y), 0)
            end

          elseif dim == 3
            Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
            brickrange = (range(DFloat(0); length=Ne[1]+1, stop=1),
                          range(DFloat(0); length=Ne[2]+1, stop=1),
            range(DFloat(0); length=Ne[2]+1, stop=1))
            topl = BrickTopology(mpicomm, brickrange,
                                 periodicity = (false, false, false))
            dt = 5e-3 / Ne[1]
            warpfun = (x, y, z) -> begin
              (x + (x-1/2)*cos(2*π*y*z)/4,
               y + exp(sin(2π*(x*y+z)))/20,
              z + x/4 + y^2/2 + sin(x*y*z))
            end
          end
          timeend = integration_testing ? 1 : 2dt
          nsteps = ceil(Int64, timeend / dt)
          dt = timeend / nsteps

          @info (ArrayType, DFloat, dim)
          result[l] = run(mpicomm, ArrayType, dim, topl, warpfun,
                          polynomialorder, timeend, DFloat, dt)
          if integration_testing
            @test result[l] ≈ DFloat(expected_result[dim-1, l])
          else
            @test result[l] ≈ expected_result[dim, MPI.Comm_size(mpicomm), DFloat]
          end
        end
        if integration_testing
          @info begin
            msg = ""
            for l = 1:lvls-1
              rate = log2(result[l]) - log2(result[l+1])
              msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
            end
            msg
          end
        end
      end
    end
  end
end

isinteractive() || MPI.Finalize()

nothing

ext: jl

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
using GeoPhyInv
using Test

# Can perform inversion of synthetic scenarios.
# First, the seismic data are modeled as in the forward problem. Then the
# data are used to perform full waveform inversion (FWI). The inverse
# problem estimates
# the Earth models and the source and receiver filters
# that resulted from the data.
# This task is necessary to test the performance of the inversion algorithm
# in various geological scenarios using different acquisition parameters.

medium = Medium(:acou_homo2);
update!(medium, [:vp,:rho], randn_perc=1)

medium0 = Medium(:acou_homo2);
update!(medium0, [:vp,:rho], randn_perc=1)

ageom=AGeom(medium.mgrid, :surf)
wav, tgrid=ricker(medium, 3, 1.0)
srcwav=SrcWav(tgrid,ageom,[:p])
update!(srcwav,[:p], wav)

parameterization=[:χvp, :χrho, :null]

mgrid=medium.mgrid

@testset "test parallel implementation during gradient" begin
	for attrib_mod in [FdtdAcoustic(), FdtdAcousticBorn()]
		global pa=SeisInvExpt(attrib_mod, Migr(), srcwav=srcwav, ageom=ageom, tgrid=tgrid, mediumm=medium0,
				     mediumm_obs=medium,
				     mediumm0=medium0,
				     igrid_interp_scheme=:B2,
				     igrid=broadcast(x->range(x[1],stop=x[end],step=300.),mgrid),
				     parameterization=parameterization,   verbose=false,
				     nworker=1)

		global pa_parallel=SeisInvExpt(attrib_mod, Migr(), srcwav=srcwav, ageom=ageom, tgrid=tgrid, mediumm=medium0,
				     mediumm_obs=medium,
				     mediumm0=medium0,
				     igrid_interp_scheme=:B2,
				     igrid=broadcast(x->range(x[1],stop=x[end],step=300.),mgrid),
				     parameterization=parameterization,   verbose=false,
				     nworker=nothing)

		result=update!(pa)

		result_parallel=update!(pa_parallel)

		@test result[2] ≈ result_parallel[2]
	end
end

@testset "Testing Born Modeling and its gradient" begin
	global expt=x->SeisInvExpt(FdtdAcousticBorn(), x,
			     srcwav=srcwav, ageom=ageom, tgrid=tgrid, mediumm=medium0,
	     		     mediumm0=medium0,
			     mediumm_obs=medium,
			     igrid_interp_scheme=:B2,
			     igrid=broadcast(x->range(x[1],stop=x[end],step=350.),mgrid),
			     parameterization=parameterization,   verbose=false)

	global pa=expt(LS())

	update!(pa, bounded_flag=true, solver=:ipopt,
			ipopt_options=[["max_iter", 0],["derivative_test", "first-order"]])
	pa=expt(Migr())

	result=update!(pa)

	pa_fd=expt(Migr_FD());
	result_fd=update!(pa_fd)

	f=squared_euclidean!(nothing, result[2], result_fd[2], nothing, norm_flag=true)

	@test f<1e-15
end

@testset "Testing gradient LS FWI" begin
	expt=x->SeisInvExpt(FdtdAcoustic(), x,
			     srcwav=srcwav, ageom=ageom, tgrid=tgrid, mediumm=medium0,
	     		     mediumm0=medium0,
			     mediumm_obs=medium,
			     igrid_interp_scheme=:B2,
			     igrid=broadcast(x->range(x[1],stop=x[end],step=350.),mgrid),
			     parameterization=parameterization,   verbose=false)
	pa=expt(LS())

	update!(pa, bounded_flag=true, solver=:ipopt,
			ipopt_options=[["max_iter", 0],["derivative_test", "first-order"]])
	pa=expt(Migr())

	result=update!(pa)

	pa_fd=expt(Migr_FD());
	result_fd=update!(pa_fd)

	f=squared_euclidean!(nothing, result[2], result_fd[2], nothing, norm_flag=true)

	@test f<1e-15
end


text:
"""
    struct VariableSparsity <: Sparsity end

Variable or correlative sparsity as developed in [WSMM06].

[WSMM06] Waki, Hayato, Sunyoung Kim, Masakazu Kojima, and Masakazu Muramatsu. "Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity." SIAM Journal on Optimization 17, no. 1 (2006): 218-242.
"""
struct VariableSparsity <: Sparsity end

const CEG = ChordalExtensionGraph

function csp_graph(poly::MP.APL, ::FullSpace)
    G = CEG.LabelledGraph{MP.variable_union_type(poly)}()
    for mono in MP.monomials(poly)
        CEG.add_clique!(G, MP.effective_variables(mono))
    end
    return G
end

function csp_graph(poly::MP.APL, domain::AbstractAlgebraicSet)
    G = csp_graph(poly, FullSpace())
    for p in equalities(domain)
        CEG.add_clique!(G, MP.effective_variables(p))
    end
    return G
end

function csp_graph(poly::MP.APL, domain::BasicSemialgebraicSet)
    G = csp_graph(poly, domain.V)
    for p in inequalities(domain)
        CEG.add_clique!(G, MP.effective_variables(p))
    end
    return G
end

function chordal_csp_graph(poly::MP.APL, domain::AbstractBasicSemialgebraicSet)
    H, cliques = CEG.chordal_extension(csp_graph(poly, domain), CEG.GreedyFillIn())
    for clique in cliques
        sort!(clique, rev=true)
        unique!(clique)
    end
    return H, cliques

end

function sparsity(poly::MP.AbstractPolynomial, domain::BasicSemialgebraicSet,
                  sp::VariableSparsity, certificate::Putinar)
    H, cliques = chordal_csp_graph(poly, domain)
    function bases(q)
        return [
            maxdegree_gram_basis(certificate.basis, clique,
                                 multiplier_maxdegree(certificate.maxdegree, q))
            for clique in cliques if variables(q) ⊆ clique
        ]
    end
    return bases(poly), map(bases, domain.p)
end


text:
#=##############################################################################
# DESCRIPTION
    Tools for formatting PSU-WOPWOP's inputs/outputs.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez and Tyler Critchfield
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Nov 2019
  * License   : MIT
=###############################################################################

"""
    `read_wopwopoutput(outputname; read_path="", tec=false)`

Reads an output PSU-WOPWOP field of observer data. It returns `(header, field)`
where `header[k]` is a string with the name of the k-th field, and
`field[i, j, t, k]` is the value at the node of coordinates (i, j) at time t
in the k-th field.

Give it `tec=true` to read a TecPlot file (output of a single observer), and it
will return `(header, field)` where `field[t, k]` is the value at time t in the
k-th field.

# ARGUMENTS
* `outputname::String`           : PSU-WOPWOP output file to read. For instance,
                                    `outputname=pressure` will look for the files
                                    pressure.nam and pressure.fn under path
                                    `read_path`.

# OPTIONAL ARGUMENTS
* `read_path::String`           : Path from where to read the inputs files.
"""
function read_wopwopoutput(outputname::String; read_path="",
                            verbose=true, verbose_level=1,
                            v_lvl=0, tec=false)

    if tec
        if verbose && verbose_level>=0
            println("\t"^v_lvl*"Reading output file $(outputname*".tec") ...")
        end

        f = open(joinpath(read_path, outputname*".tec"), "r")

        readline(f) # Title
        header = split(replace(readline(f)[13:end], "\"" => ""), ", ")
        readline(f) # Info
        readline(f) # Empty line

        # Data
        data = []

        for ln in eachline(f)

            # Get rid of space at start of line
            while ln[1] == ' '
                ln = ln[2:end]
            end

            # Get rid of space at end of line
            while ln[end] == ' '
                ln = ln[1:end-1]
            end

            # Get rid of double spaces
            while occursin("  ", ln)
                ln = replace(ln, "  " => " ")
            end

            # Parse data
            push!(data, Meta.parse.(split(ln, " ")))
        end

        field = hcat(data...)'

        close(f)
    else
        if verbose && verbose_level>=0
            println("\t"^v_lvl*"Reading output file $(outputname*".nam") ...")
        end

        # Read header file
        f = open(joinpath(read_path, outputname*".nam"), "r")
        header = readlines(f)
        close(f)

        # Remove invalid characters from headers
        for str in [" ", "(", ")"]
            for i in 1:length(header)
                header[i] = replace(header[i], str => "")
            end
        end

        # Read function file
        f = open(joinpath(read_path, outputname*".fn"), "r")
        imax, jmax, tmax, fieldmax = [Meta.parse(elem) for elem in split(readline(f))]

        if verbose && verbose_level>=0
            println("\t"^(v_lvl+1)*"Found $(imax)x$(jmax) grid in WOPWOP output with"*
                        " $tmax time steps.")
        end

        # Read fields
        field = zeros(imax, jmax, tmax, fieldmax)
        for fieldi in 1:fieldmax # Iterate over fields

            if verbose && verbose_level>=1
                println("\t"^(v_lvl+1)*"Reading field $(header[fieldi])...")
            end

            for t in 1:tmax # Iterate over time steps

                if verbose && (t-1)%ceil(Int, tmax/4)==0 && verbose_level>=1
                    println("\t"^(v_lvl+2)*"Reading time t=$t out of $tmax")
                end

                for j in 1:jmax
                    for i in 1:imax
                        # NOTE: Shouldn't I be iterating over in the inner loop?
                        field[i, j, t, fieldi] = Meta.parse(readline(f))
                    end
                end
            end
        end

        close(f)
    end

    return header, field
end

"""
    `read_wopwoploading(filename; read_path="", verbose=false, v_lvl=0)`

Read a PSU-WOPWOP functional file containing loading data (v1.0). Returns
`(hdr, names, timeinfo, dims, time, data)` where data[k][:, i, j] is the value at
the i-th node (or face) at the j-th time in the k-th zone
"""
function read_wopwoploading(filename; read_path="",
                                        verbose=false, verbose_lvl=0, v_lvl=0)

    # Open PSU-WOPWOP file
    f = open(joinpath(read_path, filename), "r")

    # Read format header
    header = []
    push!(header, read(f, Int32, 1))                #1 Magic number
    push!(header, read(f, Int32, 2))                #2 Version
    chars = [Char(c) for c in read(f, Int8, 1024)]  #3 Comments
    str = ""; for c in chars; str *= c; end;
    push!(header, str)
    push!(header, read(f, Int32, 1))                #4 Functional file
    nz = read(f, Int32, 1)[1]                       #5 Number of zones
    push!(header, nz)
    push!(header, read(f, Int32, 1))                #6 1==structured, 2==unstructured
    push!(header, read(f, Int32, 1))                #7 1==Constant, 2==Periodic, 3==Aperiodic
    ncnt = read(f, Int32, 1)[1]                     #8 1==node-centered, 2==face-centered
    push!(header, ncnt)
    push!(header, read(f, Int32, 1))                #9 1==pressure, 2==loading, 3==flow
    push!(header, read(f, Int32, 1))                #10 1==stationary ground fixed, 2==rotating ground fixed, 3==patch fixed
    push!(header, read(f, Int32, 1))                #11 1==single, 2==double
    push!(header, read(f, Int32, 1))                #12 something
    push!(header, read(f, Int32, 1))                #13 something else

    if verbose
        println("\t"^(v_lvl)*"$(header[1])\t # Magic number")
        println("\t"^(v_lvl)*"$(header[2])\t # Version")
        println("\t"^(v_lvl)*"$(header[3])\t # Comment")
        println("\t"^(v_lvl)*"$(header[4])\t # Functional file")
        println("\t"^(v_lvl)*"$(header[5])\t\t # Number of zones")
        println("\t"^(v_lvl)*"$(header[6])\t # 1==structured, 2==unstructured")
        println("\t"^(v_lvl)*"$(header[7])\t # 1==Constant, 2==Periodic, 3==Aperiodic")
        println("\t"^(v_lvl)*"$(header[8])\t\t # 1==node-centered, 2==face-centered")
        println("\t"^(v_lvl)*"$(header[9])\t # 1==pressure, 2==loading, 3==flow")
        println("\t"^(v_lvl)*"$(header[10])\t # 1==stationary ground, 2==rotating ground, 3==patch")
        println("\t"^(v_lvl)*"$(header[11])\t # 1==single, 2==double")
        println("\t"^(v_lvl)*"$(header[12])\t # something")
        println("\t"^(v_lvl)*"$(header[13])\t # something else")
    end

    # Error cases: Current implementation only takes unstructured, constant,
    # face centered patches
    Int(header[2][1]) != 1 ? error("Only v1.0 is supported") :
    Int(header[4][1]) != 2 ? error("File is not flagged as a functional file") :
    # Int(header[6][1]) != 2 ? error("Only unstructured patches are supported") :
    # Int(header[7][1]) != 1 ? error("Only constant patches are supported") :
    # Int(header[8][1]) != 2 ? error("Only face-centered patches are supported") :
    Int(header[9][1]) != 2 ? error("Only loading vectors are supported") :
    Int(header[11][1]) != 1 ? error("Only single-precision floats are supported") :
    nothing;

    # Zone specification
    zone_specs = zeros(Int, read(f, Int32))     # Read number of zones
    for i in 1:length(zone_specs)               # Read each zone number
        zone_specs[i] = Int(read(f, Int32))
    end

    if verbose
        println("\t"^(v_lvl+1)*"Zone specification: $zone_specs")
    end

    names = String[]
    timeinfo = []
    dims = []
    ndata = zeros(Int, length(zone_specs))
    nt = nothing

    # Data header
    for zi in 1:length(zone_specs)

        # Name
        chars = [Char(c) for c in read(f, Int8, 32)]
        str = ""; for c in chars; str *= c; end;
        push!(names, str)

        # Time information
        if Int(header[7][1])==2         # Periodic case
            push!(timeinfo, Any[read(f, Float32), read(f, Int32)])
            ntimes = timeinfo[end][2]

        elseif Int(header[7][1])==3     # Aperiodic case
            push!(timeinfo, read(f, Int32))
            ntimes = timeinfo[end]

        else                            # Constant case
            push!(timeinfo, 1)
            ntimes = 1
        end

        if nt != nothing && nt != ntimes
            error("Logic error: Got different times ($(nt)!=$(ntimes))")
        end
        nt = ntimes

        # Dimensions
        if Int(header[6][1])==1         # Structured case
            push!(dims, read(f, Int32, 2))
            ndata[zi] = dims[end][1]*dims[end][2]
        else
            push!(dims, read(f, Int32))
            ndata[zi] = dims[end]
        end

        if verbose
            println("\t"^(v_lvl+1)*"$(replace(names[end], " " => ""))"*
                    "\t$(timeinfo[end])\t$(dims[end]) # Name, Time info, dims")
        end
    end

    # Time of every data block of every zone
    time = [zeros(nt) for i in 1:length(zone_specs)]

    # Data of every zone, where data[k][:, i, j] is the value at the i-th node
    # (or face) at the j-th time in the k-th zone
    data = [zeros(3, nd, nt) for nd in ndata]

    for ti in 1:nt                          # Iterate over time entries
        for zi in 1:length(zone_specs)      # Iterate over zones

            if Int(header[7][1])!=1
                time[zi][ti] = Float64(read(f, Float32))
            end

            # NOTE: Here I assume is loading data (three-dimensional)
            for dim in 1:3                  # Iterate over spatial dimension
                for di in 1:ndata[zi]       # Iterate over nodes / faces

                    data[zi][dim, di, ti] = Float64(read(f, Float32))

                end
            end
        end
    end

    hdr = Dict(
                "version"   => header[2],
                "comment"   => header[3],
                "zones"     => Int(header[5][1]),
                "structured"=> Int(header[6][1])==1,
                "time"      => Int(header[7][1]),
                "ntimes"    => nt,
                "node-centered"=> Int(header[8][1])==1,
                "data-type" => Int(header[9][1]),
                "frame"     => Int(header[10][1])
              )

    return hdr, names, timeinfo, dims, time, data
end

"""
    `wopwop2vtk(grid, outputname, save_path; read_path="", paraview=false)`

Converts an output field from PSU-WOPWOP generated over a grid of observers to
vtk format. The user must make sure that the `grid` matches the exact grid that
WOPWOP used as an observer to generate the output fields.

# ARGUMENTS
* `grid::gt.AbstractGrid`        : Grid.
* `outputname::String`           : PSU-WOPWOP output file to read. For instance,
                                    `outputname=pressure` will look for the files
                                    pressure.nam and pressure.fn under path
                                    `read_path`.
* `save_path::String`           : Where to save vtk files.

# OPTIONAL ARGUMENTS
* `read_path::String`           : Path from where to read the inputs files.

"""
function wopwop2vtk(grid::gt.AbstractGrid, outputname::String, save_path::String;
                        pref="",
                        read_path="", prompt=true,
                        verbose=true, v_lvl=0, createpath=true,
                        paraview=false)

    if createpath
        gt.create_path(save_path, prompt)
    end

    header, field = read_wopwopoutput(outputname; read_path=read_path,
                                                verbose=verbose, v_lvl=v_lvl)

    tmax = size(field, 3)
    fieldmax = size(field, 4)

    # Add field to grid
    if verbose
        println("\t"^v_lvl*"Generating vtk files...")
    end

    for t in 1:tmax # Iterates over time steps

        if verbose && (t-1)%ceil(Int, tmax/10)==0
            println("\t"^(v_lvl+2)*"Saving vtk time t=$t out of $tmax")
        end

        for fieldi in 1:fieldmax # Iterates over fields

            field_name = header[fieldi]
            field_type = "scalar"
            entry_type = "node"
            field_data = reshape(field[:, :, t, fieldi], grid.nnodes)

            gt.add_field(grid, field_name, field_type, field_data,
                            entry_type; raise_warn=false)
        end

        # Save vtk of this time step
        gt.save(grid, pref*outputname; path=save_path, num=t-1)

    end

    str = pref*outputname*"...vtk;"

    if paraview
        run(`paraview --data=$(save_path)/$(str)`)
    end

    return str
end

"""
    `geomwopwop2vtk(filename; read_path="", save_path=nothing,
pref_save_fname::String="automatic", paraview=false, verbose=false, v_lvl=0)`

Converts an input patch file in PSU-WOPWOP format to VTK format.
"""
function geomwopwop2vtk(filename; read_path="", save_path=nothing,
                        pref_save_fname::String="automatic", paraview=false,
                        verbose=false, v_lvl=0, loading_file=nothing)

    # String carrying vtk files
    if save_path!=nothing && paraview
        vtk_str = save_path*"/"
    end

    # Build output file name prefix
    if pref_save_fname=="automatic"
        doti = findfirst(x->x=='.', filename)
        if length(doti)!=0
            # Get rid of extension
            preff = filename[1:doti[end]-1]
        else
            preff = filename
        end
    else
        preff = pref_save_fname
    end

    # Open PSU-WOPWOP file
    f = open(joinpath(read_path, filename), "r")

    # Read format header
    header = []
    push!(header, read(f, Int32, 1))                #1 Magic number
    push!(header, read(f, Int32, 2))                #2 Version
    chars = [Char(c) for c in read(f, Int8, 32)]    #3 Units
    str = ""; for c in chars; str *= c; end;
    push!(header, str)
    chars = [Char(c) for c in read(f, Int8, 1024)]  #4 Comments
    str = ""; for c in chars; str *= c; end;
    push!(header, str)
    push!(header, read(f, Int32, 1))                #5 Geometry file
    nz = read(f, Int32, 1)[1]                       #6 Number of zones
    push!(header, nz)
    push!(header, read(f, Int32, 1))                #7 1==structured, 2==unstructured
    push!(header, read(f, Int32, 1))                #8 1==Constant, 2==Periodic, 3==Aperiodic
    ncnt = read(f, Int32, 1)[1]                     #9 Normal centered 1==node, 2==face
    push!(header, ncnt)
    push!(header, read(f, Int32, 1))                #10 1==single, 2==double
    push!(header, read(f, Int32, 1))                #11 iblank
    push!(header, read(f, Int32, 1))                #12 something

    # Open loading file
    if loading_file != nothing
        if loading_file == "automatic"
            # Construct loading file identifier
            iden = filename[1:findfirst("compact.wop", filename)[1]-1]
            iden *= "loading_"

            # Find files that contain the identifier
            loading_files = [fname for fname in readdir(read_path)
                                                       if occursin(iden, fname)]
            # Identify loading file
            if length(loading_files)==0
                @warn("No loading file was found for $(fname).")
                lfile = nothing
            elseif length(loading_files)>1
                @warn("More than one loading file for $(fname) was found."*
                        " The first one will be used: $(loading_files).")
                lfile = loading_files[1]
            else
                lfile = loading_files[1]
            end

        else
            lfile = loading_file
        end

        if lfile == nothing
            loading_file = nothing
        else
            (load_hdr, load_names,
            load_timeinfo, load_dims,
            load_time,
            load_data) = read_wopwoploading(lfile; read_path=read_path,
                                            verbose=verbose, v_lvl=v_lvl+1)
        end
    end

    if verbose
        println("\t"^(v_lvl)*"$(header[1])\t # Magic number")
        println("\t"^(v_lvl)*"$(header[2])\t # Version")
        println("\t"^(v_lvl)*"$(header[3])\t # Units")
        println("\t"^(v_lvl)*"$(header[4])\t # Comment")
        println("\t"^(v_lvl)*"$(header[5])\t # Geometry file")
        println("\t"^(v_lvl)*"$(header[6])\t\t # Number of zones")
        println("\t"^(v_lvl)*"$(header[7])\t # 1==structured, 2==unstructured")
        println("\t"^(v_lvl)*"$(header[8])\t # 1==Constant, 2==Periodic, 3==Aperiodic")
        println("\t"^(v_lvl)*"$(header[9])\t\t # Normal centered 1==node, 2==face")
        println("\t"^(v_lvl)*"$(header[10])\t # 1==single, 2==double")
        println("\t"^(v_lvl)*"$(header[11])\t # iblank")
        println("\t"^(v_lvl)*"$(header[12])\t # something")
    end

    # Error cases:
    Int(header[2][1]) != 1 ? error("Only v1.0 is supported") :
    # Int(header[7][1]) != 2 ? error("Only unstructured patches are supported") :
    # Int(header[8][1]) != 1 ? error("Only constant patches are supported") :
    # Int(header[9][1]) != 2 ? error("Only face-centered patches are supported") :
    Int(header[10][1]) != 1 ? error("Only single-precision floats are supported") :
    nothing;

    if loading_file != nothing
        if load_hdr["structured"] == (Int(header[7][1]) == 1)
            error("Geometry and loading must both be either structured or unstructured!")
        end
    end

    Tflag = Int(header[8][1])       # 1==Constant, 2==Periodic, 3==Aperiodic

    # Read patch headers
    names = []
    nbnodes = []
    nbfaces = []
    imaxs = []
    jmaxs = []
    connectivity = []
    periods = []
    nts = []
    celltypes = []

    for zonei in 1:nz
        # NOTE: Here I assume it is UNSTRUCTURED and constant

        chars = [Char(c) for c in read(f, Int8, 32)]
        str = ""; for c in chars; str *= c; end;

        push!(names, str)
        push!(connectivity, [])

        if Tflag==2
            # Period
            push!(periods, Float64(read(f, Float32)))
        end
        if Tflag==2 || Tflag==3
            # nt
            push!(nts, read(f, Int32))
        end

        # Unstructured case
        if Int(header[7][1]) == 2
            push!(nbnodes, read(f, Int32))
            push!(nbfaces, read(f, Int32))
            if verbose
                println("\t"^(v_lvl)*"$(names[end])\t$(nbnodes[end]) $(nbfaces[end])\t# Name, nbNodes, nbFaces")
            end

            for facei in 1:nbfaces[end]       # Iterate over faces (cells)

                this_nbnodes = read(f, Int32) # Number of node indices to read
                push!(connectivity[end], [])

                for nodei in 1:this_nbnodes   # Collect node indices
                    push!(connectivity[end][end], read(f, Int32))
                end
            end

            push!(celltypes, -1)

        # Structured case
        else
            push!(imaxs, read(f, Int32))
            push!(jmaxs, read(f, Int32))
            push!(nbnodes, imaxs[end]*jmaxs[end])

            if verbose
                println("\t"^(v_lvl)*"$(names[end])\t$(imaxs[end]) $(jmaxs[end])\t# Name, iMax, jMax")
            end

            # Create connectivity: Quadrilateral-faces case
            if imaxs[end]!=1 && jmaxs[end]!=1
                push!(nbfaces, (imaxs[end]-1)*(jmaxs[end]-1))
                for ni in 1:imaxs[end]-1
                    for nj in 1:jmaxs[end]-1
                        # Subscripts of every node
                        subs = ( (ni, nj), (ni+1, nj), (ni+1, nj+1), (ni, nj+1) )
                        # Convert subscripts to linear index of every node
                        this_cell = [Base._sub2ind((imaxs[end], jmaxs[end]), sub...) for sub in subs]
                        push!(connectivity[end], this_cell)
                    end
                end

                push!(celltypes, -1)

            # Create connectivity: Compact-patch (line) case
            else
                push!(nbfaces, imaxs[end]*jmaxs[end]-1)
                for ni in 1:(imaxs[end]*jmaxs[end]-1)
                    push!(connectivity[end], [ni, ni+1])
                end

                push!(celltypes, 4)
            end
        end

        if verbose && Tflag==2
            println("\t"^(v_lvl+1)*"$(periods[end])\t$(nts[end])\t# Period, nt")
        elseif verbose && Tflag==3
            println("\t"^(v_lvl+1)*"$(nts[end])\t# nt")
        end

    end

    # -------------------- READ POINTS AND NORMALS -----------------------------

    Nflag = Int(header[9][1]) # Normal centered 1==node, 2==face

    if Tflag != 1
        points = [zeros(3, nbnodes[zi], nts[zi]) for zi in 1:nz]
        normals = [zeros(3, Nflag==1 ? nbnodes[zi] : nbfaces[zi], nts[zi]) for zi in 1:nz]
        times = [zeros(nts[zi]) for zi in 1:nz]
    else
        points = [zeros(3, nbnodes[zi]) for zi in 1:nz]
        normals = [zeros(3, Nflag==1 ? nbnodes[zi] : nbfaces[zi]) for zi in 1:nz]
        times = []
    end

    # Format cells as vtks (0-indexed)
    vtk_cells = [[Int[ind-1 for ind in cell] for cell in connectivity[zi]] for zi in 1:nz]

    # Zone of every point
    point_zones = [[zi for pi in 1:nbnodes[zi]] for zi in 1:nz]

    str = ""

    # NOTE: I assume that all zones have the same number of time steps, so here
    # I'm checking that assumption
    if Tflag != 1
        for zi in 1:nz-1
            if nts[zi]!=nts[zi+1]
                error("Found that not all zones have the same number of time"*
                        " steps ($(nts[zi])!=$(nts[zi+1])).")
            end
        end
    end

    for ti in 1:(Tflag != 1 ? nts[1] : 1)   # Iterate over time steps
        for zonei in 1:nz                   # Iterate over zones

            # Read time signature
            if Tflag != 1
                times[zonei][ti] = Float64(read(f, Float32))
            end

            # Read points
            for k in 1:3
                for pi in 1:nbnodes[zonei]
                    points[zonei][k, pi, ti] = Float64(read(f, Float32))
                end
            end

            # Read normals
            for k in 1:3
                for ni in 1:(Nflag==1 ? nbnodes[zonei] : nbfaces[zonei])
                    normals[zonei][k, ni, ti] = Float64(read(f, Float32))
                end
            end

            if save_path != nothing

                # Formats points as vtks
                vtk_points = [ points[zonei][:, pi, ti] for pi in 1:nbnodes[zonei] ]

                # Format normals as a vtk field
                vtk_normals = [ normals[zonei][:, pi, ti] for pi in 1:(
                                    Nflag==1 ? nbnodes[zonei] : nbfaces[zonei]) ]

                # Data fields
                point_data = [Dict( "field_name" => "zone",
                                    "field_type" => "scalar",
                                    "field_data" => point_zones[zonei])]

                if Nflag == 2
                    cell_data = [Dict( "field_name"  => "normals",
                                        "field_type" => "vector",
                                        "field_data" => vtk_normals)]
                else
                    cell_data = nothing
                    push!(point_data, Dict( "field_name" => "normals",
                                            "field_type" => "vector",
                                            "field_data" => vtk_normals))
                end

                # Loading data
                if loading_file != nothing

                    # NOTE: Here I assume that the loading file has data for
                    # the same number of zones that the geometry file (even
                    # though that is not necessarely true, but it is for
                    # files generated in my framework)
                    loading_vectors = [load_data[zonei][:, ni, ti] for ni in 1:size(load_data[zonei], 2)]

                    loading_field = Dict(   "field_name"  => load_names[zonei],
                                            "field_type" => "vector",
                                            "field_data" => loading_vectors)

                    if load_hdr["node-centered"]
                        push!(point_data, loading_field)
                    else
                        if cell_data == nothing
                            cell_data = []
                        end
                        push!(cell_data, loading_field)
                    end
                end

                filename = preff*"_"*replace(names[zonei], " " => "")

                gt.generateVTK(filename, vtk_points; cells=vtk_cells[zonei],
                                point_data=point_data, cell_data=cell_data,
                                path=save_path, num=Tflag!=1 ? ti : nothing,
                                override_cell_type=celltypes[zonei])

                if ti==1
                    str *= filename*(Tflag==1 ? ".vtk;" : "...vtk;")
                end
            end

        end
    end

    close(f)

    if paraview && save_path != nothing
        vtk_str *= str
        run(`paraview --data=$(vtk_str)`)
    end

    return names, points, vtk_cells, periods, times, str
end

function save_geomwopwop2vtk(read_path::String, save_path::String; prompt=true,
                                         verbose=true, verbose_level=0, v_lvl=0)
    # Create save path
    gt.create_path(save_path, prompt)

    # Identify loft and compact PSW files
    wopfiles = [fname for fname in readdir(read_path)
                if occursin("loft.wop", fname) || occursin("compact.wop", fname)]

    vtk_str = "$save_path/"

    for fname in wopfiles

        println("\t"^(v_lvl)*"*"^73)
        println("\t"^(v_lvl)*"*\t\tREADING $fname")
        println("\t"^(v_lvl)*"*"^73)

        # Loading file if compact patch
        lfname = occursin("compact.wop", fname) ? "automatic" : nothing

        # Convert from PSW to VTK
        vtk_str *= geomwopwop2vtk(fname; read_path=read_path,
                                        loading_file=lfname,
                                        verbose=verbose && verbose_level>0,
                                        v_lvl=v_lvl+1,
                                        save_path=save_path)[end]
    end

    return vtk_str
end

"""
Reads a time-sequency or single legacy ASCII unstructured VTK file and formats
them into a single file in WOPWOP's geometry format v1.0.
"""
function vtk2wopwop(in_filename, out_filename; read_path="", nums=nothing,
                            save_path="", wopbin=true,
                            t0=0.0, tf=1.0, period=nothing, compact=false)

    _in_filename = in_filename*(in_filename[end-3:end]==".vtk" ? "" : ".vtk")

    # Time type: 1==Constant, 2==periodic, 3==aperiodic
    nums==nothing   ? Tflag = 1 :
    period!=nothing ? Tflag = 2 :
                      Tflag = 3

    # Read first file
    (points, cells, cell_types,
        data) = gt.read_vtk(Tflag==1 ? _in_filename :
                                replace(_in_filename, ".vtk" => ".$(nums[1]).vtk");
                                path=read_path)

    nnodes = size(points, 2)
    ncells = size(cells, 1)

    # Create output file
    f = open(joinpath(save_path, out_filename), "w")

    # Binary / ASCII printing
    prnt(x) = wopbin ? write(f, x) : print(f, x)
    prntln(x) = wopbin ? write(f, x) : print(f, x, "\n")

    # Convertion to 4-bytes numbers
    # NOTE: 4 bytes = 4*8 bites = 32 bites
    fl(x) = Float32(x)
    nt(x) = Int32(x)
    # Convertion to n-bytes string
    st(x::String, n) = x * " "^(n-length(x))

    # Magic number
    prntln(nt(42))
    # Version number
    prnt(nt(1))
    prntln(nt(0))
    # Units for Tecplot
    prntln(st("N/m^2", 32))
    # Comments
    prntln(st("Geometry input file for PSU-WOPWOP (Format v1.0)\n"*
              "------------------------------------------------\n"*
              "Created by FLOWNoise (written by Eduardo Alvarez)\n"*
              "https://github.com/byuflowlab/FLOWNoise\n"*
              "Creation date: $(Dates.now())\n"*
              "Units: SI\n"*
              (
                compact ? "Format: Structured grid, node-centered" :
                          "Format: Unstructured grid, face-centered"
              ), 1024))

    # Format string
    prntln(nt(1))               # Geometry file flag
    prntln(nt(1))               # Number of zones
    prntln(nt(2^!compact))      # 1==structured, 2==unstructured
    prntln(nt(Tflag))           # Geometry 1==constant, 2==periodic, 3==aperiodic
    prntln(nt(2^!compact))      # Normal vectors 1==node, 2==face
    prntln(nt(1))               # Floating point 1==single, 2==double
    prntln(nt(0))               # iblank values 1==included, 0==not
    prntln(nt(0))               # WOPWOP secret conspiracy

    # ----------------- PATCH ----------------------------------------
    # Name
    prntln(st(compact ? "vtkcompactpatch" : "vtkpatch", 32))
    # TimeInformation
    if Tflag==2
        # Period
        prntln(fl(period))
        # nt
        prntln(nt( length(nums) ))
    elseif Tflag==3
        # nt
        prntln(nt( length(nums) ))
    end

    if compact
        # iMax
        prntln(nt( 1 ))
        # jMax
        prntln(nt( nnodes ))
    else
        # nbNodes
        prntln(nt( nnodes ))
        # nbFaces
        prntln(nt( ncells ))
        # Connectivity
        # NOTE: Here I assume that connectivity is the same in all VTK files
        for cell in cells
            prnt(nt( size(cell, 1) ))         # Number of nodes in this cell
            # for pi in reverse(cell)         # Clockwise node ordering
            for pi in cell                    # NOTE: Turns out that FLOWVLM's VTKs are already clockwise
                prnt(nt( pi+1 ))              # 1-indexed node index
            end
            if !wopbin; prntln(""); end;
        end
    end

    # ----------------- DATA OF PATCH -------------------------------

    Ns = zeros(3, compact ? nnodes : ncells)

    for ti in 1:(Tflag==1 ? 1 : length(nums))
        # Read vtk
        if ti != 1
            # NOTE: Here I assume that `nums` is given sequentally
            (points, cells, cell_types,
             data) = gt.read_vtk(replace(_in_filename, ".vtk" => ".$(nums[ti]).vtk");
                                    path=read_path)

            Ns[:] .= 0
        end

        # Calculate normals
        if compact
            # NOTE: Here I assume that the lifting line is defined sequentially
            for pi in 1:nnodes-1
                # NOTE: Here I give it the vector that connect the points
                # as the normal, crossing my fingers that PSU-WOPWOP will never
                # use the actual normal direction for anything (I think that
                # it only expects the magnitude to be the length of this segment
                # of the line)
                D = points[:, pi+1] - points[:, pi]

                # Here I split the length in two for each node
                Ns[:, pi] += D/2
                Ns[:, pi+1] += D/2
            end

            # # Here I'm normalizing the normal by the length
            # for pi in 1:nnodes
            #     Ns[:, pi] /= norm(Ns[:, pi])
            # end

        else
            for (ci, cell) in enumerate(cells)

                if length(cell) != 4
                    error("Only quadrilateral cells are supported!")
                end

                crss1 = -cross( points[:, cell[2]+1] - points[:, cell[1]+1],
                                points[:, cell[3]+1] - points[:, cell[1]+1])
                crss2 = -cross( points[:, cell[4]+1] - points[:, cell[3]+1],
                                points[:, cell[1]+1] - points[:, cell[3]+1])

                Ns[:, ci] .= crss1/2 + crss2/2
            end
        end

        # Time signature
        if Tflag != 1
            prntln(fl( (tf-t0)*(ti-1)/(length(nums)-1) ))
        end

        # nbNodes floating point x coordinates
        # nbNodes floating point y coordinates
        # nbNodes floating point z coordinates
        for k in 1:3
            for pi in 1:nnodes
                prntln(fl(points[k, pi]))
            end
        end

        if compact
            # imax × jmax floating point normal vector x coordinates
            # imax × jmax floating point normal vector y coordinates
            # imax × jmax floating point normal vector z coordinates
            for k in 1:3
                for j in 1:nnodes
                    prntln(fl(Ns[k, j]))
                end
            end

        else
            # nbFaces floating point normal vector x coordinates
            # nbFaces floating point normal vector y coordinates
            # nbFaces floating point normal vector z coordinates
            for k in 1:3
                for j in 1:ncells
                    prntln(fl(Ns[k, j]))
                end
            end
        end
    end

    close(f)

    return nothing
end

"""
Read the solution field created by PSU-WOPWOP named `fieldname` (i.e., pressure,
spl_spectrum, OASPLdB, or OASPLdBA) under `read_path`.
"""
function read_pswfield(fieldname::String, read_path; re20=false, tec=false,
                                                                     optargs...)
    header, field = read_wopwopoutput(fieldname; read_path=read_path, tec=tec,
                                                                    optargs...)

    if re20
        header_hash = Dict( (occursin(",re20", h) ? h[1:findfirst(",re20", h)[1]-1] : h, i)
                                                            for (i,h) in enumerate(header))
    else
        header_hash = Dict((h, i) for (i,h) in enumerate(header))
    end

    return field, header, header_hash
end

function read_pswfield(fieldnames::Array{String, 1}, read_path;
                            re20crit=["OASPLdB", "OASPLdBA"], tec=false,
                            optargs...)

    datas = Dict()

    for fieldname in fieldnames
        data = read_pswfield(fieldname, read_path; re20=fieldname in re20crit,
                                                       tec=tec, optargs...)

        datas[fieldname] = Dict("field"=>data[1], "hd"=>data[2], "hs"=>data[3])
    end

    return datas
end

function fetch_pswfield(read_path::String, args...;
                        # fieldnames=["pressure", "spl_spectrum", "OASPLdB", "OASPLdBA"],
                        fieldnames=["spl_spectrum", "OASPLdB", "OASPLdBA"],
                        datasets=def_datasets_psw, verbose=true, v_lvl=0,
                        verbose_level=1, optargs...)

    if verbose; println("\t"^v_lvl*
                        "*"^72*"\n*\tReading dataset $read_path\n"*"*"^72); end;

    datasets[read_path] = read_pswfield(fieldnames, read_path, args...;
                                            verbose=verbose,
                                            verbose_level=verbose_level-1,
                                            v_lvl=v_lvl+1,
                                            optargs...)

    return datasets[read_path]
end

function fetch_pswdataset(read_path, args...; datasets=def_datasets_psw,
                                                verbose=true, v_lvl=0,
                                                force_read=false,
                                                optargs...)

    if !(read_path in keys(datasets)) || force_read
        return fetch_pswfield(read_path, args...; verbose=verbose, v_lvl=v_lvl,
                                                  datasets=datasets, optargs...)
    else
        if verbose; println("\t"^v_lvl*
                            "Dataset $(read_path) has already been loaded."*
                            " It will not be read again."); end;
        return datasets[read_path]
    end

end

# ------------------------- FUNCTIONS FOR BPM ----------------------------------
function read_bpmoutput(args...; datasets=def_datasets_bpm, optargs...)
    return read_wopwopoutput(args...; tec=true, datasets=datasets, optargs...)
end
function fetch_bpmdataset(args...; datasets=def_datasets_bpm, optargs...)
    return fetch_pswdataset(args...; tec=true, datasets=datasets,
                                     fieldnames=["spl_spectrum", "splA_spectrum",
                                                          "OASPLdB", "OASPLdBA"],
                                                                    optargs...)
end

max_stars_repo_name: byuflowlab/FLOWNoise

max_issues_repo_name: byuflowlab/FLOWNoise

max_forks_repo_name: byuflowlab/FLOWNoise


text:
g = NEAT.Global(config) # reset globals
s1 = NEAT.Species(g, NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene))
@test s1.id == 1
@test length(s1) == 1
@test s1.representant == s1.subpopulation[1]
@test s1.representant.species_id == 1

NEAT.add(s1, NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene))
@test s1.id == 1
@test length(s1) == 2
@test s1.subpopulation[2].species_id == 1

# tests of species id increment && of no parameter sharing
s2 = NEAT.Species(g, NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene))
@test s2.id == 2
@test length(s2) == 1
@test s2.representant == s2.subpopulation[1]
@test s2.representant.species_id == 2

# setup
g = NEAT.Global(config) # reset globals
s3 = NEAT.Species(g, NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene))
NEAT.add(s3, NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene))
NEAT.add(s3, NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene))
s3.subpopulation[1].fitness = 1.
s3.subpopulation[2].fitness = 2.
s3.subpopulation[3].fitness = 3.

# avg fitness & no_improvement_age tests
@test NEAT.average_fitness(s3) == 2.
map(x->NEAT.average_fitness(s3),1:100) # run 100 more times wiht no change
@test s3.no_improvement_age == 100

# tournament selection test
best = NEAT.tournamentSelection(s3)
@test all(ch->ch.id !=1, map(i -> NEAT.tournamentSelection(s3),[1:100]))

# reproduce() tests
g = NEAT.Global(config) # reset globals
s = NEAT.Species(g,NEAT.create_fully_connected(g))
NEAT.add(s,NEAT.create_fully_connected(g))
# NEAT.add(s,NEAT.create_fully_connected(g))
# NEAT.add(s,NEAT.create_fully_connected(g))
spawn_amount = 3
s.spawn_amount = spawn_amount
offspring = NEAT.reproduce(g,s)
@test length(offspring) == spawn_amount
@test offspring[1].id == 3
@test offspring[2].id == 4
@test offspring[3].id == 5


text:
abstract type Gravity end

get_gravity_horizon(grav::Gravity) = grav.hor
get_gravity_accel(grav::Gravity) = norm(get_gravity_horizon(grav))
get_gravity_body(grav::Gravity, att::Attitude) =
    hor2body(get_gravity_horizon(grav)..., get_quaternions(att)...)

"""
    EarthConstantGravity(value, body_vector)

Earth Constant gravity model with value 9.80665 m/s² and oriented towards
local horizon Z axis.

"""
struct EarthConstantGravity<:Gravity
    hor::Array{T, 1} where T<:Number
end

EarthConstantGravity() = EarthConstantGravity([0, 0, GRAVITY_ACCEL])

calculate_gravity(grav::EarthConstantGravity, pos::Position) = grav


text:
module DojoEnvironments

using LinearAlgebra
using Random

using MeshCat
using StaticArrays
using Dojo

import Dojo: string_to_symbol, add_limits, RotX, RotY, RotZ, build_robot, initialize!, set_minimal_coordinates_velocities!, velocity_index, vector_rotate, Prototype, rotation_vector, gray_light

export
    Environment,
    dynamics,
    dynamics_jacobian_state,
    dynamics_jacobian_input,
    get_environment,
    get_mechanism,
    get_observation,
    cost,
    is_done,
    render,
    seed,
    close,
    Space,
    BoxSpace

include("mechanisms.jl")
include("environment.jl")
include("dynamics.jl")
include("utilities.jl")
include("include.jl")

end # module


text:
using StaticArrays: SVector
using LinearAlgebra: dot, cross, norm

# ---------------------------------------------------------------------------------------- #
# EXPR UTILITIES

# extract all the `Symbol`s referenced in an `Expr`
function get_expr_symbols!(syms::Vector{Symbol}, ex::Expr)
    for arg in ex.args
        if arg isa Symbol
            arg ∈ syms || push!(syms, arg)
        elseif arg isa Expr
            get_expr_symbols!(syms, arg)
        end
    end
    return syms
end
get_expr_symbols(ex::Expr) = get_expr_symbols!(Symbol[], ex)

# ---------------------------------------------------------------------------------------- #
# LOAD CONSTANT GLOBAL `pointsd_Nd`, & `paramsd_Nd` DICTIONARIES (for N={2,3})
include(joinpath(pkgdir(KPaths), "assets", "data-SeeK.jl"))
include(joinpath(pkgdir(KPaths), "assets", "data-2d.jl"))

# ---------------------------------------------------------------------------------------- #
# CODEGEN TO CREATE FUNCTIONS FOR EXTENDED BRAVAIS TYPES

# Generate method definitions `$(bt)_points` (`bt` denotes the extended Bravais type) for
# each dataset using metaprogramming. If the datasets require arguments, the arguments are
# listed in the order :a, :b, :c, :β and includes only the featured set of arguments.
for (D, pointsd, paramsd) in zip((2,3), (pointsd_2d, pointsd_3d), (paramsd_2d, paramsd_3d))
    for (bt, points) in pointsd
        fn = Symbol(string(bt)*"_points")
        params = get(paramsd, bt, nothing)

        # --- build up an expression for the dictionary of k-point labels and vectors ---
        ex = :(Dict())
        for (klab, point_expr) in points
            point_expr′ = Expr(:call, SVector{D,Float64}, point_expr.args...) # convert to SVector
            push!(ex.args, Expr(:call, :(=>), QuoteNode(klab), point_expr′))
        end

        if params === nothing
            # --- only constant points; define "constant" function of `Rs` ---
            @eval $fn(Rs::Union{Nothing, AbstractVector{<:SVector{$D, <:Real}}}) = $ex

        else
            # --- free arguments; define constant dependent on `Rs` ---

            # write out the all the parameter definitions as expressions
            params_ex = quote end
            for (paramlab, param_expr) in params
                push!(params_ex.args, Expr(:(=), :($(paramlab)), param_expr))
            end
            # find out what the arguments are; can only be {:a, :b, :c} and {:sinβ, :cosβ}
            syms = get_expr_symbols(params_ex) # get all referenced symbol in `params_ex`
            # create list of typed arguments, e.g. (:(a::Real), :(b::Real), ...)
            setup_ex = quote end
            for (i, sym) in enumerate((:a, :b, :c))
                if sym ∈ syms
                    push!(setup_ex.args, Expr(:(=), sym, :(norm(Rs[$i]))))
                end
                #sym ∈ syms && push!(args, Expr( :(::), sym, :Real))
            end
            # define "cosβ = cos(β)" and "sinβ = sin(β)" if we can find :cosβ or :sinβ in the
            # argument list; in practice, they only ever occur together in SeeK, so we can just
            # check for one and add both.
            if :cosβ ∈ syms
                for (i, sym) in zip((1,3), (:a, :c))
                    # if :a or :c are not in argslist, then we need to add them, since we need
                    # them to compute :cosβ and :sinβ as dot and cross products; β is ∠(R₃, R₁)
                    sym ∉ syms && push!(setup_ex.args, Expr(:(=), sym, :(norm(Rs[$i]))))
                end
                push!(setup_ex.args, Expr(:(=), :cosβ, :(dot(Rs[3], Rs[1])/(c*a))))
                push!(setup_ex.args, Expr(:(=), :sinβ, :(sqrt(1 - cosβ^2)))) # acos(sin(x)) = √(1-x²)
            end

            # --- make function of the featured subset of :a, :b, :c, and :β, in that order ---
            @eval function $fn(Rs::AbstractVector{<:SVector{$D, <:Real}})
                $setup_ex
                $params_ex
                return $ex
            end
        end

        #=
        # see what methods we generated and check that they "work":
        fn′ = getfield(Main, fn) # get the actual function, not just its symbol
        mth = only(methods(fn′))
        println("\n", mth)
        display(fn′([rand(SVector{D,Float64}) for _ in 1:3]))
        =#
    end
end

# ---------------------------------------------------------------------------------------- #
# CODEGEN FOR BRANCHTABLE TO ABOVE FUNCTIONS

# generate a branch-table function (just a lot of if-statements, that simply dispatch to
# `($ext_bt)_points(Rs)`, depending on the value of `ext_bt`:
for (D, pointsd) in zip((2,3), (pointsd_2d, pointsd_3d))

    branchtable = Expr(:if)
    let current = branchtable
        for (i,ext_bt) in enumerate(keys(pointsd))
            if i ≠ 1
                push!(current.args, Expr(:elseif))
                current = current.args[end]
            end
            push!(current.args, :(ext_bt == $(QuoteNode(ext_bt))))
            fn = Symbol(string(ext_bt)*"_points")
            push!(current.args, :($fn(Rs)))
        end
        # the final entry in the :if call block is implicitly an `else` block:
        push!(current.args, :(throw(DomainError(ext_bt, "invalid extended Bravais type"))))
    end

    f = Symbol(:get_points_, D, :d)
    @eval begin
        @doc """
            $($f)(ext_bt, Rs)

        Return the labels and points in the $($D)D **k**-path associated with an extended
        Bravais type `ext_bt` and a (conventional) direct basis `Rs` (can be `nothing` if
        there is no dependence on the basis) as a `Dict{Symbol, SVector{$($D), Float64}}`.
        """
        function $f(ext_bt::Symbol,
                    Rs::Union{Nothing, AVec{<:AVec{<:Real}}})
            $branchtable
        end
    end
end


text: compute(n::Int)::Int = (n - 5) ÷ 7 * 3 + 2

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
shape_map(x,upper,lower) = min(max(x,lower), upper)

"""
    reshape_values(xvals::Vector{<:Real}, yvals::Vector{<:Real}; increasing::Bool = true,
                       concave::Bool = true, shape_map::Function = shape_map)
This reshapes a vector of yvalues. For instance if we are doing fixed point acceleration that
should result in a monotonic concave function (ie the consumption smoothing problem from the documentation examples)
then we may end up with occasional non monotonic/concave values due to a dodgy optimiser or some other numerical issue.
So we can use reshape_values to adjust the values that cannot be true.
### Inputs
* `xvals` - A vector of x coordinates.
* `yvals` - A vector of y coordinates
* `increasing` - Should the y values be increasing. If false then they must be decreasing
* `concave` - Should the y values be concave. If false then they must be convex
* `shape_map` - A function used to adjust values to be increasing-concave (or whatever settings)
### Returns
* An updated vector of y values.
"""
function reshape_values(xvals::Vector{<:Real}, yvals::Vector{<:Real}; increasing::Bool = true,
                        concave::Bool = true, shape_map::Function = shape_map)
    lenlen = length(yvals)
    new_shape_vec = Array{Float64,1}(undef, lenlen)
    if lenlen > 0 new_shape_vec[1] = yvals[1] end
    if lenlen > 1
        upper = increasing ? Inf : new_shape_vec[1]
        lower = increasing ? new_shape_vec[1] : -Inf
        new_shape_vec[2] = shape_map(yvals[2],upper, lower)
    end
    for i in 3:lenlen
        stepp = xvals[i] - xvals[i-1]
        previous_grad = (new_shape_vec[i-1]-new_shape_vec[i-2])/(xvals[i-1] - xvals[i-2])
        upper =  new_shape_vec[i-1] + (increasing ? (concave ? previous_grad * stepp : Inf) : (concave ? previous_grad * stepp : 0))
        lower =  new_shape_vec[i-1] + (increasing ? (concave ? 0 :  previous_grad * stepp) : (concave ? -Inf : previous_grad * stepp))
        new_shape_vec[i] = shape_map(yvals[i],upper, lower)
    end
    return new_shape_vec
end


text:
using SnoopCompile

println("tests infer benchmark")

@snoopi_bench "Juno"


text:
module ProtStructs

using StaticArrays
import LinearAlgebra: norm, ⋅, ×
import Chemfiles    ## don't bring out names for now

# Write your package code here.
include("constants.jl")

include("Model.jl")
include("Parser.jl")
include("Assignments.jl")

include("GeomUtils.jl")
include("FrameGeom.jl")
# something on contact maps/solvent accessibility

end


text:
using DataFrames

"""
    THD(u)

Calculate the Total Harmonic Distortion at all nodes.

THD is a measure for the amount of voltage distortion present at a node.
This function calculates two alternative definitions

- THD\\_F: relative to fundamental, more commonly used and
- THD\\_R: RMS, relative to all frequencies
- THD\\_40: relative to fundamental, only until h = 40 (or 39) as used by
    standard EN 61000.
"""
function THD(u)
    harmonics = sort(collect(keys(u)))
    THD = DataFrame(THD_F = zeros(size(u[1], 1)),
                    THD_R = zeros(size(u[1], 1)),
                    THD_40 = zeros(size(u[1], 1)))
    for ID in 1:size(u[1], 1)
        THD[ID, "THD_F"] = sqrt(sum([u[h].v[ID]^2 for h in harmonics[2:end]]))./u[1].v[ID]
        THD[ID, "THD_R"] = sqrt(sum([u[h].v[ID]^2 for h in harmonics[2:end]]))./sqrt(sum([u[h].v[ID]^2 for h in harmonics]))
        harmonics_40 = harmonics[harmonics.<=40]
        THD[ID, "THD_40"] = sqrt(sum([u[h].v[ID]^2 for h in harmonics_40[2:end]]))./u[1].v[ID]
    end
    return THD
end

"""
    limits(u, nodeID)

Check if harmonic limits according to EN 50160 and EN 61000 are met at node.
"""
function limits(u, nodeID)
    h = [3, 5, 7, 9, 11, 13, 15]
    v = [u[i].v[nodeID] for i in h]
    limit_EN50160 = [0.05, 0.06, 0.05, 0.015, 0.035, 0.03, 0.005]
    limit_EN61000 = [0.05, 0.06, 0.05, 0.015, 0.035, 0.03, 0.003]
    met_EN50160 = collect(v .< limit_EN50160)
    met_EN61000 = collect(v .< limit_EN61000)

    if all(met_EN50160) && all(met_EN61000)
        println("Harmonic limits met at node "*string(nodeID))
    else
        println("Harmonic limits not met at node "*string(nodeID)*"! Check DataFrame for details.")
    end

    THD_40 = HarmonicPowerFlow.THD(u).THD_40[nodeID]

    if THD_40 < 0.08
        println("THD_40 at node "*string(nodeID)*" within limits (< 8%).")
    else
        println("THD_40 at node "*string(nodeID)*" exceeds limits (> 8%)!")
    end

    DataFrame(v=v, within_EN50160 = met_EN50160, within_EN61000=met_EN61000)
end

# function TDD()

# function PF()


text:
function basictest(lib::PolyhedraLibrary)
    @testset "Polyhedron eltype conversion tests with $(typeof(lib))" begin
        A = [1 0; 0 1]
        b = [0, 0]
        linset = BitSet([2])
        R = [1 1]
        h = hrep(A, b, linset)
        v = vrep(A, R)

        p = copy(polyhedron(h, lib))
        @test hrepiscomputed(p)
        @test !vrepiscomputed(p)
        p = Polyhedron{2, Float64}(p)
        @test hrepiscomputed(p)
        @test !vrepiscomputed(p)
        inequality_fulltest(hrep(p), A, b, linset)
        p = copy(polyhedron(v, lib))
        @test !hrepiscomputed(p)
        @test vrepiscomputed(p)
        p = Polyhedron{2, Float64}(p)
        @test !hrepiscomputed(p)
        @test vrepiscomputed(p)
        generator_fulltest(vrep(p), A, R)
    end
    @testset "Polyhedron cartesian product tests with $(typeof(lib))" begin
        v1 = vrep([1 2; 3 4], [5 6])
        v2 = vrep([5 6; 7 8], [9 1; 2 3])
        v = v1 * v2
        expv = vrep([1 2 5 6; 1 2 7 8; 3 4 5 6; 3 4 7 8], [5 6 0 0; 0 0 9 1; 0 0 2 3])
        generator_fulltest(v, expv)
        pv1 = polyhedron(v1, lib)
        pv2 = polyhedron(v2, lib)
        qv = pv1 * pv2
        generator_fulltest(vrep(qv), expv)
        h1 = hrep([1 2; 3 4], [5, 6])
        h2 = hrep([7 8; 9 1], [2, 3])
        h = h1 * h2
        exph = hrep([1 2 0 0; 3 4 0 0; 0 0 7 8; 0 0 9 1], [5, 6, 2, 3])
        inequality_fulltest(h, exph)
        ph1 = polyhedron(h1, lib)
        ph2 = polyhedron(h2, lib)
        qh = ph1 * ph2
        inequality_fulltest(hrep(qh), exph)
    end
end


text:
using MyPkg
using Documenter

DocMeta.setdocmeta!(MyPkg, :DocTestSetup, :(using MyPkg); recursive=true)

makedocs(;
    modules=[MyPkg],
    authors="Gudongyangg",
    repo="https://github.com/Gudongyangg/MyPkg.jl/blob/{commit}{path}#{line}",
    sitename="MyPkg.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Gudongyangg.github.io/MyPkg.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Gudongyangg/MyPkg.jl",
    devbranch="master",
)


text:
# TODO: handle multidimensional case

"""
    product_integral(f, g, [, solver])

Compute the integral of `f .* g` over [0,1].
"""
function product_integral(f, g, solver=HCubatureJL())
    function fun(u,_)
        f(u) .* g(u)
    end
    prob = QuadratureProblem{false}(fun, 0, 1, ())
    solve(prob, solver).u
end

"""
    kernel_integral(k, f, t [, solver])

Compute the kernel integral `∫ k(t,s)f(s)ds`.
"""
function kernel_integral(k, f, t, solver=HCubatureJL())
    fun(u) = k(t,u)
    return product_integral(fun, f, solver)
end

"""
    kernel_transpose_integral(k, f, t [, solver])

Compute the kernel integral `∫ k(s,t)f(s)ds`.
"""
function kernel_transpose_integral(k, f, t, solver=HCubatureJL())
    fun(u) = k(u,t)
    return product_integral(fun, f, solver)
end

"""
    kernel_double_integral(k, f, g, [, solver])

Compute the kernel integral `∫∫ k(t,s)f(t)g(s) ds dt`.
"""
function kernel_double_integral(k, f, g, solver=HCubatureJL())
    function fun(u,_)
        f(u[1]) .* k(u[1], u[2]) .* g(u[2])
    end
    prob = QuadratureProblem{false}(fun, [0,0], [1,1], ())
    solve(prob, solver).u
end


text:
module Spectra

export spectrum

# AbstractSpectrum and common functionality
include("common.jl")

# Spectrum types and basic arithmetic
include("Spectrum.jl")
include("EchelleSpectrum.jl")

"""
    spectrum(wave, flux; kwds...)

Construct a spectrum given the spectral wavelengths and fluxes. This will automatically dispatch the correct spectrum type given the shape and element type of the given flux. Any keyword arguments will be accessible from the spectrum as properties.

# Examples
```jldoctest
julia> wave = range(1e4, 4e4, length=1000);

julia> flux = 100 .* ones(size(wave));

julia> spec = spectrum(wave, flux)
Spectrum(Float64, Float64)

julia> spec = spectrum(wave, flux, name="Just Noise")
Spectrum(Float64, Float64)
  name: Just Noise

julia> spec.name
"Just Noise"

```

There is easy integration with [Unitful.jl](https://github.com/painterqubits/unitful.jl) and its sub-projects and [Measurements.jl](https://github.com/juliaphysics/measurements.jl)

```jldoctest
julia> using Unitful, UnitfulAstro, Measurements

julia> wave = range(1, 4, length=1000)u"μm";

julia> sigma = randn(size(wave));

julia> flux = (100 .± sigma)u"erg/cm^2/s/angstrom";

julia> spec = spectrum(wave, flux)
Spectrum(Quantity{Float64,𝐋,Unitful.FreeUnits{(μm,),𝐋,nothing}}, Quantity{Measurement{Float64},𝐌*𝐋^-1*𝐓^-3,Unitful.FreeUnits{(Å^-1, erg, cm^-2, s^-1),𝐌*𝐋^-1*𝐓^-3,nothing}})

```

For a multi-order spectrum, all orders must have the same length, so be sure to pad any ragged orders with NaN.

```jldoctest
julia> wave = reshape(range(100, 1e4, length=1000), 100, 10)';

julia> flux = ones(10, 100) .* collect(1:10);

julia> spec = spectrum(wave, flux)
EchelleSpectrum(Float64, Float64)

```
"""
function spectrum(wave::AbstractVector{<:Real}, flux::AbstractVector{<:Real}; kwds...)
    @assert size(wave) == size(flux) "wave and flux must have equal size"
    Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
end

function spectrum(wave::AbstractVector{<:Quantity}, flux::AbstractVector{<:Quantity}; kwds...)
    @assert size(wave) == size(flux) "wave and flux must have equal size"
    @assert dimension(eltype(wave)) == u"𝐋" "wave not recognized as having dimensions of wavelengths"
    Spectrum(wave, flux, Dict{Symbol,Any}(kwds))
end

function spectrum(wave::AbstractMatrix{<:Real}, flux::AbstractMatrix{<:Real}; kwds...)
    @assert size(wave) == size(flux) "wave and flux must have equal size"
    EchelleSpectrum(wave, flux, Dict{Symbol,Any}(kwds))
end

# tools
include("utils.jl")
include("transforms/transforms.jl")
include("plotting.jl")
include("fitting/fitting.jl")

end # module


text:
import OldOnly1
import ScientificTypes
import Test

Test.@testset "OldOnly1.jl" begin
    Test.@test OldOnly1.f(1) == 2
end


text:
"""
    check_method_ambiguity_binary(op;
                                  [print_results]::Bool=false,
                                  [print_warnings::Bool=print_results])::Bool

Check that a given binary operation does not have method ambiguities for some
combination of set types.

### Input

- `op`             -- binary set function name
- `print_results`  -- (optional, default: `false`) flag for printing
                      intermediate results
- `print_warnings` -- (optional, default: `print_results`) flag for printing
                      warnings
- `ignore_errors`  -- (optional, default: `true`) flag for returning `false` if
                      any error occurred

### Output

`true` iff there is no method ambiguity for any set type combination.

### Notes

It is not required that every set type combination is implemented.
The check is only applied to implemented combinations.

This function can also print all intermediate results to STDOUT.
Both method ambiguities/other errors and missing implementations are then
reported as warnings.
Optionally only the warnings can be printed.

### Examples

```julia
julia> check_method_ambiguity_binary(issubset)
```
"""
function check_method_ambiguity_binary(op;
                                       print_results::Bool=false,
                                       print_ambiguity_warnings::Bool=true,
                                       print_warnings::Bool=print_results,
                                       ignore_errors::Bool=true)::Bool
    H1 = LinearConstraint([1.], 1.)
    polytope_constraints_1D = [H1,
                               LinearConstraint([-1.], 0.)]
    sets_1D = LazySet{Float64}[
        Ball1(zeros(1), 1.),
        Ball2(zeros(1), 1.),
        BallInf(zeros(1), 1.),
        Ballp(1.5, zeros(1), 1.),
        Ellipsoid(hcat([1.])),
        EmptySet(2),
        H1,
        HPolytope(polytope_constraints_1D),
        HPolyhedron([H1]),
        Hyperplane(ones(1), 1.),
        Hyperrectangle(zeros(1), ones(1)),
        Interval(0., 1.),
        Singleton(zeros(1)),
        VPolytope([zeros(1)]),
        Universe(1),
        ZeroSet(1),
        Zonotope(zeros(1), [1. 2.;])
    ]
    H2 = LinearConstraint(ones(2), 1.)
    polytope_constraints_2D = [H2,
                               LinearConstraint([-1., 0.], 0.),
                               LinearConstraint([0., -1.], 0.)]
    sets_2D = LazySet{Float64}[
        Ball1(zeros(2), 1.),
        Ball2(zeros(2), 1.),
        BallInf(zeros(2), 1.),
        Ballp(1.5, zeros(2), 1.),
        Ellipsoid([1. 0.; 0. 1.]),
        EmptySet(2),
        H2,
        HPolygon(polytope_constraints_2D),
        HPolygonOpt(polytope_constraints_2D, 1),
        HPolytope(polytope_constraints_2D),
        HPolyhedron([H2]),
        Hyperplane(ones(2), 1.),
        Hyperrectangle(zeros(2), ones(2)),
        Line(ones(2), 1.),
        LineSegment(zeros(2), ones(2)),
        Singleton(zeros(2)),
        Universe(2),
        VPolygon([zeros(2)]),
        VPolytope([zeros(2)]),
        Universe(2),
        ZeroSet(2),
        Zonotope(zeros(2), [1. 0.; 0. 1.])
    ]

    if startswith("$op", "LazySets.")
        ops = "$op"[10:end] # remove prefix 'LazySets.'
    else
        ops = op
    end
    has_ambiguities = false
    has_other_problems = false
    print_results && println("testing operation $ops")
    for types in [sets_1D, sets_2D]
        for i1 in types
            t1 = typeof(i1)
            print_results && println("type 1: $t1")
            for i2 in types
                t2 = typeof(i2)
                print_results && println("type 2: $t2")
                try
                    op(i1, i2)
                    print_results && println("test succeeded")
                catch ex
                    # print error message
#                     showerror(STDOUT, ex); println()

                    # get error message
                    buf = IOBuffer()
                    showerror(buf, ex)
                    msg = String(take!(buf))

                    if occursin("is ambiguous", msg)
                        # ambiguity error
                        has_ambiguities = true
                        print_warnings && @warn "ambiguous method " *
                            "'$ops(::$(typeof(i1)), ::$(typeof(i2)))' found"
                        print_ambiguity_warnings && @warn msg
                    elseif occursin("no method matching", msg)
                        # not implemented
                        print_results && println("no method " *
                            "'$ops(::$(i1), ::$(i2)' exists")
                    else
                        # other error
                        has_other_problems = true
                        if print_warnings
                            @warn "unknown error message:"
                            @warn msg
                        end
                    end
                end
            end
        end
    end
    result = !has_ambiguities
    if !ignore_errors
        result &= !has_other_problems # also report other errors
    end
    return result
end


text:
using Test
using Hyperspecialize
push!(LOAD_PATH, ".")

@testset "Simple Concretization Pregame" begin
  # Concretization is rather persistent
  @test (@concretize Pregame String) == Set{Type}([String])
end

include("simple_concretization.jl")
include("module_concretization.jl")
include("simple_replicable.jl")
include("module_replicable.jl")


text:
using Logging

macro transition_tests_state(name)
    return :(
        mutable struct $name <: HSM.AbstractHsmState
            state_info::HSM.HsmStateInfo
            entered::Bool
            exited::Bool

            $name(parent_state) = new(HSM.HsmStateInfo(parent_state), false, false)
        end
    )
end

function reset_transition_test_state!(s)
    s.entered = false
    s.exited = false
end


text:
# This file is a part of JuliaPackageTemplate.jl, licensed under the MIT License (MIT).

import Test
import JuliaPackageTemplate
import Documenter

Test.@testset "Package JuliaPackageTemplate" begin
    include("test_hello_world.jl")

    # doctests
    Documenter.DocMeta.setdocmeta!(
        JuliaPackageTemplate,
        :DocTestSetup,
        :(using JuliaPackageTemplate);
        recursive=true,
    )
    Documenter.doctest(JuliaPackageTemplate)
end # testset


text:
using BenchmarkTools

function to_pari(x::fmpz_mat, varname = "x")
  s = "$varname = ["
  y = Hecke._eltseq(x)
  n = ncols(x)
  for i in 1:length(y)
    s = s * "$(y[i])"

    if i % n != 0
      s = s * ","
    else
      if i != length(y)
        s = s * ";"
      end
    end
  end

  s = s * "]"
  return s
end

function to_magma(x::fmpz_mat, varname = "x")
  n = nrows(x)
  m = ncols(x)
  s = "$varname := Matrix(Integers(), $n, $m, ["
  y = Hecke._eltseq(x)
  for i in 1:length(y)
    s = s * "$(y[i])"
    if i < length(y)
      s = s * ","
    end
  end
  s = s * "]);"
end

function to_magma(x::fmpq_mat, varname = "x")
  n = nrows(x)
  m = ncols(x)
  s = "$varname := Matrix(Rationals(), $n, $m, ["
  y = Hecke._eltseq(x)
  for i in 1:length(y)
    s = s * "$(y[i])"
    if i < length(y)
      s = s * ","
    end
  end
  s = s * "]);"
  return s
end

function to_sage(x::fmpq_mat, varname = "x")
  n = nrows(x)
  m = ncols(x)
  s = "$varname = Matrix(QQ, $n, $m, ["
  y = Hecke._eltseq(x)
  for i in 1:length(y)
    s = s * "$(y[i])"
    if i < length(y)
      s = s * ","
    end
  end
  s = s * "]);"
  return s
end

function benchmark_code_sage(L::Vector{Hecke.ZLat}, repetitions::Int = 10; out_file::String = "sage_timings")
  s = """
  from sage.all import *
  import timeit
  res = []
  out_file = \"$(out_file)\"
  """
  for i in 1:length(L)
    G = gram_matrix(L[i])
    s = s * to_sage(G, "G$i") * "\n"
    s = s * """
    print(\"Doing\", $(i))
    t = 0.0
    for i in range($repetitions):
      L = IntegralLattice(G$(i))
      s = timeit.repeat(lambda: L.automorphisms(), repeat = 1, number = 1)
      t = t + s[0]
    t = t/$(repetitions)
    t = t/1000
    res.append(t)
    """
  end
  s = s * """
f = open(out_file, 'w')
f.write(\"%s\\n\" % res)
f.close()
quit()

  """
  return s
end

function benchmark_code_magma(L::Vector{Hecke.ZLat}, repetitions::Int = 10; out_file::String = "magma_timings")
  default_repetitions = repetitions
  s = """
  res := [];
  """
  for i in 1:length(L)
    if rank(L[i]) <= 15
      repititions = 100
    else
      repititions = default_repetitions
    end
    G = gram_matrix(L[i])
    s = s * to_magma(G, "G$i") * "\n"
    s = s * """
    t0 := 0.0;
    \"Doing $(i)\";
    for i in [1..$(repetitions)] do
      L := LatticeWithGram(G$(i));
      t := Cputime();
      _ := AutomorphismGroup(L);
      t0 := t0 + Cputime(t);
    end for;
    t0av := RealField()!t0/$(repetitions);
    Append(~res, [t0av, RealField()!t0]);
    """
  end
  s = s * """
  PrintFile(\"$(out_file)\", res);
  quit;
  """
  return s
end

function benchmark_magma(L::Vector{Hecke.ZLat})
  t, fp = mktemp()
  @show t
  out_file, fp_out = mktemp()
  close(fp_out)
  mcode = benchmark_code_magma(L, out_file = out_file)
  print(fp, mcode)
  close(fp)

  run(`magma $(t)`)
  rm(t)

  s = read(out_file, String)

  rm(out_file)

  t = Meta.eval(Meta.parse(s))
  return Float64[x[1] for x in t]
end

function benchmark_sage(L::Vector{Hecke.ZLat})
  t, fp = mktemp()
  @show t
  out_file, fp_out = mktemp()
  close(fp_out)
  scode = benchmark_code_sage(L, out_file = out_file)
  print(fp, scode)
  close(fp)

  run(`/home/thofmann/SageMath/sage $(t)`)
  rm(t)

  s = read(out_file, String)

  rm(out_file)

  return Meta.eval(Meta.parse(s))
end

function benchmark_hecke(L::Vector{Hecke.ZLat})
  res = Float64[]
  for i in 1:length(L)
    LL = L[i]
    println("Doing $i")
    t = @belapsed Hecke.assert_has_automorphisms($LL, redo = true) samples=10
    push!(res, t)
  end
  return res
end

function benchmark_all(L::Vector{Hecke.ZLat})
  bM = benchmark_magma(L)
  bH = benchmark_hecke(L)
  bS = benchmark_sage(L)
end