text:
# Functions in this file generate rel code with integrity constraints
#   to implement ontologies that are modeled in ORM.
#
# At a high level, we treat certain diagrams in the ORM model as
#   designating central concepts and "constellations" of relations,
#   where the relations are inferred from names and uniqueness
#   constraints that pertain to roles of ORM fact types.
#

import Pkg

# Pkg.activate((ENV["ORM_REL_GENERATOR"]))
Pkg.activate(("."))

using EzXML

include("./orm_api.jl")
include("./rel_ast.jl")

"""
Returns the Rel type of the object that plays the
  role identified by r_id
"""
function rel_type_of_role_player(r_id, M)::RelType

   etype = entity_that_plays(r_id, M)
   RelType( etype !== nothing ?
              name(etype) :
              data_type_of_player(r_id, M) .|>
                 ( dt_id ->
                     maps_to_rel_string(dt_id, M) ? "String" :
                     maps_to_rel_int(dt_id, M) ? "Int" :
                     maps_to_rel_date(dt_id, M) ? "Date" :
                     maps_to_rel_datetime(dt_id, M) ? "DateTime" :
                     "String" ) )
end

"""
Maps an ORM name to a Rel name using our naming conventions
"""
rel_name_of(orm_name::String) = pascal_to_snake_case(orm_name)

"""
Maps a role to the the name of the object that plays that role
  using Rel naming conventions
"""
rel_name_of_role_player(role, M) =
   rel_name_of(name(role_player_by_id(id(role), M)))

"""
Given an entity type etype, its preferred identifying uc,
   and the roles of the entity and its identifier in the
   identifying fact type, generate the refmode predicate
   (one that maps the entity to its identifier)
"""
function refmode_predicate(etype, uc, e_role, i_role, M)::Predicate

   function gen_refmode_pred_name()::String
      e_name = pascal_to_snake_case(name(etype))
      mode = refmode(etype)
      if mode != ""
         e_name * "_" * mode
      else
         v_name = rel_name_of_role_player(i_role, M)
         e_name * "_" * v_name
      end
   end

   i_role_name = name(i_role)
   pred_name = i_role_name != "" ? i_role_name : gen_refmode_pred_name()

   Predicate(pred_name, uc, id(e_role), id(i_role))
end

"""
Given an entity type etype, the uc that is *not* the preferred
   identifier but that applies to the identifying fact type,
   and the roles of the entity and its identifier from that
   fact type, generate the refmode inverse predicate (one that
   maps the identifier to the entity)
"""
function refmode_inverse_predicate(e, uc, e_role, i_role, M)::Predicate

   function gen_refmode_inverse_pred_name()::String
      e_name = pascal_to_snake_case(name(e))
      mode = refmode(e)
      if mode != ""
         return mode * "_to_" * e_name
      else
         v_name = rel_name_of_role_player(i_role, M)
         v_name * "_to_" * e_name
      end
   end

   e_role_name = name(e_role)
   pred_name = e_role_name != "" ? e_role_name : gen_refmode_inverse_pred_name()

   Predicate(pred_name, uc, id(i_role), id(e_role))
end

"""
Generate relations for ORM fact type ft using system generated
   relation names
"""
function system_named_relations_for(ft, M)::Vector{Predicate}
   rel_seq = Relation[]
   uc_ids = extract_uc_ids_of(ft)
   if size(uc_ids)[1] > 1
      println("[Warning]: Cannot generate relations for fact type '" *
              fact_name(ft) *
              "' because there is more than one UC and none of the roles are named.")
   else
      uc_id = uc_ids[1]
      uc = uc_by_id(uc_id, M)
      value_role_id_seq = roles_excluded_by_uc_in_fact_type(uc, ft, M)
      key_role_id_seq = extract_role_ids_from_simple(uc)
      push!(rel_seq,
            Predicate(pascal_to_snake_case(fact_name(ft)),
            uc_id,
            key_role_id_seq,
            value_role_id_seq))
   end
   return rel_seq
end

"""
Generate relations for ORM fact type ft using user-provided
   relation names. User preferences are provided by naming
   one or more of the roles of ft.
"""
function user_named_relations_for(ft, M)::Vector{Predicate}
   rel_seq = Relation[]
   # Look for relations designated by role names
   for role in extract_roles(ft)
      role_name = name(role)
      if length(role_name) == 0 continue
      else
         key_uc = uc_that_excludes_role_in_fact_type(role, ft, M)
         if key_uc === nothing
            println("[Warning]: Cannot generate relation for role with name " *
                    role_name *
                    " in fact type " * fact_name(ft) *
                    " because model lacks a non-spanning UC that excludes this role.")
         else
            roles = extract_role_ids_from_simple(key_uc)
            rel = Predicate(role_name, id(key_uc), roles, id(role))
            push!(rel_seq, rel)
         end
      end
   end
   rel_seq
end

"""
Given a Rel predicate rel and a model M, construct a
   Rel type constraint for rel.
"""
function rel_type_constraint(rel::Predicate, M)::IC
   constraint_name = "$(rel.relname)_types"
   rseq = role_seq(rel)
   num_roles = size(rseq, 1)
   vars = fresh_variables(num_roles)
   conjuncts = zip( rseq .|> (r_id -> rel_type_of_role_player(r_id, M)),
                    vars ) .|> (t -> RelAtom(t[1], t[2]))

   IC(constraint_name, vars, Implication(RelAtom(rel, vars), Conjunction(conjuncts)))
end

"""
Given an ORM mandatory constraint id mc_id, a Rel predicate rel,
   and a model M, construct a Rel mandatory constraint.
"""
function rel_mandatory_constraint(mc_id, rel, M)::IC
   var = "v"
   mc = mc_by_id(mc_id, M)
   role_ids = extract_role_ids_from_simple(mc)
   role_id = first_occurrence_from(role_ids, rel)
   atom = bind_role_to(rel, role_id, var)

   role_player = rel_type_of_role_player(role_id, M)
   constraint_name = "$(lowercase(role_player.relname))_mandatory_in_$(rel.relname)"
   IC(constraint_name,
      [var],
      Implication(RelAtom(role_player, var), atom))
end

"""
Given an ORM uniqueness constraint id uc_id, a Rel predicate rel,
   and a model M, construct a Rel many-to-one constraint.
"""
function rel_many_one_constraint(uc_id, rel, M)::Clause
   # If the uc refers to only key roles of rel, then we can
   #   implement the constraint by declaring rel to be functional
   if rel.uc_id == uc_id
      constraint_name = "$(rel.relname)_many_one"
      IC(constraint_name, IsFunctional(rel))
   else
      Unsupported("TBD: Many to one constraint for non-generating uniqueness constraint of relation " * rel.relname)
   end
end

"""
For each relation that is determined by some fact type,
   and for each ORM constraint that pertains to any of
   the roles of that relation, generate a Rel integrity
   constraint (IC)
"""
function rel_constraints_for(ft, M)::Vector{Clause}

   relations_for_non_refmode(ft, M) =
      has_named_roles(ft) ?
         user_named_relations_for(ft, M) :
         system_named_relations_for(ft, M)

   relations = Predicate[]
   constraints = Clause[]

   if is_refmode_fact_type(ft, M)
      e_role_id = i_role_id = "<undefined-role-id>"
      e_uc = i_uc = "<undefined-uc>"

      for uc_id in extract_uc_ids_of(ft)
         r_id = lone_role_of(uc_id, M)
         if entity_identified_by(uc_id, M) === nothing
            e_uc = uc_id
            e_role_id = r_id
         else
            i_uc = uc_id
            i_role_id = r_id
         end
      end

      e = entity_that_plays(e_role_id, M)
      e_role = extract_role_by_id(e_role_id, ft)
      i_role = extract_role_by_id(i_role_id, ft)

      refmode_pred = refmode_predicate(e, e_uc, e_role, i_role, M)
      refmode_inverse_pred = refmode_inverse_predicate(e, i_uc, e_role, i_role, M)

      push!(relations, refmode_pred)
      push!(constraints, rel_many_one_constraint(i_uc, refmode_inverse_pred, M))
      push!(relations, refmode_inverse_pred)
      push!(constraints, rel_many_one_constraint(e_uc, refmode_pred, M))
   else
      for rel in relations_for_non_refmode(ft, M)
         push!(relations, rel)
         if size(rel.value_role_seq)[1] > 0
            # Then there must be at least one non-spanning UC
            for uc_id in extract_uc_ids_of(ft)
               push!(constraints, rel_many_one_constraint(uc_id, rel, M))
            end
         end
      end
   end
   for rel in relations
      push!(constraints, rel_type_constraint(rel, M))
      for mc_id in extract_mc_ids_of(ft)
         push!(constraints, rel_mandatory_constraint(mc_id, rel, M))
      end
   end
   constraints
end

"""
Generates a Rel type constraint for the top-level concept
"""
function rel_constraint_for_top_level_entity(concept)::Clause
   var = "v"
   etype_name = name(concept)
   constraint_name = "$(rel_name_of(etype_name))_is_entity"
   rel_entity_base_atom = RelAtom(RelType("Entity"), var)

   IC(constraint_name,
      var,
      Implication(RelAtom(RelType(etype_name), var),
                  rel_entity_base_atom))
end

"""
Generates a Rel type constraint for concept, when it is known
   to be a subtype of some supertype
"""
function rel_constraint_for_entity_subtype(concept, supertype)::Clause
   var = "v"
   sub_name = name(concept)
   super_name = name(supertype)
   constraint_name = "$(rel_name_of(sub_name))_is_subtype_of_$(rel_name_of(super_name))"

   IC(constraint_name,
      var,
      Implication(RelAtom(RelType(sub_name), var),
                  RelAtom(RelType(super_name), var)))
end

"""
Generates and adds to 'concepts' all Rel type constraints for
   concept and its subtypes that are drawn on diagram D,
"""
function rel_constraints_for_entity_types(concept, concepts, D, M)
   if is_drawn_on(concept, D)
      if is_top_level_concept(concept, M)
         push!( concepts, rel_constraint_for_top_level_entity(concept) )
      else
         for super_type in super_types_of(concept, M)
            push!( concepts, rel_constraint_for_entity_subtype(concept, super_type) )
         end
      end
   end
   for sub_concept in sub_types_of(concept, M)
      rel_constraints_for_entity_types(sub_concept, concepts, D, M)
   end
end

function generate(orm_model::String, output_folder::String)
   println("Generating schemas for: " * orm_model * " in " * output_folder)
   !isdir(output_folder) && mkpath(output_folder)

   doc = root(readxml(orm_model))
   M = findfirst("orm:ORMModel", doc)

   for D in findall("ormDiagram:ORMDiagram", doc)
      concept = central_concept(D, M)
      concept === nothing && continue

      filename = lowercase(name(concept)) * "_schema.rel"
      open("$output_folder/$filename", "w") do io
         println(io)
         println(io, "// Constraints for the '" * name(concept) * "' concept and its subtypes")
         concept_constraints = Clause[]
         rel_constraints_for_entity_types(concept, concept_constraints, D, M)
         for constraint in concept_constraints
            println(io, emit(constraint))
         end
         fact_types = fact_types_in_diagram(D, M)
         if refmode(concept) != ""
            # Add the fact type that identifies the central concept,
            #   even though it is not displayed on the diagram
            ft = identifying_fact_type_for(concept, M)
            push!(fact_types, ft)
         end
         for ft in fact_types
            println(io)
            println(io, "// Constraints for relations modeled by the '" * name(concept) * "' fact type")
            for ft_constraint in rel_constraints_for(ft, M)
               println(io, emit(ft_constraint))
            end
         end
      end
   end
end

max_stars_repo_name: RelationalAI-oss/kse-tools

max_issues_repo_name: RelationalAI-oss/kse-tools

max_forks_repo_name: RelationalAI-oss/kse-tools


text:
export compute_relative_feasibility

"""
Compute transform-domain relative feasibility w.r.t. a set
r_feas =  ||P(A*x)-A*x||/||A*x||
Helper function that is intended for the parallel version only, where just a single
operator is sent to a worker where this script runs
"""
function compute_relative_feasibility(m::Vector{TF}, feasibility_initial::Vector{TF},
                                      TD_OP::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}},
                                       P_sub) where {TF<:Real,TI<:Integer}
# function compute_relative_feasibility(m::Vector{TF},feasibility_initial::TF,
#     TD_OP::Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}},
#      P_sub) where {TF<:Real,TI<:Integer}

    #note that P_sub mutates the input in-place
    s_temp = TD_OP[1]*m
    norm(P_sub[1](deepcopy(s_temp)) .- s_temp)
     # feasibility_initial[1]=norm(P_sub[1](pm) .- pm) ./ (norm(pm)+(100*eps(TF)))
    feasibility_initial[1] = norm(P_sub[1](deepcopy(s_temp)) .- s_temp) ./ (norm(s_temp)+(100*eps(TF)))
    s_temp = []

    #return feasibility_initial
end

ext: jl

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# This file is a modified from Julia's LinearAlgebra/qr.jl. License is MIT: https://julialang.org/license

# QL and Hessenberg Factorizations
"""
    QL <: Factorization

A QL matrix factorization stored in a packed format, typically obtained from
[`ql`](@ref). If ``A`` is an `m`×`n` matrix, then

```math
A = Q L
```

where ``Q`` is an orthogonal/unitary matrix and ``L`` is lower triangular.
The matrix ``Q`` is stored as a sequence of Householder reflectors ``v_i``
and coefficients ``\\tau_i`` where:

```math
Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T).
```

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

The object has two fields:

* `factors` is an `m`×`n` matrix.

  - The lower triangular part contains the elements of ``L``, that is `L =
    tril(F.factors)` for a `QL` object `F`.

  - The subdiagonal part contains the reflectors ``v_i`` stored in a packed format where
    ``v_i`` is the ``i``th column of the matrix `V = I + tril(F.factors, -1)`.

* `τ` is a vector  of length `min(m,n)` containing the coefficients ``\tau_i``.

"""
struct QL{T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} <: Factorization{T}
    factors::S
    τ::Tau

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

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

function generic_qlfactUnblocked!(A::AbstractMatrix{T}) where {T}
    require_one_based_indexing(A)
    m, n = size(A)
    τ = zeros(T, min(m,n))
    for k = n:-1:max((n - m + 1 + (T<:Real)),1)
        μ = m+k-n
        x = view(A, μ:-1:1, k)
        τk = reflector!(x)
        τ[k-n+min(m,n)] = τk
        reflectorApply!(x, τk, view(A, μ:-1:1, 1:k-1))
    end
    QL(A, τ)
end

@inline qlfactUnblocked!(A::AbstractMatrix) = generic_qlfactUnblocked!(A)
@inline qlfactUnblocked!(A::StridedMatrix{T}) where T<:BlasFloat = QL(LAPACK.geqlf!(A)...)

# Generic fallbacks

"""
    ql!(A, pivot=Val(false))

`ql!` is the same as [`ql`](@ref) when `A` is a subtype of
`StridedMatrix`, but saves space by overwriting the input `A`, instead of creating a copy.
An [`InexactError`](@ref) exception is thrown if the factorization produces a number not
representable by the element type of `A`, e.g. for integer types.

# Examples
```jldoctest
julia> a = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> ql!(a)
LinearAlgebra.QLCompactWY{Float64,Array{Float64,2}}
Q factor:
2×2 LinearAlgebra.QLCompactWYQ{Float64,Array{Float64,2}}:
 -0.316228  -0.948683
 -0.948683   0.316228
 L factor:
2×2 Array{Float64,2}:
 -3.16228  -4.42719
  0.0      -0.632456

julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> ql!(a)
ERROR: InexactError: Int64(-3.1622776601683795)
Stacktrace:
[...]
```
"""
@inline _ql!(layout, axes, A, ::Val{false}) = qlfactUnblocked!(A)
@inline _ql!(layout, axes, A) = ql!(A, Val(false))
@inline ql!(A::AbstractMatrix; kwds...) = _ql!(MemoryLayout(typeof(A)), axes(A), A; kwds...)
@inline ql!(A::AbstractMatrix, val; kwds...) = _ql!(MemoryLayout(typeof(A)), axes(A), A, val; kwds...)

@inline _qleltype(::Type{T}) where T = typeof(zero(T)/sqrt(abs2(one(T))))

"""
    ql(A, pivot=Val(false)) -> F

Compute the QL factorization of the matrix `A`: an orthogonal (or unitary if `A` is
complex-valued) matrix `Q`, and a lower triangular matrix `L` such that

```math
A = Q L
```

The returned object `F` stores the factorization in a packed format:

 - `F` is a [`QL`](@ref) object.

The individual components of the decomposition `F` can be retrieved via property accessors:

 - `F.Q`: the orthogonal/unitary matrix `Q`
 - `F.L`: the lower triangular matrix `L`

Iterating the decomposition produces the components `Q`, `L`, and if extant `p`.

The following functions are available for the `QL` objects: [`inv`](@ref), [`size`](@ref),
and [`\\`](@ref). When `A` is rectangular, `\\` will return a least squares
solution and if the solution is not unique, the one with smallest norm is returned. When
`A` is not full rank, factorization with (column) pivoting is required to obtain a minimum
norm solution.

Multiplication with respect to either full/square or non-full/square `Q` is allowed, i.e. both `F.Q*F.L`
and `F.Q*A` are supported. A `Q` matrix can be converted into a regular matrix with
[`Matrix`](@ref).  This operation returns the "thin" Q factor, i.e., if `A` is `m`×`n` with `m>=n`, then
`Matrix(F.Q)` yields an `m`×`n` matrix with orthonormal columns.  To retrieve the "full" Q factor, an
`m`×`m` orthogonal matrix, use `F.Q*Matrix(I,m,m)`.  If `m<=n`, then `Matrix(F.Q)` yields an `m`×`m`
orthogonal matrix.

# Examples
```jldoctest
julia> A = [3.0 -6.0; 4.0 -8.0; 0.0 1.0]
3×2 Array{Float64,2}:
 3.0  -6.0
 4.0  -8.0
 0.0   1.0

julia> F = ql(A)
LinearAlgebra.QLCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QLCompactWYQ{Float64,Array{Float64,2}}:
 -0.6   0.0   0.8
 -0.8   0.0  -0.6
  0.0  -1.0   0.0
L factor:
2×2 Array{Float64,2}:
 -5.0  10.0
  0.0  -1.0

julia> F.Q * F.L == A
true
```
"""
@inline ql(A::AbstractMatrix, args...; kwds...) = _ql(MemoryLayout(typeof(A)), axes(A), A, args...; kwds...)

function _ql(layout, axes, A, args...; kwds...)
    require_one_based_indexing(A)
    AA = similar(A, _qleltype(eltype(A)), size(A))
    copyto!(AA, A)
    return ql!(AA, args...; kwds...)
end
@inline ql(x::Number, args...; kwds...) = ql(fill(x,1,1), args...; kwds...)
@inline function ql(v::AbstractVector, args...; kwds...)
    require_one_based_indexing(v)
    ql(reshape(v, (length(v), 1)), args...; kwds...)
end

# Conversions
QL{T}(A::QL) where {T} = QL(convert(AbstractMatrix{T}, A.factors), convert(AbstractVector{T}, A.τ))
Factorization{T}(A::QL{T}) where {T} = A
Factorization{T}(A::QL) where {T} = QL{T}(A)
AbstractMatrix(F::QL) = F.Q * F.L
AbstractArray(F::QL) = AbstractMatrix(F)
Matrix(F::QL) = Array(AbstractArray(F))
Array(F::QL) = Matrix(F)

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::QL)
    summary(io, F); println(io)
    println(io, "Q factor:")
    show(io, mime, F.Q)
    println(io, "\nL factor:")
    show(io, mime, F.L)
end

@inline function getL(F::QL, _)
    m, n = size(F)
    tril!(getfield(F, :factors)[end-min(m,n)+1:end, 1:n], max(n-m,0))
end
@inline getQ(F::QL, _) = QLPackedQ(getfield(F, :factors), F.τ)

getL(F::QL, ::Tuple{AbstractVector,AbstractVector}) = getL(F, size(F))
getQ(F::QL, ::Tuple{AbstractVector,AbstractVector}) = getQ(F, size(F))

getL(F::QL) = getL(F, axes(F.factors))
getQ(F::QL) = getQ(F, axes(F.factors))

function getproperty(F::QL, d::Symbol)
    if d == :L
        return getL(F)
    elseif d == :Q
        return getQ(F)
    else
        getfield(F, d)
    end
end

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

"""
    QLPackedQ <: AbstractMatrix

The orthogonal/unitary ``Q`` matrix of a QL factorization stored in [`QL`](@ref).
"""
struct QLPackedQ{T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} <: LayoutQ{T}
    factors::S
    τ::Tau

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

QLPackedQ{T}(Q::QLPackedQ) where {T} = QLPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(AbstractVector{T}, Q.τ))
AbstractMatrix{T}(Q::QLPackedQ{T}) where {T} = Q
AbstractMatrix{T}(Q::QLPackedQ) where {T} = QLPackedQ{T}(Q)

size(Q::QLPackedQ, dim::Integer) = size(getfield(Q, :factors), dim == 2 ? 1 : dim)

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

## Multiplication by Q
function _mul(A::QLPackedQ, B::AbstractMatrix)
    TAB = promote_type(eltype(A), eltype(B))
    Anew = convert(AbstractMatrix{TAB}, A)
    if size(A.factors, 1) == size(B, 1)
        Bnew = copy_oftype(B, TAB)
    elseif size(A.factors, 2) == size(B, 1)
        Bnew = [zeros(TAB, size(A.factors, 1) - size(B,1), size(B, 2)); B]
    else
        throw(DimensionMismatch("first dimension of matrix must have size either $(size(A.factors, 1)) or $(size(A.factors, 2))"))
    end
    lmul!(Anew, Bnew)
end

(*)(A::QLPackedQ, B::StridedMatrix) = _mul(A, B)

### QB
abstract type AbstractQLLayout <: MemoryLayout end

"""
    QLPackedLayout{SLAY,TLAY}()

represents a Packed QL factorization whose
factors are stored with layout SLAY and τ stored with layout TLAY
"""
struct QLPackedLayout{SLAY,TLAY} <: AbstractQLLayout end
struct QLPackedQLayout{SLAY,TLAY} <: AbstractQLayout end
struct AdjQLPackedQLayout{SLAY,TLAY} <: AbstractQLayout end

MemoryLayout(::Type{<:QL{<:Any,Mat,Tau}}) where {Mat,Tau} =
    QLPackedLayout{typeof(MemoryLayout(Mat)),typeof(MemoryLayout(Tau))}()

adjointlayout(::Type, ::QLPackedQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQLPackedQLayout{SLAY,TLAY}()

MemoryLayout(::Type{<:QLPackedQ{<:Any,Mat,Tau}}) where {Mat,Tau} =
    QLPackedQLayout{typeof(MemoryLayout(Mat)),typeof(MemoryLayout(Tau))}()

function materialize!(M::Lmul{<:QLPackedQLayout})
    A,B = M.A, M.B
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    @inbounds begin
        for k = max(nA - mA + 1,1):nA
            μ = mA+k-nA
            for j = 1:nB
                vBj = B[μ,j]
                for i = 1:μ-1
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = A.τ[k-nA+min(mA,nA)]*vBj
                B[μ,j] -= vBj
                for i = 1:μ-1
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end

function materialize!(M::Lmul{<:AdjQLPackedQLayout})
    adjA,B = M.A, M.B
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    @inbounds begin
        for k = nA:-1:max(nA - mA + 1,1)
            μ = mA+k-nA
            for j = 1:nB
                vBj = B[μ,j]
                for i = 1:μ-1
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = conj(A.τ[k-nA+min(mA,nA)])*vBj
                B[μ,j] -= vBj
                for i = 1:μ-1
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end

### QBc/QcBc
function materialize!(M::Rmul{<:Any,<:QLPackedQLayout})
    A,Q = M.A, M.B
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    begin
        for k = nQ:-1:max(nQ - mQ + 1,1)
            μ = mQ+k-nQ
            for i = 1:mA
                vAi = A[i,μ]
                for j = 1:μ-1
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*Q.τ[k-nQ+min(mQ,nQ)]
                A[i,μ] -= vAi
                for j = 1:μ-1
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end

### AQc
function materialize!(M::Rmul{<:Any,<:AdjQLPackedQLayout})
    A,adjQ = M.A, M.B
    Q = adjQ.parent
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    @inbounds begin
        for k = max(nQ - mQ + 1,1):nQ
            μ = mQ+k-nQ
            for i = 1:mA
                vAi = A[i,μ]
                for j = 1:μ-1
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*conj(Q.τ[k-nQ+min(mQ,nQ)])
                A[i,μ] -= vAi
                for j = 1:μ-1
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end

# Julia implementation similar to xgelsy
ldiv!(F::QL, B::AbstractVecOrMat) = ArrayLayouts.ldiv!(F, B)

function materialize!(Ldv::Ldiv{<:QLPackedLayout,<:Any,<:Any,<:AbstractMatrix{T}}) where T
    A,B = Ldv.A,Ldv.B
    m, n = size(A)
    minmn = min(m,n)
    mB, nB = size(B)
    ArrayLayouts.lmul!(adjoint(A.Q), view(B, 1:m, :))
    L = A.L
    @inbounds begin
        if n > m # minimum norm solution
            τ = zeros(T,m)
            for k = m:-1:1 # Trapezoid to triangular by elementary operation
                x = view(L, k, [k; m + 1:n])
                τk = reflector!(x)
                τ[k] = conj(τk)
                for i = 1:k - 1
                    vLi = L[i,k]
                    for j = m + 1:n
                        vLi += L[i,j]*x[j - m + 1]'
                    end
                    vLi *= τk
                    L[i,k] -= vLi
                    for j = m + 1:n
                        L[i,j] -= vLi*x[j - m + 1]
                    end
                end
            end
        end
        ArrayLayouts.ldiv!(LowerTriangular(view(L, 1:minmn, :)), view(B, 1:minmn, :))
        if n > m # Apply elementary transformation to solution
            B[m + 1:mB,1:nB] .= zero(T)
            for j = 1:nB
                for k = 1:m
                    vBj = B[k,j]
                    for i = m + 1:n
                        vBj += B[i,j]*L[k,i]'
                    end
                    vBj *= τ[k]
                    B[k,j] -= vBj
                    for i = m + 1:n
                        B[i,j] -= L[k,i]*vBj
                    end
                end
            end
        end
    end
    return B
end
materialize!(Ldv::Ldiv{<:QLPackedLayout,<:Any,<:Any,<:AbstractVector{T}}) where T =
    ldiv!(Ldv.A, reshape(Ldv.B, length(Ldv.B), 1))[:]

function (\)(A::QL{TA}, B::AbstractVecOrMat{TB}) where {TA,TB}
    require_one_based_indexing(B)
    S = promote_type(TA,TB)
    m, n = size(A)
    m == size(B,1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(B,1)) rows"))

    AA = Factorization{S}(A)

    X = _zeros(S, B, n)
    X[1:size(B, 1), :] = B
    ldiv!(AA, X)
    return _cut_B(X, 1:n)
end

function (\)(A::QL{T}, BIn::VecOrMat{Complex{T}}) where T<:BlasReal
    require_one_based_indexing(BIn)
    m, n = size(A)
    m == size(BIn, 1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(BIn,1)) rows"))

# |z1|z3|  reinterpret  |x1|x2|x3|x4|  transpose  |x1|y1|  reshape  |x1|y1|x3|y3|
# |z2|z4|      ->       |y1|y2|y3|y4|     ->      |x2|y2|     ->    |x2|y2|x4|y4|
#                                                 |x3|y3|
#                                                 |x4|y4|
    B = reshape(copy(transpose(reinterpret(T, reshape(BIn, (1, length(BIn)))))), size(BIn, 1), 2*size(BIn, 2))

    X = _zeros(T, B, n)
    X[1:size(B, 1), :] = B

    ldiv!(A, X)

# |z1|z3|  reinterpret  |x1|x2|x3|x4|  transpose  |x1|y1|  reshape  |x1|y1|x3|y3|
# |z2|z4|      <-       |y1|y2|y3|y4|     <-      |x2|y2|     <-    |x2|y2|x4|y4|
#                                                 |x3|y3|
#                                                 |x4|y4|
    XX = reshape(collect(reinterpret(Complex{T}, copy(transpose(reshape(X, div(length(X), 2), 2))))), _ret_size(A, BIn))
    return _cut_B(XX, 1:n)
end


text:
function equivalent_to_wkt(geom::GeoInterface.AbstractGeometry, wkt::String)
    test_geom = readgeom(wkt)
    @test writegeom(geom) == writegeom(test_geom)
end

function factcheck(f::Function, geom::String, expected::String)
    result = f(readgeom(geom))
    equivalent_to_wkt(result, expected)
end

function factcheck(f::Function, geom::String, expected::Bool)
    @test f(readgeom(geom)) == expected
end

function factcheck(f::Function, g1::String, g2::String, expected::String)
    result = f(readgeom(g1),readgeom(g2))
    equivalent_to_wkt(result, expected)
end

function factcheck(f::Function, g1::String, g2::String, expected::Bool)
    @test f(readgeom(g1),readgeom(g2)) == expected
    @test f(prepareGeom(readgeom(g1)),readgeom(g2)) == expected
end

@testset "GEOS operations" begin
    ls = LineString(Vector{Float64}[[8,1],[9,1],[9,2],[8,2]])
    pt = interpolate(ls, 2.5)
    @test GeoInterface.coordinates(pt) ≈ [8.5, 2.0] atol=1e-5
    for (pt,dist,dest) in [(Point(10.0,1.0), 1.0, Point(9.0,1.0)),
                           (Point( 9.0,1.0), 1.0, Point(9.0,1.0)),
                           (Point(10.0,0.0), 1.0, Point(9.0,1.0)),
                           (Point( 9.0,2.0), 2.0, Point(9.0,2.0)),
                           (Point( 8.7,1.5), 1.5, Point(9.0,1.5))]
        test_dist = project(ls, pt)
        @test test_dist ≈ dist atol=1e-5
        @test equals(interpolate(ls, test_dist), dest)
    end

    # GEOSConvexHullTest
    input = readgeom("MULTIPOINT (130 240, 130 240, 130 240, 570 240, 570 240, 570 240, 650 240)")
    expected = readgeom("LINESTRING (130 240, 650 240)")
    output = convexhull(input)
    @test !isEmpty(output)
    @test writegeom(output) == writegeom(expected)

    # LibGEOS.delaunayTriangulationTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = delaunayTriangulationEdges(g1)
    @test isEmpty(g1)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :MultiLineString

    g1 = readgeom("POINT(0 0)")
    g2 = delaunayTriangulation(g1)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :GeometryCollection

    g1 = readgeom("MULTIPOINT(0 0, 5 0, 10 0)")
    g2 = delaunayTriangulation(g1, 0.0)
    @test isEmpty(g2)
    @test GeoInterface.geotype(g2) == :GeometryCollection
    g2 = delaunayTriangulationEdges(g1, 0.0)
    equivalent_to_wkt(g2, "MULTILINESTRING ((5 0, 10 0), (0 0, 5 0))")

    g1 = readgeom("MULTIPOINT(0 0, 10 0, 10 10, 11 10)")
    g2 = delaunayTriangulationEdges(g1, 2.0)
    equivalent_to_wkt(g2, "MULTILINESTRING ((0 0, 10 10), (0 0, 10 0), (10 0, 10 10))")

    # GEOSDistanceTest
    g1 = readgeom("POINT(10 10)")
    g2 = readgeom("POINT(3 6)")
    @test distance(g1, g2) ≈ 8.06225774829855 atol=1e-12

    # GEOSGeom_extractUniquePointsTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = uniquePoints(g1)
    @test isEmpty(g2)

    g1 = readgeom("MULTIPOINT(0 0, 0 0, 1 1)")
    g2 = uniquePoints(g1)
    @test equals(g2, readgeom("MULTIPOINT(0 0, 1 1)"))

    g1 = readgeom("GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))")
    @test LibGEOS.equals(uniquePoints(g1), readgeom("MULTIPOINT(0 0, 1 1, 2 2, 5 5, 0 2)"))

    # GEOSGetCentroidTest
    test_centroid(geom::String, expected::String) = factcheck(centroid, geom, expected)
    test_centroid("POINT(10 0)", "POINT (10 0)")
    test_centroid("LINESTRING(0 0, 10 0)", "POINT (5 0)")
    test_centroid("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))", "POINT (5 5)")
    test_centroid("LINESTRING EMPTY", "POINT EMPTY")

    # GEOSIntersectionTest
    test_intersection(g1::String, g2::String, expected::String) = factcheck(intersection, g1, g2, expected)
    test_intersection("POLYGON EMPTY", "POLYGON EMPTY", "GEOMETRYCOLLECTION EMPTY")
    test_intersection("POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", "POINT(2 2)")
    test_intersection("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", "POLYGON((-1 1,-1 2,2 2,2 1,-1 1))", "POLYGON ((0 1, 0 2, 2 2, 2 1, 0 1))")
    test_intersection("MULTIPOLYGON(((0 0,5 10,10 0,0 0),(1 1,1 2,2 2,2 1,1 1),(100 100,100 102,102 102,102 100,100 100)))",
                      "POLYGON((0 1,0 2,10 2,10 1,0 1))",
                      "GEOMETRYCOLLECTION (LINESTRING (1 2, 2 2), LINESTRING (2 1, 1 1), POLYGON ((0.5 1, 1 2, 1 1, 0.5 1)), POLYGON ((9 2, 9.5 1, 2 1, 2 2, 9 2)))")

    # LineString_PointTest
    g1 = readgeom("LINESTRING(0 0, 5 5, 10 10)")
    @test !isClosed(g1)
    @test GeoInterface.geotype(g1) == :LineString
    @test numPoints(g1) == 3
    @test geomLength(g1) ≈ sqrt(100 + 100) atol=1e-5
    @test GeoInterface.coordinates(startPoint(g1)) ≈ [0,0] atol=1e-5
    @test GeoInterface.coordinates(endPoint(g1)) ≈ [10,10] atol=1e-5

    # GEOSNearestPointsTest
    g1 = readgeom("POLYGON EMPTY")
    g2 = readgeom("POLYGON EMPTY")
    @test length(nearestPoints(g1, g2)) == 0

    g1 = readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
    g2 = readgeom("POLYGON((8 8, 9 9, 9 10, 8 8))")
    points = nearestPoints(g1, g2)
    @test length(points) == 2
    @test GeoInterface.coordinates(points[1])[1:2] == [5.0,5.0]
    @test GeoInterface.coordinates(points[2])[1:2] == [8.0,8.0]

    # GEOSNodeTest
    g1 = node(readgeom("LINESTRING(0 0, 10 10, 10 0, 0 10)"))
    LibGEOS.normalize!(g1)
    equivalent_to_wkt(g1, "MULTILINESTRING ((5 5, 10 0, 10 10, 5 5), (0 10, 5 5), (0 0, 5 5))")

    g1 = node(readgeom("MULTILINESTRING((0 0, 2 0, 4 0),(5 0, 3 0, 1 0))"))
    LibGEOS.normalize!(g1)
    equivalent_to_wkt(g1, "MULTILINESTRING ((4 0, 5 0), (3 0, 4 0), (2 0, 3 0), (1 0, 2 0), (0 0, 1 0))")

    g1 = node(readgeom("MULTILINESTRING((0 0, 2 0, 4 0),(0 0, 2 0, 4 0))"))
    LibGEOS.normalize!(g1)
    equivalent_to_wkt(g1, "MULTILINESTRING ((2 0, 4 0), (0 0, 2 0))")

    # GEOSPointOnSurfaceTest
    test_pointonsurface(geom::String, expected::String) = factcheck(pointOnSurface, geom, expected)
    test_pointonsurface("POINT(10 0)", "POINT (10 0)")
    test_pointonsurface("LINESTRING(0 0, 5 0, 10 0)", "POINT (5 0)")
    test_pointonsurface("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))", "POINT (5 5)")
    test_pointonsurface("LINESTRING EMPTY", "POINT EMPTY")
    test_pointonsurface("LINESTRING(0 0, 0 0)", "POINT (0 0)")

    g1 = readgeom(
    """POLYGON((
    56.528666666700 25.2101666667,
    56.529000000000 25.2105000000,
    56.528833333300 25.2103333333,
    56.528666666700 25.2101666667))""")
    @test GeoInterface.coordinates(pointOnSurface(g1)) ≈ [56.528917,25.210417] atol=1e-5

    # GEOSSharedPathsTest
    factcheck(sharedPaths,
              "LINESTRING (-30 -20, 50 60, 50 70, 50 0)",
              "LINESTRING (-29 -20, 50 60, 50 70, 51 0)",
              "GEOMETRYCOLLECTION (MULTILINESTRING ((50 60, 50 70)), MULTILINESTRING EMPTY)")

    # GEOSSimplifyTest
    g1 = readgeom("POLYGON EMPTY")
    @test isEmpty(g1)
    g2 = simplify(g1, 43.2)
    @test isEmpty(g2)
    g1 = readgeom(
    """POLYGON((
    56.528666666700 25.2101666667,
    56.529000000000 25.2105000000,
    56.528833333300 25.2103333333,
    56.528666666700 25.2101666667))""")
    equivalent_to_wkt(simplify(g1, 0.0), "POLYGON EMPTY")
    @test equals(g1, topologyPreserveSimplify(g1, 43.2))

    # GEOSSnapTest
    function test_snap(g1::String, g2::String, expected::String, tol::Float64=0.0)
        equivalent_to_wkt(snap(readgeom(g1), readgeom(g2), tol), expected)
    end
    test_snap("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))", "POINT(0.5 0)",
              "POLYGON ((0.5 0, 10 0, 10 10, 0 10, 0.5 0))", 1.0)
    test_snap("LINESTRING (-30 -20, 50 60, 50 0)", "LINESTRING (-29 -20, 40 60, 51 0)",
              "LINESTRING (-29 -20, 50 60, 51 0)", 2.0)
    test_snap("LINESTRING (-20 -20, 50 50, 100 100)", "LINESTRING (-10 -9, 40 20, 80 79)",
              "LINESTRING (-20 -20, -10 -9, 50 50, 80 79, 100 100)", 2.0)
    test_snap("LINESTRING(0 0, 10 0)", "LINESTRING(0 0, 9 0)",
              "LINESTRING(0 0, 9 0)", 2.0)
    test_snap("LINESTRING(0 0, 10 0)", "LINESTRING(0 0, 9 0, 10 0, 11 0)",
              "LINESTRING(0 0, 9 0, 10 0, 11 0)", 2.0)
    test_snap("LINESTRING(0 3,4 1,0 1)", "MULTIPOINT(5 0,4 1)",
              "LINESTRING (0 3, 4 1, 5 0, 0 1)", 2.0)
    test_snap("LINESTRING(0 3,4 1,0 1)", "MULTIPOINT(4 1,5 0)",
              "LINESTRING (0 3, 4 1, 5 0, 0 1)", 2.0)
    test_snap("LINESTRING(0 0,10 0,10 10,0 10,0 0)", "MULTIPOINT(0 0,-1 0)",
              "LINESTRING (-1 0, 0 0, 10 0, 10 10, 0 10, -1 0)", 3.0)
    test_snap("LINESTRING(0 2,5 2,9 2,5 0)", "POINT(5 0)",
              "LINESTRING (0 2, 5 2, 9 2, 5 0)", 3.0)
    test_snap("LINESTRING(-71.1317 42.2511,-71.1317 42.2509)", "MULTIPOINT(-71.1261 42.2703,-71.1257 42.2703,-71.1261 42.2702)",
              "LINESTRING (-71.1257 42.2703, -71.1261 42.2703, -71.1261 42.2702, -71.1317 42.2509)", .5)

    # GEOSUnaryUnionTest
    test_unaryunion(geom::String, expected::String) = factcheck(unaryUnion, geom, expected)
    test_unaryunion("POINT EMPTY", "GEOMETRYCOLLECTION EMPTY")
    test_unaryunion("POINT (6 3)", "POINT (6 3)")
    test_unaryunion("POINT (4 5 6)", "POINT Z (4 5 6)")
    test_unaryunion("MULTIPOINT (4 5, 6 7, 4 5, 6 5, 6 7)", "MULTIPOINT (4 5, 6 5, 6 7)")
    test_unaryunion("GEOMETRYCOLLECTION (POINT(4 5), MULTIPOINT(6 7, 6 5, 6 7), LINESTRING(0 5, 10 5), LINESTRING(4 -10, 4 10))",
                    "GEOMETRYCOLLECTION (POINT (6 7), LINESTRING (4 -10, 4 5), LINESTRING (4 5, 4 10), LINESTRING (0 5, 4 5), LINESTRING (4 5, 10 5))")
    test_unaryunion("GEOMETRYCOLLECTION (POINT(4 5), MULTIPOINT(6 7, 6 5, 6 7), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)))",
                    "GEOMETRYCOLLECTION (POINT (6 7), POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (5 6, 7 6, 7 8, 5 8, 5 6)))")
    test_unaryunion("GEOMETRYCOLLECTION (MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5)), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)))",
                    "GEOMETRYCOLLECTION (LINESTRING (5 7, 7 7), LINESTRING (10 7, 12 7), LINESTRING (5.5 7.5, 6.5 7.5), POLYGON ((10 7, 10 0, 0 0, 0 10, 10 10, 10 7), (5 6, 7 6, 7 7, 7 8, 5 8, 5 7, 5 6)))")
    test_unaryunion("GEOMETRYCOLLECTION (MULTILINESTRING((5 7, 12 7), (4 5, 6 5), (5.5 7.5, 6.5 7.5)), POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 6, 7 6, 7 8, 5 8, 5 6)), MULTIPOINT(6 6.5, 6 1, 12 2, 6 1))",
                    "GEOMETRYCOLLECTION (POINT (6 6.5), POINT (12 2), LINESTRING (5 7, 7 7), LINESTRING (10 7, 12 7), LINESTRING (5.5 7.5, 6.5 7.5), POLYGON ((10 7, 10 0, 0 0, 0 10, 10 10, 10 7), (5 6, 7 6, 7 7, 7 8, 5 8, 5 7, 5 6)))")

    for (g1, g2, testvalue) in (
            ("POLYGON EMPTY", "POLYGON EMPTY", false),
            ("POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", false),
            ("POINT(2 2)", "POLYGON((1 1,1 5,5 5,5 1,1 1))", true),
            ("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", "POLYGON((1 1,1 2,2 2,2 1,1 1))", false),
            ("POLYGON((1 1,1 2,2 2,2 1,1 1))", "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", true),
            ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 201, 300 300)", false),
            ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 200, 300 300)", false),
            ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 300 100)", false)
        )
        for f in (within, coveredby)
            factcheck(f, g1, g2, testvalue)
        end
        for f in (contains, covers)
            factcheck(f, g2, g1, testvalue)
        end
    end

    for (g1, g2, testvalue) in (
            ("POLYGON EMPTY", "POLYGON EMPTY", false),
            ("POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", true),
            ("POINT(2 2)", "POLYGON((1 1,1 5,5 5,5 1,1 1))", true),
            ("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", "POLYGON((1 1,1 2,2 2,2 1,1 1))", true),
            ("POLYGON((1 1,1 2,2 2,2 1,1 1))", "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", true),
            ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 201, 300 300)", false),
            ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 200 200, 300 300)", true),
            ("LINESTRING (100 200, 300 200)", "LINESTRING (100 300, 300 100)", true)
        )
        factcheck(intersects, g1, g2, testvalue)
        factcheck(disjoint, g1, g2, !testvalue)
    end

    # GEOSisClosedTest
    @test !isClosed(readgeom("LINESTRING(0 0, 1 0, 1 1)"))
    @test isClosed(readgeom("LINESTRING(0 0, 0 1, 1 1, 0 0)"))

    # Buffer should return Polygon or MultiPolygon
    @test buffer(MultiPoint([[1.0, 1.0], [2.0, 2.0], [2.0, 0.0]]), 0.1) isa LibGEOS.MultiPolygon
    @test buffer(MultiPoint([[1.0, 1.0], [2.0, 2.0], [2.0, 0.0]]), 10) isa LibGEOS.Polygon

    # bufferWithStyle
    g1 = bufferWithStyle(readgeom("LINESTRING(0 0,0 1,1 1)"), 0.1, endCapStyle=LibGEOS.GEOSBUF_CAP_FLAT, joinStyle=LibGEOS.GEOSBUF_JOIN_BEVEL)
    g2 = readgeom("POLYGON((-0.1 0.0,-0.1 1.0,0.0 1.1,1.0 1.1,1.0 0.9,0.1 0.9,0.1 0.0,-0.1 0.0))")
    @test equals(g1, g2)

    g1 = bufferWithStyle(readgeom("LINESTRING(0 0,0 1,1 1)"), 0.1, endCapStyle=LibGEOS.GEOSBUF_CAP_SQUARE, joinStyle=LibGEOS.GEOSBUF_JOIN_MITRE)
    g2 = readgeom("POLYGON((-0.1 -0.1,-0.1 1.1,1.1 1.1,1.1 0.9,0.1 0.9,0.1 -0.1,-0.1 -0.1))")
    @test equals(g1, g2)

    g1 = bufferWithStyle(readgeom("POLYGON((-1 -1,1 -1,1 1,-1 1,-1 -1))"), 0.2, joinStyle=LibGEOS.GEOSBUF_JOIN_MITRE)
    g2 = readgeom("POLYGON((-1.2 1.2,1.2 1.2,1.2 -1.2,-1.2 -1.2,-1.2 1.2))")
    @test equals(g1, g2)

    #GEOSWKBReadWriteTest
    x, y = 2.0, 3.0
    srid = 4326
    wkt = "Point($x $y)"
    hexwkb = "010100000000000000000000400000000000000840"
    bytewkb = hex2bytes(hexwkb) |> String
    hexwkb_srid = "0101000020E610000000000000000000400000000000000840"
    bytewkb_srid = hex2bytes(hexwkb_srid) |> String

    pt = readgeom(wkt)
    @test writewkb(pt) == bytewkb
    @test writewkb(pt, hex=true) == hexwkb
    @test writewkb(pt, 4326) == bytewkb_srid
    @test writewkb(pt, 4326, hex=true) == hexwkb_srid

    pt = readwkb(bytewkb)
    @test isa(pt, GeoInterface.AbstractPoint)
    @test GeoInterface.coordinates(pt) == [x, y]

    pt = readwkb(hexwkb, hex=true)
    @test GeoInterface.geotype(pt) == :Point
    @test GeoInterface.coordinates(pt) == [x, y]

    pt = readwkb(bytewkb_srid)
    @test GeoInterface.geotype(pt) == :Point
    @test GeoInterface.coordinates(pt) == [x, y]
    @test getSRID(pt) == srid

    pt = readwkb(hexwkb_srid, hex=true)
    @test GeoInterface.geotype(pt) == :Point
    @test GeoInterface.coordinates(pt) == [x, y]
    @test getSRID(pt) == srid

    # GEOSSRIDGetSetTest
    pt = readgeom(wkt)
    @test_broken getSRID(pt)
    setSRID(pt, srid)
    @test getSRID(pt) == srid
end


text:
#=
Philosophy and convention:
A node is defined as the center of a grid cell. Grid cells may be those
in a primary grid or those in a dual grid.

The definitions "dual" and "primal" are made based on their conventional use
in holding fluid dynamic data. For example, pressure is held in primal nodes,
vorticity and streamfunction at dual nodes.

The definitions here are made with a "dual grid" in mind. That is, the grid is
defined by an integer number of dual cells in each direction. If a "primal grid"
is needed, then all of the defintions can be swapped (primal -> dual, dual -> primal).

Also, note that there might be dual cells that are "ghosts" (i.e. lie outside
the grid), but these are not distinguished in these basic definitions and operators.
=#

#=
Notes on data types:
All data types can be reduced to one of two different types in each direction.
(1) Lying on dual cell centers (C), numbered 1 through N
(2) Lying on dual cell edges (E), numbered 1 through N-1 (i.e. lying midway between C)

For example, Nodes{Dual} by definition lie on dual cell centers in both directions,
so they would be C x C. Nodes{Primal} are aligned with the corners of the dual
cells, so they lie along dual cell edges in both directions, E x E.

Nodes{Dual}    -- C x C
Nodes{Primal}  -- E x E
XEdges{Dual}   -- E x C
YEdges{Dual}   -- C x E
XEdges{Primal} -- C x E
YEdges{Primal} -- E x C

This is important when considering the interpolation and differentiation operations,
both of which must account for the relative indexing between E and C. The
key to remember is that E[1] is bounded by C[1] and C[2] in our indexing convention,
so that, for either interpolation or differentiation:
      E[i] <- C[i], C[i+1] for i in 1:N-1
and
      C[i] <- E[i-1], E[i] for i in 2:N-1
and
      C[i] <- C[i] for i in 1:N
and
      E[i] <- E[i] for i in 1:N-1
=#

module Fields

import Base: @propagate_inbounds, show, summary, fill!

#using Compat
using FFTW
using SpecialFunctions
using Statistics

using LinearAlgebra
using SparseArrays
using Interpolations

import LinearAlgebra: mul!, ldiv!, cross, ×, dot, ⋅

import Base: parentindices
const GAMMA = MathConstants.γ

export Primal, Dual, ScalarGridData, VectorGridData, GridData,
       Points, ScalarData, VectorData, TensorData,
       celltype, griddatatype, indexshift,
       diff!,grid_interpolate!,
       curl, curl!, Curl, divergence, divergence!, Divergence,
       grad, grad!, Grad,
       laplacian, laplacian!, laplacian_symm!, plan_laplacian, plan_laplacian!,
       helmholtz, helmholtz!, plan_helmholtz, plan_helmholtz!,
       plan_intfact,plan_intfact!,Identity,
       product, product!, ∘,
       directional_derivative!, directional_derivative_conserve!, curl_cross!,
       convective_derivative!, convective_derivative_rot!,
       DDF, GradDDF,
       Regularize, RegularizationMatrix, InterpolationMatrix,
       CircularConvolution

abstract type CellType end
abstract type Primal <: CellType end
abstract type Dual <: CellType end

abstract type GridData{NX,NY,T} <: AbstractMatrix{T} end

abstract type ScalarGridData{NX,NY,T} <: GridData{NX,NY,T} end

abstract type VectorGridData{NX,NY,T} <: GridData{NX,NY,T} end

# List of scalar grid types. Each pair of numbers specifies
# the number of grid points in each direction for this data type, relative
# to the reference grid. The two pairs of numbers correspond to Primal
# and Dual versions of this grid data type.
const SCALARLIST = [ :Nodes, (-1,-1), (0,0)],
                   [ :XEdges, (0,-1), (-1,0)],
                   [ :YEdges, (-1,0), (0,-1)]

# List of collection grid types. This specifies the list of different types
# to create regularization and coordinates routines for.
const VECTORLIST = [:Edges, [:Primal]],
                   [:Edges, [:Dual]],
                   [:NodePair, [:Primal,:Dual]],
                   [:NodePair, [:Dual,:Primal]]

const TENSORLIST = [:EdgeGradient, [:Dual,:Primal]],
                   [:EdgeGradient, [:Primal,:Dual]]

function othertype end

macro othertype(celltype, k)
    esc(quote
        Fields.othertype(::$celltype) = $k
        Fields.othertype(::Type{$celltype}) = $k
    end)
end

@othertype Primal Dual
@othertype Dual Primal
@othertype CellType CellType

#@wraparray ScalarGridData data 2

include("fields/fieldmacros.jl")
include("fields/scalargrid.jl")

# Generate the scalar grid field types and associated functions
for (wrapper,primaldn,dualdn) in SCALARLIST
    @eval @scalarfield $wrapper $primaldn $dualdn
end

include("fields/collections.jl")

include("fields/basicoperations.jl")
include("fields/points.jl")

#CollectedData = Union{EdgeGradient{R,S,NX,NY,T},NodePair{R,S,NX,NY,T}} where {R,S,NX,NY,T}

include("fields/physicalgrid.jl")
include("fields/operators.jl")

end


text:
module NavAbilityCaesarExtensions
using RoME
using IncrementalInference

# Primary Struct Extensions
include("./entities/NavAbilityBlobStore.jl")
include("./entities/NavAbilityDFG.jl")
export NavAbilityDFG
export NavAbilityBlobStore

# Function Extensions
import DistributedFactorGraphs: addVariable!, getVariable, getVariables # Import required to extend
import DistributedFactorGraphs: addFactor!, getFactor, getFactors # Import required to extend
import DistributedFactorGraphs: exists # Import required to extend
import DistributedFactorGraphs: addDataEntry!, addDataBlob!, getDataEntry, getDataBlob # Import required to extend
import IncrementalInference: solveGraph! # Import required to extend
export addVariable!, getVariable, getVariables
export addFactor!, getFactor, getFactors, packFactor
export solveGraph!
export exists
export addDataEntry!, addDataBlob!, getDataEntry, getDataBlob
include("./services/BlobStore.jl")
include("./services/Variables.jl")
include("./services/Factor.jl")
include("./services/Solve.jl")
include("./services/Node.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 679551a2-320f-11ec-14eb-db9ace5d282c
md"Many ideas for the style guide here are borrowed fom the great [FastAI library](https://docs.fast.ai/dev/style.html). So, some of the code style may seem similar.

There is good reason to borrow coding style ideas from the fastai library. Some of them are 👇

* Avoid too much scrolling while reading the code.
* Convey information as fast as possible.
* Keep things simple.
* Spend less time on thinking what a particular piece of code does i.e. the name of a code piece should make it clear about the intention of that code."

# ╔═╡ 9504a064-0376-4eb5-8600-9d7f36830ade
md"## Naming"

# ╔═╡ 524fb413-e364-43da-b56e-41b14a90c04e
md"""
* Module names are single words with the starting letter as uppercase.
* Functions and variables are lowercase single words.
  * For example,`makedocs`. This is to keep thing simple and easier for the user to remember the names. For example, the end user should not get confused about wether the name had an underscore or not.
* Struct names are single word and start with an upper case.

* Commonly used things should be shorter.
  * For example,
    * "object" is named `o`.
    * "image" is named `img`.
    * input is denoted with x.
  * Look in the "Abbr" section for a list of abbrevations used.
    * Feel free to add any abbrevations which are missed."""

# ╔═╡ 95bb2db3-a85d-4ec1-9ae5-588aac0ff197
md"## Layout"

# ╔═╡ f6fffc14-20d2-482a-97a9-32c25ce91e96
md"""
* One line of code should implement one complete idea, where possible
Generally therefore an if part and its 1-line statement should be on one line
Using the ternary operator x = y if a else b can help with this guideline.

* If a 1-line function body comfortably fits on the same line as the function section, feel free to put them together.
  * In such cases, wherever possible, the functions can also be denoted in the "expression" format like `func(x) = x+1`.

* If you’ve got a bunch of 1-line functions doing similar things, they don’t need a blank line between them

  * These can be grouped together."""

# ╔═╡ fe9ef97b-48d5-472d-948b-3648ac9438fd
begin
f(x:: Int8) = x+1
function f(x:: Float16) x*2 end
end

# ╔═╡ 3fdec90a-4b42-400e-8e8c-93624716b015
md"""
* Aim to align statement parts that are conceptually similar. It allows the reader to quickly see how they’re different.

  * For example
    ```
    if foo == 0 x = f(x)
       else     x = f(x) end
    ```"""

# ╔═╡ Cell order:
# ╠═679551a2-320f-11ec-14eb-db9ace5d282c
# ╠═9504a064-0376-4eb5-8600-9d7f36830ade
# ╠═524fb413-e364-43da-b56e-41b14a90c04e
# ╠═95bb2db3-a85d-4ec1-9ae5-588aac0ff197
# ╠═f6fffc14-20d2-482a-97a9-32c25ce91e96
# ╠═fe9ef97b-48d5-472d-948b-3648ac9438fd
# ╠═3fdec90a-4b42-400e-8e8c-93624716b015


text:
using BenchmarkTools

function calculate_color(diff::Float64, title::Symbol)
    threshold = improvement_thresholds(title)
    if diff > threshold
        color = :red
    elseif diff < -threshold
        color = :green
    else
        color = :normal
    end
    return color
end

function trial_to_dict(trial::BenchmarkTools.Trial)
    d = Dict{Symbol, Union{Vector, Int}}()
    d[:times] = trial.times
    d[:gctimes] = trial.gctimes
    d[:memory] = trial.memory
    d[:allocs] = trial.allocs

    return d
end

function construct_trial_group(trials::Vector{BenchmarkTools.Trial},
                               component_names::Vector{Symbol})
    d = Dict{Symbol, Dict}()
    for (trial, component) in zip(trials, component_names)
        d[component] = trial_to_dict(trial)
    end
    return d
end

function improvement_thresholds(title::Symbol)
    title_str = string(title)
    if contains(title_str, "gctimes")
        threshold = 100.
    elseif contains(title_str, "memory")
        threshold = 5.
    elseif contains(title_str, "allocs")
        threshold = 5.
    else # if the title is for times
        threshold = 10.
    end
    return threshold
end

function stat_maps()
    d = Dict()
    d[:min] = minimum
    d[:median] = median
    d[:mean] = mean
    d[:max] = maximum
    d[:identity] = identity

    return d
end


text:
module Libxc
using Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Libxc not installed properly, run Pkg.build(\"Libxc\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

include("version.jl")
include("xc.jl")
include("functionals_info.jl")

export Functional
include("Functional.jl")

export evaluate_lda!
export evaluate_gga!
include("evaluate.jl")

end # module


text:
using SoftRobots
using MultiPoly
# using ProfileView

function test_snake()
    world = SoftRobots.World3D()
    world_state = Dict{SoftRobots.Object, SoftRobots.ObjectState}()

    robot, robot_state = SoftRobots.tetrahedron()
    push!(world.objects, robot)
    world_state[robot] = robot_state

    x, y, z = generators(MPoly{Float64}, :x, :y, :z)
    terrain = SoftRobots.FixedObject()
    terrain_state = SoftRobots.FixedObjectState(-(0.1 + 0.5x^2  + 0.5y^2 - 1z))
    push!(world.objects, terrain)
    world_state[terrain] = terrain_state

    # Profile.clear()
    SoftRobots.update!(world, world_state, 0.001)
    @time @profile for j = 1:1000
        SoftRobots.update!(world, world_state, 0.001)
        SoftRobots.draw(robot, world_state[robot])
    end
    # ProfileView.view()
end

test_snake()


text:
# Polyharmonic functions used below.
function polyharmonic_phi1(r) return r; end
function polyharmonic_phi2(r) return r>1e-16 ? r^2 * log(r) : r; end
function polyharmonic_phi3(r) return r^3; end
function polyharmonic_phi4(r) return r>1e-16 ? r^4 * log(r) : r^3; end

# Uses polyharmonic interpolation to interpolate at x using values u at centers c.
# sizes are:
# x is d x nx
# c is d x nc
# u is nc
# (optional) ord determines the type of basis function. default is 3. 1,2,3,4 are available
function polyharmonic_interp(x, c, u, ord=3)
    dim = size(x,1);
    nx = size(x,2); # number of places to evaluate

    # If no centers are provided, this is meaningless, but return zero.
    if length(c)==0
        return zeros(nx);
    end

    nc = size(c,2); # number of centers

    # Constant if only one center is given
    if nc == 1
        return u[1] .* ones(nx);
    end

    if ord == 1
        phi = polyharmonic_phi1;
    elseif ord == 2
        phi = polyharmonic_phi2;
    elseif ord == 3
        phi = polyharmonic_phi3;
    else
        phi = polyharmonic_phi4;
    end

    # build A matrix
    A = zeros(nc + dim+1, nc + dim+1);
    for i=1:nc
        for j=1:nc
            A[i,j] = phi(norm(c[:,j] - c[:,i]));
        end
    end
    # build B matrix (within A matrix)
    A[1:nc, nc+1] .= 1;
    A[nc+1, 1:nc] .= 1;
    A[1:nc, (nc+2):(nc+dim+1)] = c';
    A[(nc+2):(nc+dim+1), 1:nc] = c;

    # build RHS vector
    b = vcat(u, zeros(dim+1));

    wv = A\b; # coefficients for the polyharmonic

    sol = zeros(nx);
    for i=1:nx
        # evaluate at x[i]
        sol[i] = wv[nc+1];
        for k=1:dim
            sol[i] += wv[nc+k+1] * x[k,i];
        end
        for j=1:nc
            sol[i] += wv[j]*phi(norm(c[:,j] - x[:,i]));
        end
    end

    return sol;
end


text:
function count_upto_one()
    counter = 0
    accumulated = 0.0
    while true
        accumulated += rand()
        counter += 1
        if accumulated >= 1.0
            break
        end
    end
    return counter
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using PastaQ
using ITensors
using Test
using LinearAlgebra
using Random

@testset "randomstate: MPS" begin
  N = 5
  χ = 3

  # Real-valued with randpars
  ψ = randomstate(Float64, N; χ=χ)
  @test length(ψ) == N
  @test maxlinkdim(ψ) == χ
  for j in 1:length(ψ)
    @test eltype(ψ[j]) == Float64
  end
  # Real-valued with circuit
  ψ = randomstate(Float64, N; χ=χ, alg="circuit")
  @test length(ψ) == N
  @test maxlinkdim(ψ) == χ
  for j in 1:length(ψ)
    @test eltype(ψ[j]) == Float64
  end
  # Complex-valued with randpars
  ψ = randomstate(N; χ=χ)
  @test length(ψ) == N
  @test maxlinkdim(ψ) == χ
  for j in 1:length(ψ)
    @test eltype(ψ[j]) == Complex{Float64}
  end

  # TODO: complex-valued with circuit: not implemented

  # Normalization
  ψ = randomstate(Float64, N; χ=χ, normalize=true)
  @test norm(ψ)^2 ≈ 1
  ψ = randomstate(Float64, N; χ=χ, alg="circuit", normalize=true)
  @test norm(ψ)^2 ≈ 1
  ψ = randomstate(N; χ=χ, normalize=true)
  @test norm(ψ)^2 ≈ 1
end

@testset "randomstate: LPDO" begin
  N = 5
  χ = 3
  ξ = 2

  # Real-valued with randpars
  ρ = randomstate(Float64, N; χ=χ, ξ=ξ, mixed=true)
  @test length(ρ) == N
  @test maxlinkdim(ρ.X) == χ
  for j in 1:length(ρ)
    @test eltype(ρ.X[j]) == Float64
  end
  ρ_mat = PastaQ.array(ρ)
  @test sum(abs.(imag(diag(ρ_mat)))) ≈ 0.0 atol = 1e-10
  @test all(real(eigvals(ρ_mat)) .≥ 0)
  # Complex-valued with randpars
  ρ = randomstate(N; χ=χ, ξ=ξ, mixed=true)
  @test length(ρ) == N
  @test maxlinkdim(ρ.X) == χ
  for j in 1:length(ρ)
    @test eltype(ρ.X[j]) == Complex{Float64}
  end
  ρ_mat = PastaQ.array(ρ)
  @test sum(abs.(imag(diag(ρ_mat)))) ≈ 0.0 atol = 1e-10
  @test all(real(eigvals(ρ_mat)) .≥ 0)

  # TODO: randomLPDO with circuit + thermal state not implemented

  # Normalization
  ρ = randomstate(Float64, N; χ=χ, normalize=true, mixed=true)
  @test tr(ρ) ≈ 1
  ρ = randomstate(N; χ=χ, normalize=true, mixed=true)
  @test tr(ρ) ≈ 1

  ρ = randomstate(N; χ=χ, ξ=ξ)
  ρ_mat = PastaQ.array(ρ)
  @test sum(abs.(imag(diag(ρ_mat)))) ≈ 0.0 atol = 1e-10
  @test all(real(eigvals(ρ_mat)) .≥ 0)
  @test real(tr(ρ_mat * ρ_mat)) < 1
end

@testset "randomprocess: MPO" begin
  N = 5
  χ = 3

  # Real-valued with randpars
  U = randomprocess(Float64, N; χ=χ, mixed=false)
  @test length(U) == N
  @test maxlinkdim(U) == χ
  for j in 1:length(U)
    @test eltype(U[j]) == Float64
  end
  # Complex-valued with randpars
  U = randomprocess(N; χ=χ, mixed=false)
  @test length(U) == N
  @test maxlinkdim(U) == χ
  for j in 1:length(U)
    @test eltype(U[j]) == Complex{Float64}
  end
end

@testset "randomprocess: LPDO" begin
  N = 4
  χ = 3
  ξ = 2

  # Real-valued with randpars
  Λ = randomprocess(Float64, N; χ=χ, ξ=ξ, mixed=true)
  @test length(Λ) == N
  @test maxlinkdim(Λ.X) == χ
  for j in 1:length(Λ)
    @test eltype(Λ.X[j]) == Float64
  end
  # Complex-valued with randpars
  Λ = randomprocess(N; χ=χ, ξ=ξ, mixed=true)
  @test length(Λ) == N
  @test maxlinkdim(Λ.X) == χ
  for j in 1:length(Λ)
    @test eltype(Λ.X[j]) == Complex{Float64}
  end
  # Complex-valued with randpars
  Λ = randomprocess(N; χ=χ, ξ=ξ)
  @disable_warn_order begin
    Λmat = PastaQ.array(Λ)
    @test real(tr(Λmat * Λmat)) < 1
  end
end

@testset "initialization given a state" begin
  # Complex-valued with randpars
  N = 3
  χ = 3

  Ψ0 = randomstate(N; χ=χ)
  Ψ = randomstate(Ψ0)
  for j in 1:length(N)
    @test firstind(Ψ[j]; tags="Site") == firstind(Ψ0[j]; tags="Site")
  end

  ρ = randomstate(Ψ0; mixed=true)
  for j in 1:length(N)
    @test firstind(ρ.X[j]; tags="Site", plev=0) == firstind(Ψ0[j]; tags="Site")
  end

  ## 1. Given a LPDO
  ρ0 = randomstate(N; mixed=true, χ=χ)

  # to MPS
  Ψ = randomstate(ρ0)
  for j in 1:length(N)
    @test firstind(Ψ[j]; tags="Site") == firstind(ρ0.X[j]; tags="Site", plev=0)
  end
  # to LPDO
  ρ = randomstate(ρ0; mixed=true)
  for j in 1:length(N)
    @test firstind(ρ.X[j]; tags="Site", plev=0) == firstind(ρ0.X[j]; tags="Site", plev=0)
  end

  ψ = randomstate(N; χ = χ, normalize = true)
  @test norm(ψ) ≈ 1.0
  ρ = randomstate(N; χ = χ, ξ = 2, normalize = true)
  @test tr(ρ) ≈ 1.0

end

@testset "initialization given a process" begin
  # Complex-valued with randpars
  N = 3
  χ = 3

  # 1. Given a MPO
  U0 = randomprocess(N; mixed=false, χ=χ)

  # to MPO
  U = randomprocess(U0; mixed=false)
  for j in 1:length(N)
    @test firstind(U[j]; tags="Site", plev=0) == firstind(U0[j]; tags="Site", plev=0)
  end
  # to Choi
  Λ = randomprocess(U0; mixed=true)
  Γ0 = LPDO(PastaQ.unitary_mpo_to_choi_mps(U0))
  for j in 1:length(N)
    @test firstind(Λ.X[j]; tags="Input") == firstind(Γ0.X[j]; tags="Input")
    @test firstind(Λ.X[j]; tags="Output") == firstind(Γ0.X[j]; tags="Output")
  end

  # 1. Given a Choi
  Λ0 = randomprocess(N; mixed=true, χ=χ)
  # to MPO
  U = randomprocess(Λ0; mixed=false)
  for j in 1:length(N)
    @test firstind(U[j]; tags="Site", plev=0).id == firstind(Λ0.X[j]; tags="Output").id
    @test firstind(U[j]; tags="Site", plev=1).id == firstind(Λ0.X[j]; tags="Input").id
  end

  # to Choi
  Λ = randomprocess(Λ0; mixed=true)
  for j in 1:length(N)
    @test firstind(Λ.X[j]; tags="Input") == firstind(Λ0.X[j]; tags="Input")
    @test firstind(Λ.X[j]; tags="Output") == firstind(Λ0.X[j]; tags="Output")
  end

  U = randomprocess(N; χ = χ, normalize = true)
  Λ = MPO(PastaQ.unitary_mpo_to_choi_mps(U))
  @test tr(Λ) ≈ 1 << N
  Λ = randomprocess(N; χ = χ, ξ = 2, normalize = true)
  @test tr(Λ) ≈ 1 << N
end


text:
function chebyshev_quadrature(node_type::Function,n::S,domain = [1.0,-1.0]) where {S <: Integer}

  if n <= 0

    error("The number of nodes must be positive.")

  end

  nodes   = zeros(n)
  weights = zeros(n)

  if node_type == chebyshev_nodes # weighting function is (1.0 .- nodes.^2).^(-1/2) -- (first kind)

    for i = 1:n
      nodes[i] = (domain[1]+domain[2])/2.0 - cos((i-0.5)*pi/n)*(domain[1]-domain[2])/2.0
      weights[i] = pi/n
    end

    return nodes, weights

  elseif node_type == chebyshev_extrema # weighting function is (1.0 .- nodes.^2).^(1/2) -- (second kind)

    for i = 1:n
      nodes[i] = (domain[1]+domain[2])/2.0 - cos((i-1)*pi/(n-1))*(domain[1]-domain[2])/2.0
      weights[i] = (pi/(n-1))*sin(((i-1)/(n-1))*pi)^2
    end

    return nodes, weights

  else

    error("This node type is not supported")

  end

end


text:
using LinearAlgebra
using Plots

struct parameters
    gamma::Float64
    x_range::Vector{Float64}
    y_range::Vector{Float64}
    max_tree_nodes::Int64
    nominal_step::Float64
    expansion_time::Float64
    RewireTime::Float64
    obstacle_grid_n::Int64
    cell_size_x::Float64
    cell_size_y::Float64
    function parameters(gamma::Float64, x_range::Vector{Float64}, y_range::Vector{Float64},
        max_tree_nodes::Int64, nominal_step::Float64, expansion_time::Float64,
        RewireTime::Float64, obstacle_grid_n::Int64)
        cell_size_x = diff(x_range)[1] / obstacle_grid_n
        cell_size_y = diff(y_range)[1] / obstacle_grid_n
        new(gamma, x_range, y_range, max_tree_nodes, nominal_step, expansion_time,
            RewireTime, obstacle_grid_n, cell_size_x, cell_size_y)
    end
end

function calc_edge_cost(p1::Vector{Float64}, p2::Vector{Float64})
    "Euclidean distance edge cost"
    return norm(p2[1:2] - p1[1:2])
end

function sample_point(parameters)
    "Sample uniformly in workspace"
    x = parameters.x_range[1] + (parameters.x_range[2] - parameters.x_range[1]) * rand()
    y = parameters.y_range[1] + (parameters.y_range[2] - parameters.y_range[1]) * rand()
    return [x, y]
end

function plot_tree(tree)
    x_data = []
    y_data = []

    for node in tree
        if node.parent_id != 0
            push!(x_data, [tree[node.parent_id].pose[1], node.pose[1]])
            push!(y_data, [tree[node.parent_id].pose[2], node.pose[2]])
        end
    end
    plot!(x_data, y_data, marker=(:circle,2), legend=false, color="grey")
    plot!([tree[1].pose[1]], [tree[1].pose[2]], marker=(:circle,5), legend=false, color="red")
end

function intersected_cells(points::Vector{Vector{Float64}}, params::parameters)
    # Min point p1, max point x2
    min_x = params.x_range[2]
    max_x = params.x_range[1]
    min_y = params.y_range[2]
    max_y = params.y_range[1]

    for pt in points
        pt[1] < min_x ? min_x = pt[1] : nothing
        pt[1] > max_x ? max_x = pt[1] : nothing
        pt[2] < min_y ? min_y = pt[2] : nothing
        pt[2] > max_y ? max_y = pt[2] : nothing
    end

    # Get indices of occ cells in obstacle grid
    id_x_min = 1 + clamp(floor(Int64, (min_x - params.x_range[1]) / params.cell_size_x), 0, params.obstacle_grid_n - 1);
    id_y_min = 1 + clamp(floor(Int64, (min_y - params.y_range[1]) / params.cell_size_y), 0, params.obstacle_grid_n - 1);
    id_x_max = 1 + clamp(floor(Int64, (max_x - params.x_range[1]) / params.cell_size_x), 0, params.obstacle_grid_n - 1);
    id_y_max = 1 + clamp(floor(Int64, (max_y - params.y_range[1]) / params.cell_size_y), 0, params.obstacle_grid_n - 1);

    occ_cells = Vector{Vector{Int64}}()
    for i = id_x_min:id_x_max, j = id_y_min:id_y_max
        push!(occ_cells, [i, j])
    end
    return occ_cells
end


text:
using LinearAlgebra

# coefficients from PD of π functions
π = [70 -30 0;
     -30 140 -50;
     0 -50 50]
# Force Vector
f = Float64[0 100 60]
F = transpose(f)

# Displacement vector
u = π\F
println("Nodal displacement vector (mm):")
display(u)


text:
"""
    shotlines(target,rays)

Returns the number of hit and a list of rays that interserct with the target.
"""
function shotlines(target::CSGGenerator{Float64},rays::Vector)

    pang =0
    shots=[]
    #Vis.draw(t)
    #Vis.draw(t, debug=true, resolution = (800, 600))
    for r in rays
        if OpticSim.surfaceintersection(target(),r) !=EmptyInterval{Float64}()
            pang +=1
            push!(shots,r)
            #Vis.draw!(r,rayscale=0.2)
            #Vis.draw!(surfaceintersection(target(),r), color=:red)
        end
    end
    return pang,shots
end


text:
function play_mh_clean{T<:Integer}(ncur::T=3, ncar::T=1)
    ncar < ncur || throw(DomainError())
    pick = rand(1:ncur)
    isstickwin = pick <= ncar
    pick = rand(1:(ncur-2))
    if isstickwin # remove initially picked car from consideration
        pick += 1
    end
    isswitchwin = pick <= ncar
    return (isstickwin, isswitchwin)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module LogisticRegression

using LinearAlgebra

"""
    sigmoid(z)

Compute the sigmoid of z.
"""
@inline function sigmoid(z::Number)
    return 1.0 / (1.0 + exp(-z))
end

"""
    logistic_loss(y, a)

Calculate logistic loss or cross-entropy loss for sample.
"""
@inline function logistic_loss(y::Number, a::Number)
    return -1.0 * (y * log(a) + (1.0 - y) * log(1.0 - a))
end

"""
    evaluate_linear!(z, W, b, x)

Evaluate linear layer, i.e. calculate `z = W*x + b`.
"""
@inline function evaluate_linear!(
    z::AbstractVector,
    W::AbstractMatrix,
    b::AbstractVector,
    x::AbstractVector)
    mul!(z, W, x)
    z .+= b
    return z
end

"""
    propagate!(J, dw, db, W, b, X, Y, z, a, dz)

Forward- and backpropagation.

Given set of samples (X, Y), calculate total cost and partial derivatives with
respect to w and b. X and Y must be iterables with length of m, where m is the
number of samples. Algoritm writes to J, dw, db, z, a, and dz. The last three
auxiliary vectors for intemediate results and they size needs to equal b.
"""
function propagate!(
    J::AbstractVector,
    dW::AbstractMatrix,
    db::AbstractVector,
    W::AbstractMatrix,
    b::AbstractVector,
    X::AbstractVector{V},
    Y::AbstractVector,
    z::AbstractVector,
    a::AbstractVector,
    dz::AbstractVector
) where {V<:AbstractVector}
    @assert length(X) == length(Y)
    m = length(X)
    @assert m > 0
    fill!(J, 0.)
    fill!(dW, 0.0)
    fill!(db, 0.0)
    for (x, y) in zip(X, Y)
        z = evaluate_linear!(z, W, b, x)
     @. a = sigmoid(z)
     @. dz = a - y
     @. J += logistic_loss(y, a)
     @. dW += x' * dz
     @. db += dz
    end
 @. J /= m
 @. dW /= m
 @. db /= m
    return
end

"""
    optimize!(J, W, b, X, Y; num_iterations = 1000, learning_rate = 1.0)

Given initial model (J, W, b) and data (X, Y), optimize model parameters using
gradient descent method. Model optimization is done in place.
"""
function optimize!(
    J::AbstractVector,
    W::AbstractMatrix,
    b::AbstractVector,
    X::AbstractVector{T},
    Y::AbstractVector;
    num_iterations::Integer = 1000,
    learning_rate::Float64 = 1.0,
) where {T<:AbstractVector}
    dW = zero(W)
    db = zero(b)
    z = zero(b)
    a = zero(b)
    dz = zero(b)
    for i in 1:num_iterations
        propagate!(J, dW, db, W, b, X, Y, z, a, dz)
     @. W = W - learning_rate * dW
     @. b = b - learning_rate * db
    end
    return J, W, b
end

"""
    predict(W, b, x)

Given model (W, b) and sample x, predict the class y.
"""
@inline function predict(W, b, x)
    return 1 * (sigmoid.(W*x .+ b) .> 0.5)
end

end # module


text:
using .ArcadeLearningEnvironment
using .ReinforcementLearningEnvironments

function atari_env_factory(
    name,
    state_size,
    n_frames,
    max_episode_steps = 100_000;
    seed = nothing,
    repeat_action_probability = 0.25,
)
    AtariEnv(;
        name = string(name),
        grayscale_obs = true,
        noop_max = 30,
        frame_skip = 4,
        terminal_on_life_loss = false,
        repeat_action_probability = repeat_action_probability,
        max_num_frames_per_episode = n_frames * max_episode_steps,
        color_averaging = false,
        full_action_space = false,
        seed = seed,
    ) |>
    StateOverriddenEnv(
        ResizeImage(state_size...),  # this implementation is different from cv2.resize https://github.com/google/dopamine/blob/e7d780d7c80954b7c396d984325002d60557f7d1/dopamine/discrete_domains/atari_lib.py#L629
        StackFrames(state_size..., n_frames),
    ) |>
    StateCachedEnv |>
    RewardOverriddenEnv(r -> clamp(r, -1, 1))
end

for f in readdir(@__DIR__)
    if f != splitdir(@__FILE__)[2]
        include(f)
    end
end


text:
include("Load.jl")
include("Technology.jl")

mutable struct MarketOperator <: PassiveAgent
    name::String

    Match::Function # Match technology with power flow
    ClearMarket::Function # Market clearing
    Payment::Function

    function MarketOperator(name::String)
        self = new()
        self.name = name

        self.Match = function(technology::Technology, allPowerFlows::Vector{PowerFlow})
            for powerFlow in allPowerFlows
                if powerFlow.from == technology
                    return powerFlow
                end
            end
        end

        self.ClearMarket = function(allTechnologies::Array{Technology,1}, allPowerFlows::Array{PowerFlow,1}, systemLoad::Load)
            onlineTechnologies = allTechnologies
            #NOTE:start clear the market
            marketPrice = 0
            if length(onlineTechnologies) == 0
                println("No Technology!")
                marketPrice = systemLoad.priceCap
            end

            if length(onlineTechnologies) > 0
                sort!(onlineTechnologies, by = x -> x.offerPrice)
                generatorCount = 1 #initialize generator number
                 #this is a variable indicating the energy gap of one consumer when switching generators
                maxLoad = deepcopy(systemLoad.hourlyLoad)
                residualTechnology = 0

                while generatorCount <= length(onlineTechnologies)
                    thisQuantity = onlineTechnologies[generatorCount].offerQuantity
                    maxLoad += - thisQuantity
                    if maxLoad <= 0
                        residualTechnology = -maxLoad
                        marketPrice = onlineTechnologies[generatorCount].offerPrice
                        break
                    end
                    generatorCount += 1
                end
                #println(maxLoad)
                #println("Generator count ", generatorCount)
                if maxLoad > 0
                    marketPrice = systemLoad.priceCap
                    residualTechnology = 0
                end

                for thisCount in 1:length(onlineTechnologies)
                    if thisCount < generatorCount
                        thisPowerFlow = self.Match(onlineTechnologies[thisCount], allPowerFlows)
                        if typeof(thisPowerFlow)==Nothing
                            println(onlineTechnologies[thisCount].type)
                            println(generatorCount)
                        end
                        thisPowerFlow.quantity = onlineTechnologies[thisCount].offerQuantity
                        thisPowerFlow.marketPrice = marketPrice
                    end

                    if thisCount == generatorCount
                        thisPowerFlow = self.Match(onlineTechnologies[thisCount], allPowerFlows)
                        thisPowerFlow.quantity = onlineTechnologies[thisCount].offerQuantity - residualTechnology
                        thisPowerFlow.marketPrice = marketPrice
                    end

                    if thisCount > generatorCount
                        thisPowerFlow = self.Match(onlineTechnologies[thisCount], allPowerFlows)
                        thisPowerFlow.quantity = 0.0
                        thisPowerFlow.marketPrice = marketPrice
                    end
                end
            end
            #println("price:",marketPrice)
            return marketPrice
        end

        self.Payment = function(allPowerFlows::Vector{PowerFlow}, day::Int64,hour::Int64,auctionIdentify::Symbol, systemSettings::SystemSettings)
            for powerFlow in allPowerFlows
                thisTechnology = powerFlow.from
                subsidyAmount = 0.0
                #NOTE: in the case of different auction process, update different histories
                if auctionIdentify == :virtual
                    thisTechnology.UpdateHistory!(auctionIdentify, thisTechnology.economicProperties.virtualHistory,powerFlow, hour,day,subsidyAmount,systemSettings)
                    #println("Yearly: ", size(thisTechnology.economicProperties.virtualHistory.yearlyHistory))
                end
                if auctionIdentify == :real
                    thisTechnology.UpdateHistory!(auctionIdentify,thisTechnology.economicProperties.realHistory,powerFlow, hour,day,subsidyAmount,systemSettings)
                end
            end
        end
        return self
    end
end


text:
using ValidatedNumerics
using FactCheck

facts("Tests with rational intervals") do

    a = Interval(1//2, 3//4)
    b = Interval(3//7, 9//12)

    @fact a + b --> Interval(13//14, 3//2)  # exact

    @fact sqrt(a + b) --> Interval(0.9636241116594314, 1.2247448713915892)

    X = Interval(1//3)
    @fact sqrt(X) --> Interval(0.5773502691896257, 0.5773502691896258)

end

setprecision(Interval, 64)

facts("Rounding rational intervals") do

    X = Interval(big(1)//3)
    Y = Interval(big"5.77350269189625764452e-01", big"5.77350269189625764561e-01")
    @fact sqrt(X) --> Y
end

facts("Tests with float intervals") do

    c = @floatinterval(0.1, 0.2)

    @fact isa(@floatinterval(0.1), Interval) --> true
    @fact c --> Interval(0.09999999999999999, 0.2)

    @fact widen(c) --> Interval(0.09999999999999998, 0.20000000000000004)

    @fact @floatinterval(pi) --> Interval(3.141592653589793, 3.1415926535897936)
end

facts("Testing functions of intervals") do
    f(x) = x + 0.1

    c = @floatinterval(0.1, 0.2)
    @pending f(c) --> Interval(0.19999999999999998, 0.30000000000000004)

    d = @interval(0.1, 0.2)
    @pending f(d) --> Interval(1.9999999999999998e-01, 3.0000000000000004e-01)
end

facts("Testing conversions") do
    f = @interval(0.1, 0.2)
    @pending @floatinterval(f) --> Interval(0.09999999999999999, 0.2)

    g = @floatinterval(0.1, 0.2)
    @pending @interval(g) --> Interval(9.9999999999999992e-02, 2.0000000000000001e-01)
end


text:
using LinearAlgebra, AverageShiftedHistograms, Distances
"""
    T = ind_class(N::Int)

    Given `N`, determine what class `ind` (and `rs`) should be

### Input Argument
* `N::Int`: the length of the graph signal

### Output Argument
* `T <: Unsigned`: the unsigned integer class that `ind` should be

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito) |
Translated and revised by Naoki Saito, Feb. 8, 2017
"""
function ind_class(N::Int)
    if N < 2^8-2
        T = UInt8
    elseif N < 2^16-2
        T = UInt16
    elseif N < 2^32-2
        T = UInt32
    elseif N < 2^64-2
        T = UInt64
    else
        T = UInt128             # Hope this does not happen!
    end
    return T
end

"""
    T = tag_class(jmax::Int)

Given `jmax`, determine what class `tag` should be.

### Input Argument
* `jmax::Int`: the number of levels in the recursive partitioning (j = 1:jmax)

### Output Argument
* `T <: Unsigned`: the unsigned integer class that `tag` should be

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito) |
Translated and revised by Naoki Saito, Feb. 8, 2017
"""
function tag_class(jmax::Int)
    if jmax - 1 <= 8
        T = UInt8
    elseif jmax - 1 <= 16
        T = UInt16
    elseif jmax - 1 <= 32
        T = UInt32
    elseif jmax - 1 <= 64
        T = UInt64
    else
        T = UInt128             # Hope this does not happen!
    end
    return T
end

"""
    dmatrix_flatten(dmatrix::Array{Float64,3}, flatten::Any)

Flatten dmatrix using the method specified by the string "flatten"

### Input Arguments
* `dmatrix::Array{Float64,3}`: the matrix of expansion coefficients; after this function is called, it becomes the size of (~, ~, 1).
* `flatten::Any`: the method for flattening dmatrix (see the code in details)
"""
function dmatrix_flatten(dmatrix::Array{Float64,3}, flatten::Any)

    # p-norms or p-quasinorms
    if isa(flatten, Number)
        if flatten == 1
            dmatrix = sum(abs.(dmatrix), dims = 3)
        elseif flatten == 0
            dmatrix = sum(dmatrix != 0, dims = 3)
        else
            dmatrix = (sum(abs.(dmatrix).^flatten, dims = 3)).^(1 / flatten)
        end

        # histogram
        # # # elseif strcmpi(flatten,'histogram') # not yet implemented
        # What I meant here was to implement ASH-like pdf estimation
        # followed by the entropy computation.

        # sum of absolute values = 1-norm
    elseif isa(flatten, Symbol)
        if flatten == :abs
            dmatrix = sum(abs.(dmatrix), dims = 3)
            # standard deviation
        elseif flatten == :std
            # MATLAB: dmatrix = std(dmatrix, 0, 3)
            dmatrix = std(dmatrix, dims = 3)
            # variance
        elseif flatten == :var
            dmatrix = var(dmatrix, dims = 3)
            # inverse variance
        elseif flatten == :invvar
            dmatrix = var(dmatrix, dims = 3)
            dmatrix[ abs.(dmatrix) < eps() ] = eps()
            dmatrix = dmatrix.^-1;
            # max coefficient - min coefficient
        elseif flatten == :minmax || flatten == :maxmin
            # MATLAB: dmatrix = max(dmatrix, [], 3) - min(dmatrix, [], 3)
            dmatrix = maximum(dmatrix, dims = 3) - minimum(dmatrix, dims = 3)
            # sum (no absolute value)
        elseif flatten == :sum
            dmatrix = sum(dmatrix, dims = 3)
            # the sum of the differences of consecutive coeffs
        elseif flatten == :sumdiff
            dmatrix = sum(dmatrix[:, :, 2:end] - dmatrix[:, :, 1:(end - 1)], dims = 3)
            # the sum of the absolute values of the differences of consecutive coeffs
        elseif flatten == :sumabsdiff
            dmatrix = sum(abs.(dmatrix[:, :, 2:end] - dmatrix[:, :, 1:(end - 1)]), dims = 3)
            # Shannon entropy
        elseif flatten === :entropy
            # MATLAB: p = abs(dmatrix)/norm(dmatrix[:], 'fro')
            # p = abs(dmatrix)/norm(dmatrix[:])
            # dmatrix = sum(p. * log2(p), 3)
            # Normally, we flatten to get an energy matrix, and then take entropy
            dmatrix = sum(abs.(dmatrix).^2, dims = 3)
            # now, all the column vectors have the same norm, so normalize it.
            dmatrix /= sum(dmatrix[:, 1, 1])
            # now, each column has a unit norm.
            dmatrix = sum(-dmatrix .* log2.(dmatrix), dims = 3)
            # threshold coefficients below 0.5*norm(dmatrix)/length(dmatrix) and sum
        elseif flatten == :sub
            # MATLAB: t = 0.5 * norm(dmatrix[:], 'fro') / numel(dmatrix)
            t = 0.5 * norm(dmatrix[:]) / length(dmatrix)
            dmatrix[abs.(dmatrix) .< t] = 0
            dmatrix = sum(dmatrix, dims = 3)
            # ASH pdf estimation followed by the entropy computation.
        elseif flatten == :ash
            (N, jmax, fcols) = Base.size(dmatrix)
            for j = 1:jmax
                for i = 1:N
                    # ASH pdf estimation
                    h = ash(dmatrix[i, j, :])
                    # number of sampling points (default: 500)
                    npts = length(h.density)
                    # width of the sampling range (default: max - min + 1*std)
                    width = h.rng[end] - h.rng[1]
                    # normalize pdf
                    pdf = h.density / norm(h.density, 1)
                    # estimate of differential entropy by the limiting density
                    # of discrete points (uniform sampling)
                    dmatrix[i, j, 1] = entropy(pdf, 2.0) + log2(width / npts)
                end
            end
            # flatten dmatrix by slicing
            dmatrix = reshape(dmatrix[:, :, 1], N, jmax, 1)
            # default (1-norm)
        else
            @warn("the specified flatten symbol $(flatten) is not recognized; hence we assume it as 1-norm.")
            dmatrix = sum(abs.(dmatrix), dims = 3)
        end
    else
        @warn("the specified flatten argument $(flatten) is neither of number nor of symbol type; hence we assume it as 1-norm.")
        dmatrix = sum(abs.(dmatrix), dims = 3)
    end
    return dmatrix
end # of function dmatrix_flatten!

"""
function ldb_discriminant_measure(p::Vector{Float64}, q::Vector{Float64}; dm::Symbol = :KLdivergence)

Discriminat measure in LDB.

### Input Arguments
* `p,q::Vector{Float64}`: probability mass functions.
* `dm::Symbol`: discriminant measure. Options: `:KLdivergence`(default),
    `:Jdivergence`, `:l1`, `:l2`, and `:Hellinger`.
"""
function ldb_discriminant_measure(p::Vector{Float64}, q::Vector{Float64}; dm::Symbol = :KLdivergence)
    @assert all(p .>= 0) && all(q .>= 0)
    @assert length(p) == length(q)
    if dm == :KLdivergence
        ind = findall((p .> 0) .& (q .> 0))
        return Distances.kl_divergence(p[ind], q[ind])
    elseif dm == :Jdivergence
        ind = findall((p .> 0) .& (q .> 0))
        return Distances.kl_divergence(p[ind], q[ind]) + Distances.kl_divergence(q[ind], p[ind])
    elseif dm == :l1
        return norm(p - q, 1)
    elseif dm == :l2
        return norm(p - q, 2)
    elseif dm == :Hellinger
        return hellinger(p, q)
    else
        error("This discriminat measure $(dm) is not supported! ")
    end
end

"""
    function dmatrix_ldb_flatten(dmatrix::Array{Float64,3}...; dm::Symbol = :KLdivergence)

Flatten dmatrices using the LDB method; after this function is called, it returns
a matrix of size (~, ~, 1).

### Input Arguments
* `dmatrix::Array{Float64,3}`: the matrix of LDB expansion coefficients in one class.
* `dm::Symbol`: discriminant measure. Options: `:KLdivergence` (default),
    `:Jdivergence`, `:l1`, `:l2`, and `:Hellinger`.

### Example Usage:
`dmatrix_ldb_flatten(dmatrix1, dmatrix2, dmatrix3)`,
each argument is the expansion coefficient matrix of a class of signals. It uses
the default discriminant measure KL divergence to flatten these matrices.
In other words, it flattens these expansion coefficent matrices by computing and
summing "statistical distances" among them.
"""
function dmatrix_ldb_flatten(dmatrix::Array{Float64,3}...; dm::Symbol = :KLdivergence)
    C = length(dmatrix) # number of signal classes
    if C < 2
        error("Input should contain at least two classes of signals.")
    end
    N, jmax, _ = Base.size(dmatrix[1])
    res = zeros(N, jmax)
    for u = 1:(C - 1), v = (u + 1):C
        for j = 1:jmax
            for i = 1:N
                h1 = ash(dmatrix[u][i, j, :])
                p = h1.density / norm(h1.density, 1)
                h2 = ash(dmatrix[v][i, j, :])
                q = h2.density / norm(h2.density, 1)
                res[i, j] += ldb_discriminant_measure(p, q; dm = dm)
            end
        end
    end
    res = reshape(res[:, :, 1], N, jmax, 1)
    return res
end

"""
    costfun = cost_functional(cfspec::Any)

Determine the cost functional to be used by the best-basis algorithm.

### Input Argument
* `cfspec::Any`: the specification for the cost functional

### Output Argument
* `costfun::Function`: the cost functional (as a function_handle)
"""
function cost_functional(cfspec::Any)
    if isa(cfspec, Number)
        return function (x) return norm(x, cfspec) end
    elseif isa(cfspec, Function)
        return function (x) return cfspec(x) end
    else return function (x) return norm(x, 1) end # by default 1-norm
        # else return function (x) return norm(x, 0.1) end # by default 0.1-quasinorm
    end
end # of function cost_functional

# """
#     (dvec, levlist) = bbchange(dvec::Vector{Float64}, j::Int)
#
# Change to the new best basis
#
# ### Input Arguments
# * `dvec::Vector{Float64}`: an input coefficient matrix
# * `j::Int`: a level index
#
# ### Output Arguments
# * `dvec::Vector{Float64}`: an output coefficient matrix
# * `levlist::Vector{Int}`: an output levlist
# """
# function bbchange(dvec::Vector{Float64}, j::Int)
#     n = Base.size(dvec, 1) # assume each column corresponding to one signal
#     levlist = zeros(Int, n)
#     levlist[1] = j
#     return dvec, levlist        # might need deepcopy here.
# end

"""
    function rs_to_region(rs::Matrix{Any}, tag::Matrix{Any})

From the imformation of rs, tag from GP, compute the tag_r matrix, which have the same
size of dmatrix (expansion coefficient matrix). Each element indicates the place of the
coefficient in the expansion tree.

### Input Arguments
* `rs::Matrix{Any}`: rs from GP, showing information of the partition tree
* `tag::Matrix{Any}`: tag from GP, indicating coefficients tag

### Output Arguments
* `tag_r::Matrix{UInt64}`: showing information of the partition tree, same size as dmatrix
"""
function rs_to_region(rs::Matrix{Int}, tag::Matrix{Int})
  (m,n) = size(tag)
  tag_r = zeros(m,n)
  for j = 1:(n-1)
    regioncount = count(!iszero, rs[:,j]) - 1
    for r = 1:regioncount
      rs1 = rs[r,j]
      rs3 = rs[r+1,j]
      s = rs3-rs1
      if s==1 #the region and subregion have only one element
        tag_r[rs1,j+1] = 2*tag_r[rs1,j]
      elseif s > 1
        # rs2 marks the start of the second subregion
        rs2 = rs1 + 1
        while rs2 < rs3 && tag[rs2, j+1]!=0
          rs2 = rs2 + 1;
        end

        # the parent region is a copy of the subregion
        if rs2 == rs3
          tag_r[rs1:rs3-1,j+1] = 2*tag_r[rs1:rs3-1,j]

        # the parent region has 2 child regions
        else
        tag_r[rs1:rs2-1,j+1] = 2*tag_r[rs1:rs2-1,j]
        tag_r[rs2:rs3-1,j+1] = 2*tag_r[rs2:rs3-1,j] .+ 1
        end
      end
    end
  end
return Matrix{Int}(tag_r)
end

"""
function nonorth2relerror(nonorth::array{Float64,1},B::array{Float64,2})

Given a vector 'nonorth' of non-orthonormal expansion coefficients and
the matrix 'B' such that B*nonorth is the original signal, return a
vector of relative approximation errors when retaining the 1,2,...,N
largest coefficients in magnitude.

### Input argument
* `nonorth`:    a vector of non-orthonormal expansion coefficients
* `B`:          the matrix whose such that B*nonorth is the original signal

### Output argument
* `relerror`:   a vector of relative approximation errors

"""
function nonorth2relerror(nonorth::Array{Float64,1},B::Array{Float64,2})

    # sort the expansion coefficients
    IX = sortperm(nonorth, by = abs, rev = true)

    # generate a matrix where column j contains the j largest coefficients
    matrix0 = Matrix(UpperTriangular(repeat(nonorth[IX],1,length(IX))))
    matrix = deepcopy(matrix0)
    matrix[IX,:] = matrix0

    # the original signal and a matrix of reconstructions
    f = B*nonorth
    recons = B*matrix

    # the relative errors
    relerror = sum((repeat(f,1,length(IX)) - recons).^2, dims = 1)'.^(1/2)./norm(f,2)
    return relerror
end

"""
function orth2relerror(orth)

Given a vector 'orth' of orthonormal expansion coefficients, return a
vector of relative approximation errors when retaining the 1,2,...,N
largest coefficients in magnitude.

### Input argument
* `orth`:    a vector of orthonormal expansion coefficients

### Output argument
* `relerror`:   a vector of relative approximation errors

"""
function orth2relerror(orth::Array{Float64,1})
    # sort the coefficients
    orth = sort(orth.^2, rev = true)

    #compute the relative errors
    relerror = ((abs.(sum(orth) .- cumsum(orth))).^(1/2))/sum(orth).^(1/2)

    return relerror
end

## NGWP utils functions
# project vector x onto matrix A's column space, given A is full column rank
function Proj(x, A)
    y = try
        pinv(A' * A) * (A' * x)
    catch err
        if err != Nothing
            A' * x
        end
    end
    return A * y
end

function wavelet_perp_Matrix(w,A)
    N, k = size(A)
    B = zeros(N, k - 1)
    if k == 1
        return
    end
    tmp = A'*w
    if tmp[1] > 1e-6
        M = Matrix{Float64}(I, k-1, k-1)
        return A*vcat((-tmp[2:end] ./ tmp[1])', M)
    end
    if tmp[end] > 1e-6
        M = Matrix{Float64}(I, k-1, k-1)
        return A*vcat(M,(-tmp[1:end-1] ./ tmp[end])')
    end
    ind = findall(tmp .== maximum(tmp))[1]
    rest_ind = setdiff([i for i in 1:k], ind)
    M = Matrix{Float64}(I, k-1, k-1)
    B = A * vcat(vcat(M[1:ind-1,:], (-tmp[rest_ind] ./ tmp[ind])'), M[ind:end, :])
    return B
end

# Gram-Schmidt Process Orthogonalization
function gram_schmidt(A; tol = 1e-12)
    # Input: matirx A
    # Output: orthogonalization matrix of A's column vectors

    # Convert the matrix to a list of column vectors
    a = [A[:,i] for i in 1:size(A,2)]

    # Start Gram-Schmidt process
    q = []
    complement_dim = 0
    for i = 1:length(a)
        qtilde = a[i]
        for j = 1:(i - 1)
            qtilde -= (q[j]' * a[i]) * q[j]
        end
        if norm(qtilde) < tol
            complement_dim = size(A,2) - i + 1
            break
        end
        push!(q, qtilde / norm(qtilde))
    end
    Q = zeros(size(A,1), length(q))
    for i = 1:length(q)
        Q[:,i] .= q[i]
    end
    return Q, complement_dim
end

"""
    mgslp(A::Matrix{Float64}; tol::Float64 = 1e-12, p::Float64 = 1.0)

Modified Gram-Schmidt Process Orthogonalization with ℓᵖ pivoting algorithm (MGSLp)

# Input Arguments
- `A::Matrix{Float64}`: whose column vectors are to be orthogonalized.

# Output Argument
- `A::Matrix{Float64}`: orthogonalization matrix of A's column vectors based on ℓᵖ pivoting.
"""
function mgslp(A::Matrix{Float64}; tol::Float64 = 1e-12, p::Float64 = 1.0)
    n = size(A, 2)
    # Convert the matrix to a list of column vectors
    a = [A[:, i] for i in 1:n]

    # Start modified Gram-Schmidt process
    q = []
    v = copy(a)
    vec_lp_norm = [norm(v[i], p) for i in 1:n]
    complement_dim = 0
    for i = 1:n
        # Pivoting based on minimum ℓᵖ-norm
        idx = findmin(vec_lp_norm)[2]
        v[i], v[idx + i - 1] = v[idx + i - 1], v[i]
        # Check the linear dependency
        if norm(v[i]) < tol
            complement_dim = n - i + 1
            break
        end
        qtilde = v[i] / norm(v[i], 2)
        vec_lp_norm = []
        for j = (i + 1):n
            v[j] .-= (qtilde' * v[j]) * qtilde
            push!(vec_lp_norm, norm(v[j], p))
        end
        push!(q, qtilde)
    end
    Q = zeros(size(A,1), length(q))
    for i = 1:length(q)
        Q[:,i] .= q[i]
    end
    return Q, complement_dim
end


text:
module EC2

using LightXML
using AWS.Crypto
using Requests
using AWS.AWSEnv
using AWS
using Compat

import AWS: AbstractAWSType, AWSError, AWSResponse, _ename, _parse_from_xml, _prep_params, aws_execute

const EC2Error = AWSError
const EC2Response = AWSResponse

parse_from_xml{T<:AbstractAWSType}(::Type{T}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd::Union{LightXML.XMLElement,Void}) = _parse_from_xml(T, types, pd; member_pfx="item")
prep_params{T<:AbstractAWSType}(obj::T) = _prep_params(obj; member_pfx="")

ec2_execute(env::AWSEnv, action::String, params::Vector{Tuple}=Vector{Tuple}()) = aws_execute(env, action, "ec2", EC2ApiVer, params)
ec2_execute(env::AWSEnv, action::String, req::AbstractAWSType) = ec2_execute(env, action, prep_params(req))
ec2_execute(env::AWSEnv, action::String, req::Void) = ec2_execute(env, action)

include("ec2_types.jl")
include("ec2_operations.jl")
include("ec2_spec.jl")

export EC2Error, EC2Response

end # module EC2