text:
function area(box::IntervalBox{2,T}) where {T}
    return prod([box[i].hi - box[i].lo for i = 1:2])
end

function extend_edge_quadrature_to_area_quadrature(
    funcs,
    sign_conditions,
    height_dir,
    lo,
    hi,
    x0::V,
    w0,
    quad1d::ReferenceQuadratureRule{T},
) where {T,V<:AbstractVector}

    @assert length(funcs) == length(sign_conditions)

    lower_dim = length(x0)
    dim = lower_dim + 1
    quad = TemporaryQuadrature(T, dim)

    extended_funcs = [x -> f(extend(x0, height_dir, x)) for f in funcs]
    roots = roots_and_ends(extended_funcs, lo, hi)
    for j = 1:length(roots)-1
        xc = 0.5 * (roots[j] + roots[j+1])
        if sign_conditions_satisfied(extended_funcs, xc, sign_conditions)
            p, w = transform(quad1d, roots[j], roots[j+1])
            extended_points = extend(x0, height_dir, p)
            update!(quad, extended_points, w0 * w)
        end
    end
    return quad
end

function extend_edge_quadrature_to_area_quadrature(
    funcs,
    sign_conditions,
    height_dir,
    lo,
    hi,
    x0::M,
    w0::V,
    quad1d::ReferenceQuadratureRule{T},
) where {T,V<:AbstractVector,M<:AbstractMatrix}

    @assert length(funcs) == length(sign_conditions)
    lower_dim, npoints = size(x0)
    nweights = length(w0)
    @assert npoints == nweights

    dim = lower_dim + 1
    quad = TemporaryQuadrature(T, dim)

    for i = 1:npoints
        x0p = view(x0, :, i)
        w0p = w0[i]
        lq = extend_edge_quadrature_to_area_quadrature(
            funcs,
            sign_conditions,
            height_dir,
            lo,
            hi,
            x0p,
            w0p,
            quad1d,
        )
        update!(quad, lq.points, lq.weights)
    end
    return quad
end

function extend_edge_quadrature_to_area_quadrature(
    funcs,
    sign_conditions,
    height_dir,
    interval,
    edgequad,
    quad1d,
)

    return extend_edge_quadrature_to_area_quadrature(
        funcs,
        sign_conditions,
        height_dir,
        interval.lo,
        interval.hi,
        edgequad.points,
        edgequad.weights,
        quad1d,
    )
end

function cut_area_quadrature(
    func,
    grad,
    sign_condition,
    box,
    quad1d,
    recursionlevel,
    maxlevels,
    numsplits,
)
    heightdir = height_direction(grad, box)
    issuitable, gradsign = is_suitable(heightdir, grad, box)
    if !issuitable
        splitboxes = split_box(box, numsplits)
        splitquads = [
            subdivision_area_quadrature(
                func,
                grad,
                sign_condition,
                sb,
                quad1d,
                recursionlevel + 1,
                maxlevels,
                numsplits,
            ) for sb in splitboxes
        ]
        return combine_quadratures(splitquads)
    else
        lower(x) = func(extend(x, heightdir, box[heightdir].lo))
        upper(x) = func(extend(x, heightdir, box[heightdir].hi))

        lowersign = sign(gradsign, sign_condition, false, -1)
        uppersign = sign(gradsign, sign_condition, false, +1)

        lowerbox = heightdir == 1 ? box[2] : box[1]

        edgequad = one_dimensional_quadrature(
            [lower, upper],
            [lowersign, uppersign],
            lowerbox,
            quad1d,
        )

        return extend_edge_quadrature_to_area_quadrature(
            [func],
            [sign_condition],
            heightdir,
            box[heightdir],
            edgequad,
            quad1d,
        )
    end
end

function subdivision_area_quadrature(
    func,
    grad,
    sign_condition,
    box,
    quad1d::ReferenceQuadratureRule{T},
    recursionlevel,
    maxlevels,
    numsplits,
) where {T}

    @assert sign_condition == +1 || sign_condition == -1

    if recursionlevel >= maxlevels
        error("Failed to construct an appropriate area integration rule after $maxlevels subdivisions")
        # xc = reshape(mid(box), 2, 1)
        # wt = [area(box)]
        # return TemporaryQuadrature(xc, wt)
    else
        s = sign(func, box)
        if s == +1 || s == -1
            if s != sign_condition
                return TemporaryQuadrature(T, 2)
            else
                return temporary_tensor_product(quad1d, box)
            end
        elseif s == 0
            return cut_area_quadrature(
                func,
                grad,
                sign_condition,
                box,
                quad1d,
                recursionlevel,
                maxlevels,
                numsplits,
            )
        else
            splitboxes = split_box(box, numsplits)
            splitquads = [
                subdivision_area_quadrature(
                    func,
                    grad,
                    sign_condition,
                    sb,
                    quad1d,
                    recursionlevel + 1,
                    maxlevels,
                    numsplits,
                ) for sb in splitboxes
            ]
            return combine_quadratures(splitquads)
        end
    end
end

function area_quadrature(
    func,
    grad,
    sign_condition,
    xL,
    xR,
    numqp;
    maxlevels = 5,
    numsplits = 2,
)

    box = IntervalBox(xL, xR)
    quad1d = ReferenceQuadratureRule(numqp)

    tempquad = subdivision_area_quadrature(
        func,
        grad,
        sign_condition,
        box,
        quad1d,
        0,
        maxlevels,
        numsplits,
    )
    return QuadratureRule(tempquad)
end

# function area_quadrature(
#     poly::InterpolatingPolynomial{1},
#     interpgrad,
#     sign_condition,
#     xL,
#     xR,
#     numqp;
#     maxlevels = 5,
#     numsplits = 2,
# ) where {NF,B,T}
#
#     update_interpolating_gradient!(interpgrad,poly)
#     return area_quadrature(
#         poly,
#         interpgrad,
#         sign_condition,
#         xL,
#         xR,
#         numqp,
#         maxlevels = maxlevels,
#         numsplits = numsplits,
#     )
# end

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EmpiricalVariogram(X, z₁, z₂=z₁; [optional parameters])

Computes the empirical (a.k.a. experimental) omnidirectional
(cross-)variogram from data locations `X` and values `z₁` and `z₂`.

    EmpiricalVariogram(spatialdata, var₁, var₂=var₁; [optional parameters])

Alternatively, compute the (cross-)variogram for the variables
`var₁` and `var₂` stored in a `spatialdata` object.

    EmpiricalVariogram(partition, var₁, var₂=var₁; [optional parameters])

Alternatively, compute the (cross-)variogram on a `partition` of the data.

## Parameters

  * nlags - number of lags (default to 20)
  * maxlag - maximum lag (default to half of maximum lag of data)
  * distance - custom distance function (default to Euclidean distance)

See also: [`DirectionalVariogram`](@ref)
"""
struct EmpiricalVariogram
  abscissa::Vector{Float64}
  ordinate::Vector{Float64}
  counts::Vector{Int}

  function EmpiricalVariogram(sdata, var₁, var₂, nlags, maxlag, distance)
    # retrieve relevant parameters
    npts = npoints(sdata)
    hmax = maxlag ≠ nothing ? maxlag : diagonal(boundbox(sdata)) / 2

    # sanity checks
    @assert (var₁, var₂) ⊆ keys(variables(sdata)) "invalid variable names"
    @assert nlags > 0 "number of lags must be positive"
    @assert npts > 1 "variogram requires at least 2 points"
    @assert hmax > 0 "maximum lag distance must be positive"

    # lookup variables as vectors
    z₁ = sdata[var₁]; z₂ = sdata[var₂]

    # bin (or lag) size
    Δh = hmax / nlags

    # lag sums and counts
    sums   = zeros(nlags)
    counts = zeros(Int, nlags)

    # preallocate memory for coordinates
    xi = MVector{ndims(sdata),coordtype(sdata)}(undef)
    xj = MVector{ndims(sdata),coordtype(sdata)}(undef)

    @inbounds (
    for j in 1:npts
      coordinates!(xj, sdata, j)
      for i in j+1:npts
        coordinates!(xi, sdata, i)

        # evaluate lag and (cross-)variance
        h = evaluate(distance, xi, xj)
        v = (z₁[i] - z₁[j])*(z₂[i] - z₂[j])

        # bin (or lag) where to accumulate result
        lag = ceil(Int, h / Δh)

        if lag ≤ nlags && !ismissing(v) && !isnan(v)
          sums[lag] += v
          counts[lag] += 1
        end
      end
    end
    )

    # variogram abscissa
    abscissa = range(Δh/2, stop=hmax - Δh/2, length=nlags)

    # variogram ordinate
    ordinate = (sums ./ counts) / 2
    ordinate[counts .== 0] .= NaN

    new(abscissa, ordinate, counts)
  end
end

EmpiricalVariogram(sdata::AbstractData, var₁::Symbol, var₂::Symbol=var₁;
                   nlags=20, maxlag=nothing, distance=Euclidean()) =
  EmpiricalVariogram(sdata, var₁, var₂, nlags, maxlag, distance)

function EmpiricalVariogram(X::AbstractMatrix,
                            z₁::AbstractVector,
                            z₂::AbstractVector=z₁;
                            kwargs...)
  if z₁ === z₂
    EmpiricalVariogram(PointSetData(Dict(:z₁=>z₁), X), :z₁; kwargs...)
  else
    EmpiricalVariogram(PointSetData(Dict(:z₁=>z₁, :z₂=>z₂), X), :z₁, :z₂; kwargs...)
  end
end

function EmpiricalVariogram(partition::SpatialPartition,
                            var₁::Symbol, var₂::Symbol=var₁; kwargs...)
  # retain spatial data with at least 2 points
  filtered = Iterators.filter(d -> npoints(d) > 1, partition)

  @assert !isempty(filtered) "invalid partition of spatial data"

  sdata, _ = iterate(filtered)
  γ = EmpiricalVariogram(sdata, var₁, var₂; kwargs...)
  for sdata in Iterators.drop(filtered, 1)
    γiter = EmpiricalVariogram(sdata, var₁, var₂; kwargs...)
    merge!(γ, γiter)
  end

  γ
end

"""
    DirectionalVariogram(spatialdata, direction, var₁, var₂=var₁; [optional parameters])

Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in
`spatialdata` along a given `direction`.

Optional parameters include the parameters for `EmpiricalVariogram` and the parameters
for `DirectionPartitioner`.

### Notes

A `DirectionalVariogram` is just a function that first partitions the `spatialdata`
using a `DirectionPartitioner` and then passes the result to the corresponding
`EmpiricalVariogram` constructor.

See also: [`EmpiricalVariogram`](@ref), [`DirectionPartitioner`](@ref)
"""
function DirectionalVariogram(sdata::S, direction::NTuple,
                              var₁::Symbol, var₂::Symbol=var₁;
                              tol=1e-6, kwargs...) where {S<:AbstractData}
  partitioner = DirectionPartitioner(direction; tol=tol)
  EmpiricalVariogram(partition(sdata, partitioner), var₁, var₂; kwargs...)
end

"""
    merge!(γₐ, γᵦ)

Merge the empirical variogram `γₐ` with the empirical variogram `γᵦ`
assuming that both have the same abscissa.
"""
function merge!(γₐ::EmpiricalVariogram, γᵦ::EmpiricalVariogram)
  yₐ = γₐ.ordinate
  yᵦ = γᵦ.ordinate
  nₐ = γₐ.counts
  nᵦ = γᵦ.counts

  for i in eachindex(yₐ)
    if nᵦ[i] > 0
      if nₐ[i] > 0
        yₐ[i] = (yₐ[i]*nₐ[i] + yᵦ[i]*nᵦ[i]) / (nₐ[i] + nᵦ[i])
      else
        yₐ[i] = yᵦ[i]
      end
    end
  end

  @. nₐ = nₐ + nᵦ

  nothing
end

"""
    values(γ)

Returns the center of the bins, the mean squared differences divided by 2
and the number of squared differences at the bins for a given empirical
variogram `γ`.

## Examples

Plotting empirical variogram manually:

```julia
julia> x, y, n = values(γemp)
julia> plot(x, y, label="variogram")
julia> bar!(x, n, label="histogram")
```
"""
Base.values(γ::EmpiricalVariogram) = γ.abscissa, γ.ordinate, γ.counts

ext: jl


text:
module StaggeredAcoustic
using LinearAlgebra
using SparseArrays
import sbp
using sbp.Sparse: diag_block_matrix_2x2, block_matrix_2x2, block_matrix,
                  block_matrix_insert
using sbp.Staggered: build_all_operators_2d,
                   build_covariant_basis_mm,
                   build_covariant_basis_pm,
                   build_covariant_basis_mp
using sbp.Metrics: build_jacobian,
                   build_contravariant_basis,
                   contravariant_metric_tensor_to_matrix,
                   ContravariantMetricTensor

struct AcousticOperators
        nx::Int64
        ny::Int64
        np::Int64
        n1::Int64
        n2::Int64
        Jp::AbstractArray
        J1::AbstractArray
        J2::AbstractArray
        Jip::AbstractArray
        Ji1::AbstractArray
        Ji2::AbstractArray
        Pc1::AbstractArray
        P1c::AbstractArray
        Pc2::AbstractArray
        P2c::AbstractArray
        D1p::AbstractArray
        D2p::AbstractArray
        D1::AbstractArray
        D2::AbstractArray
        Hp::AbstractArray
        H1::AbstractArray
        H2::AbstractArray
        Gp::ContravariantMetricTensor
        G1::ContravariantMetricTensor
        G2::ContravariantMetricTensor
        B1::AbstractArray
        B2::AbstractArray
end

function init_operators(nx::Int64, ny::Int64,
                        build_operators::Function,
                        fx::AbstractArray, fy::AbstractArray)
        pp, mm, pm, mp = build_all_operators_2d(build_operators, nx, ny)
        np = (nx + 1) * (ny + 1)
        n1 = nx * (ny + 1)
        n2 = (nx + 1) * ny
        ap = build_covariant_basis_mm(fx, fy, mm, pm)
        a1 = build_covariant_basis_pm(fx, fy, mm, pm, mp)
        a2 = build_covariant_basis_mp(fx, fy, mm, pm, mp)
        Jp = build_jacobian(ap)
        J1 = build_jacobian(a1)
        J2 = build_jacobian(a2)
        bp = build_contravariant_basis(Jp, ap)
        b1 = build_contravariant_basis(J1, a1)
        b2 = build_contravariant_basis(J2, a2)
        Jip = spdiagm(0 => 1 ./Jp)
        Ji1 = spdiagm(0 => 1 ./ J1)
        Ji2 = spdiagm(0 => 1 ./ J2)
        Jp = spdiagm(0 => Jp)
        J1 = spdiagm(0 => J1)
        J2 = spdiagm(0 => J2)
        Gp = sbp.Metrics.build_contravariant_metric_tensor(bp)
        Gp = contravariant_metric_tensor_to_matrix(Gp)
        G1 = sbp.Metrics.build_contravariant_metric_tensor(b1)
        G1 = contravariant_metric_tensor_to_matrix(G1)
        G2 = sbp.Metrics.build_contravariant_metric_tensor(b2)
        G2 = contravariant_metric_tensor_to_matrix(G2)
        P1c = pm.Px
        P2c = mp.Py
        Pc1 = mm.Px
        Pc2 = mm.Py
        D1p = mm.Dx
        D2p = mm.Dy
        D1 = pm.Dx
        D2 = mp.Dy
        Hp = mm.Hx * mm.Hy
        H1 = pm.Hx * pm.Hy
        H2 = mp.Hx * mp.Hy
        Hxi = spdiagm(0 => 1 ./ diag(pm.Hx))
        Hyi = spdiagm(0 => 1 ./ diag(mp.Hy))
        B1 = Hxi * pm.Bx
        B2 = Hyi * mp.By
        return AcousticOperators(nx, ny, np, n1, n2, Jp, J1, J2, Jip, Ji1, Ji2,
                                     Pc1, P1c, Pc2, P2c, D1p, D2p, D1, D2, Hp,
                                     H1, H2, Gp, G1, G2, B1, B2)

end

function grid_points(field::String, nx::Int64, ny::Int64)
        if field == "p"
                return nx+1, ny+1
        elseif field == "v1"
                return nx+1, ny
        elseif field == "v2"
                return nx, ny+1
        elseif field == "node"
                return nx, ny
        end
end

function grids(field::String, nx::Int64, ny::Int64; pad=false::Bool)
        xp1, h = sbp.Grid.grid_xp(nx, pad=pad)
        yp1, h = sbp.Grid.grid_xp(ny, pad=pad)

        xm1, h = sbp.Grid.grid_xm(nx+1)
        ym1, h = sbp.Grid.grid_xm(ny+1)

        if field == "p"
                xm = sbp.Grid.grid_2d_x(xm1, length(ym1))
                ym = sbp.Grid.grid_2d_y(ym1, length(xm1))
                return xm, ym
        elseif field == "v1"
                xp = sbp.Grid.grid_2d_x(xp1, length(ym1))
                ym = sbp.Grid.grid_2d_y(ym1, length(xp1))
                return xp, ym
        elseif field == "v2"
                xm = sbp.Grid.grid_2d_x(xm1, length(yp1))
                yp = sbp.Grid.grid_2d_y(yp1, length(xm1))
                return xm, yp
        elseif field == "node"
                xp = sbp.Grid.grid_2d_x(xp1, length(yp1))
                yp = sbp.Grid.grid_2d_y(yp1, length(xp1))
                return xp, yp
        end
end

function pressure_norm(op::AcousticOperators)
        return op.Hp * op.Jp
end

function velocity_norm(op::AcousticOperators)
        H = diag_block_matrix_2x2(op.H1, op.H2)
        J = diag_block_matrix_2x2(op.J1, op.J2)
        return H * J
end

function energy_norm(HJp::AbstractArray, HJ::AbstractArray, Gp::AbstractArray)
        rows = [size(HJp,1), size(HJ,1)]
        E = block_matrix(rows, rows)
        E = block_matrix_insert(E, rows, rows, 1, 1, HJp)
        Gip = inv(Matrix(Gp))
        E = block_matrix_insert(E, rows, rows, 2, 2, HJ * Gip)
        return E
end

function contravariant_metric_tensor(op::AcousticOperators)
        G11 = op.Ji1 * op.P1c * op.Jp * op.Gp.b11 * op.Pc1
        G12 = op.Ji1 * op.P1c * op.Jp * op.Gp.b12 * op.Pc2
        G21 = op.Ji2 * op.P2c * op.Jp * op.Gp.b21 * op.Pc1
        G22 = op.Ji2 * op.P2c * op.Jp * op.Gp.b22 * op.Pc2
        return block_matrix_2x2(G11, G12, G21, G22)
end

function modified_contravariant_metric_tensor(op::AcousticOperators)
        G11 = op.G1.b11
        G12 = op.Ji1 * op.P1c * op.Jp * op.Gp.b12 * op.Pc2
        G21 = op.Ji2 * op.P2c * op.Jp * op.Gp.b21 * op.Pc1
        G22 = op.G2.b22
        return block_matrix_2x2(G11, G12, G21, G22)
end

function divergence(op::AcousticOperators)
        rows = [op.np]
        cols = [op.n1, op.n2]
        Ap = block_matrix(rows, cols)
        Ap1 = op.Jip * op.D1p * op.J1
        Ap2 = op.Jip * op.D2p * op.J2
        Ap = block_matrix_insert(Ap, rows, cols, 1, 1, Ap1)
        Ap = block_matrix_insert(Ap, rows, cols, 1, 2, Ap2)
        return Ap
end

function gradient(op::AcousticOperators, G::AbstractArray)
        rows = [op.n1, op.n2]
        cols = [op.np]
        Av = block_matrix(rows, cols)
        A1 = op.D1
        A2 = op.D2
        Av = block_matrix_insert(Av, rows, cols, 1, 1, A1)
        Av = block_matrix_insert(Av, rows, cols, 2, 1, A2)
        return G * Av

end

function spatial_discretization(Ap::AbstractArray, Av::AbstractArray)
        rows = [size(Ap,1), size(Av,1)]
        cols = [size(Av,2), size(Ap,2)]
        A = block_matrix(rows, cols)
        A = block_matrix_insert(A, rows, cols, 1, 2, Ap)
        A = block_matrix_insert(A, rows, cols, 2, 1, Av)
        return A
end

function pressure_sat(op::AcousticOperators, G::AbstractArray)
        rows = [op.n1, op.n2]
        cols = [op.np]
        S = block_matrix(rows, cols)
        S = block_matrix_insert(S, rows, cols, 1, 1, op.B1)
        S = block_matrix_insert(S, rows, cols, 2, 1, op.B2)
        return G * S
end

function pressure_source(op::AcousticOperators,
                         r1_s::Float64,
                         r2_s::Float64,
                         num_moment::Int64,
                         num_smoothness::Int64)
        rows = [op.np, op.n1, op.n2]
        cols = [1]
        d = pressure_source_disc(op, r1_s, r2_s, num_moment, num_smoothness)
        S = block_matrix(rows, cols)
        S = block_matrix_insert(S, rows, cols, 1, 1, d)
        return S
end

function pressure_source_disc(
                              op::AcousticOperators,
                              r1_s::Float64,
                              r2_s::Float64,
                              num_moment::Int64,
                              num_smoothness::Int64)
        r1, h1 = sbp.Grid.grid_xm(op.nx + 1)
        r2, h2 = sbp.Grid.grid_xm(op.ny + 1)
        d = sbp.Source.source_discretize_2d(r1_s, r2_s, num_moment,
                                            num_smoothness, r1, r2, h1,
                                            h2)
        return d
end

function pressure_receiver(op::AcousticOperators,
                         r1_s::Float64,
                         r2_s::Float64,
                         num_moment::Int64,
                         num_smoothness::Int64)
        rows = [op.np, op.n1, op.n2]
        cols = [1]
        r = pressure_receiver_disc(op, r1_s, r2_s, num_moment, num_smoothness)
        S = block_matrix(rows, cols)
        S = block_matrix_insert(S, rows, cols, 1, 1, r)
        return S
end

function pressure_receiver_disc(op::AcousticOperators,
                                r1_s::Float64,
                                r2_s::Float64,
                                num_moment::Int64,
                                num_smoothness::Int64)
        r1, h1 = sbp.Grid.grid_xm(op.nx + 1)
        r2, h2 = sbp.Grid.grid_xm(op.ny + 1)
        r = sbp.Source.source_discretize_2d(r1_s, r2_s, num_moment,
                                            num_smoothness, r1, r2, 1.0,
                                            1.0)

        return r
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "KaHyParBuilder"
version = v"1.1.0"

#Collection of sources required to build CMakeBuilder
sources = [
    "https://cmake.org/files/v3.15/cmake-3.15.0.tar.gz" =>
    "0678d74a45832cacaea053d85a5685f3ed8352475e6ddf9fcb742ffca00199b5",
    "https://dl.bintray.com/boostorg/release/1.70.0/source/boost_1_70_0.tar.gz" =>
    "882b48708d211a5f48e60b0124cf5863c1534cd544ecd0664bb534a4b5d506e9",
    "https://github.com/SebastianSchlag/kahypar.git" =>
    "328760c4fd4307be3c972bee814649863febc678"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/boost_1_70_0
./bootstrap.sh --prefix=$prefix/boost
apk add python-dev
./b2 -j$nproc --layout=system install --with-program_options

cd $WORKSPACE/srcdir
cd cmake-3.15.0/
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_USE_SYSTEM_LIBRARIES=Off -DBUILD_CursesDialog=Off ..
make
make install

cd $WORKSPACE/srcdir/kahypar
git submodule update --init --recursive
mkdir build
cd build
$WORKSPACE/srcdir/cmake-3.15.0/build/bin/cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_BUILD_TYPE=RELEASE -DBOOST_ROOT=$prefix/boost -DBoost_NO_SYSTEM_PATHS=ON -DBoost_USE_STATIC_LIBS=OFF -DBUILD_SHARED_LIBS=ON ..
make install.library
"""

# These are the platforms we will build for by default, unless further
platforms = [
    BinaryProvider.Linux(:x86_64, libc = :glibc, compiler_abi = CompilerABI(:gcc7)),
    #BinaryProvider.MacOS(compiler_abi = CompilerABI(:gcc7)),
    #BinaryProvider.Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libkahypar", :libkahypar)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jalving/KaHyParBuilder

max_issues_repo_name: jalving/KaHyParBuilder

max_forks_repo_name: jalving/KaHyParBuilder


text:
################################################################################
## Parameter models
################################################################################

abstract type ParamsProvider end
# interface:
#   tip_dimension(provider::ParamsProvider)
#   data_dimension(provider::ParamsProvider)
#   make_xmlelement(provider::ParamsProvider, tm::TreeModelXMLElement; ind::Int = 1)

const DEFAULT_PSS = 0.001

function data_dimension(provider::ParamsProvider)
    return tip_dimension(provider)
end

mutable struct DiffusionModel <: ParamsProvider
    Σ::AbstractArray{Float64, 2} # diffusion variance
    μ::AbstractVector{Float64} # diffusion mean
    pss::Float64 # prior sample size of diffusion mean

    function DiffusionModel(Σ::AbstractMatrix{Float64},
                            μ::AbstractVector{Float64},
                            pss::Float64)
        check_posdef(Σ)
        p = size(Σ, 1)
        if length(μ) != p
            error("The mean and variance have incompatible dimensions.")
        end
        if pss <= 0.0
            error("The prior sample size must be positive.")
        end
        return new(Σ, μ, pss)
    end
end

function DiffusionModel(Σ::AbstractMatrix{Float64})
    p = size(Σ, 1)
    return DiffusionModel(Σ, zeros(p), DEFAULT_PSS)
end

function DiffusionModel(p::Int)
    return DiffusionModel(Diagonal(ones(p)), zeros(p), DEFAULT_PSS)
end

function tip_dimension(provider::DiffusionModel)
    return length(provider.μ)
end

mutable struct IntegratedFactorModel <: ParamsProvider
    L::AbstractMatrix{Float64} # KxP loadings matrix
    λ::AbstractVector{Float64} # P-vector of residual variances

    function IntegratedFactorModel(L::AbstractMatrix{Float64},
                                   λ::AbstractVector{Float64})
        k, p = size(L)
        if length(λ) != p
            error("The loadings and residual variances have incompatible " *
                  "dimensions.")
        end
        return new(L, λ)
    end
end

function IntegratedFactorModel(k::Int, p::Int)
    return IntegratedFactorModel(zeros(k, p), ones(p))
end

function tip_dimension(provider::IntegratedFactorModel)
    return size(provider.L, 1)
end

function data_dimension(provider::IntegratedFactorModel)
    return size(provider.L, 2)
end

mutable struct ResidualVarianceModel <: ParamsProvider
    Γ::AbstractMatrix{Float64} # residual variance

    function ResidualVarianceModel(Γ::AbstractMatrix{Float64})
        check_posdef(Γ)
        return new(Γ)
    end
end

function ResidualVarianceModel(p::Int)
    return ResidualVarianceModel(Diagonal(ones(p)))
end

function tip_dimension(provider::ResidualVarianceModel)
    return size(provider.Γ, 1)
end

mutable struct JointProcessModel
    diffusion_model::DiffusionModel
    extensions::AbstractArray{<:ParamsProvider}

    function JointProcessModel(diff::DiffusionModel,
                               extensions::AbstractArray{<:ParamsProvider})
        p_diff = tip_dimension(diff)
        p_ext = 0
        for ext in extensions
            p_ext += tip_dimension(ext)
        end
        if p_diff != p_ext
            error("Incompatible dimensions between diffusion process and " *
                  "extensions")
        end
        return new(diff, extensions)
    end
end

function JointProcessModel(extensions::AbstractArray{<:ParamsProvider})
    p_ext = 0
    for ext in extensions
        p_ext += tip_dimension(ext)
    end
    return JointProcessModel(DiffusionModel(p_ext), extensions)
end

function tip_dimension(jpm::JointProcessModel)
    return tip_dimension(jpm.diffusion_model)
end

function data_dimension(jpm::JointProcessModel)
    dim = 0
    for ext in jpm.extensions
        dim += data_dimension(ext)
    end
    return dim
end

################################################################################
## Data models
################################################################################

const DEFAULT_TRAIT_NAME = "traits"

mutable struct DataModel
    taxa::AbstractVector{<:AbstractString}
    data::Matrix{Float64}
    partitions::Vector{UnitRange{Int}}
    trait_names::Vector{<:AbstractString}

    function DataModel(taxa::AbstractVector{<:AbstractString},
                       data::AbstractMatrix{Float64}, name::AbstractString)
        p = size(data, 2)
        return new(taxa, Matrix(data), [1:p], [name])
    end

    function DataModel(taxa::AbstractVector{<:AbstractString},
                       data::AbstractArray{<:AbstractArray{Float64}},
                       trait_names::AbstractArray{<:AbstractString})
        if length(data) != length(trait_names)
            error("The number of data matrices must equal the number of " *
                  "trait names")
        end
        n_mats = length(data)
        n = length(taxa)
        partitions = Vector{UnitRange}(undef, n_mats)
        stop = 0
        for i = 1:n_mats
            start = stop + 1
            n_sub, p = size(data[i])
            stop += p
            if n_sub != n
                error("All matrices must have the same number of rows as taxa.")
            end
            partitions[i] = start:stop
        end
        p_total = stop
        mat = zeros(n, p_total)
        for i = 1:n_mats
            mat[1:n, partitions[i]] .= data[i]
        end

        return new(taxa, mat, partitions, trait_names)
    end
end

function DataModel(taxa::AbstractArray{<:AbstractString},
                   data::AbstractMatrix{Float64})
    return DataModel(taxa, data, DEFAULT_TRAIT_NAME)
end

function trait_dimensions(dm::DataModel)
    return [length(x) for x in dm.partitions]
end

function add_factor_proportion(bx::BEASTXMLElement)
    trait_el = get_traitLikelihood(bx)
    if_el = get_integratedFactorModel(bx)
    fps = FactorProportionStatistic(if_el, trait_el, "factorProportion")
    add_loggable(bx, fps, already_made = false)
end


text:
module Diagnostics

export StateChecker, CFL, AdvectiveCFL, DiffusiveCFL

using CUDA
using Oceananigans
using Oceananigans.Operators

using Oceananigans: AbstractDiagnostic
using Oceananigans.Utils: TimeInterval, IterationInterval, WallTimeInterval

import Base: show
import Oceananigans: run_diagnostic!

include("state_checker.jl")
include("cfl.jl")

end


text:
module TestDeprecated

using Test, DataFrames

@testset "DataFrame!" begin
    x = [1,2,3]
    y = [4,5,6]
    @test DataFrame!(x=x, y=y, copycols=true) == DataFrame(x=x,y=y)
    df1 = DataFrame(x=x, y=y)
    df2 = DataFrame!(df1)
    @test df1 == df2
    @test df1.x === df2.x
    @test df1.y === df2.y

    a=[1,2,3]
    df = DataFrame!(:a=>a, :b=>1, :c=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!("a"=>a, "b"=>1, "c"=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!(Dict(:a=>a, :b=>1, :c=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!(Dict("a"=>a, "b"=>1, "c"=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!((x, y))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), (:x1, :x2))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), ("x1", "x2"))
    @test names(df) == ["x1", "x2"]
    @test df."x1" === x
    @test df."x2" === y

    @test_throws MethodError DataFrame!([1 2; 3 4], copycols=false)
    @test_throws MethodError DataFrame!([1 2; 3 4])
    @test_throws MethodError DataFrame!([Union{Int, Missing}, Union{Float64, Missing}],
                                        [:x1, :x2], 2)
end

end # module


text:
# compound expression
z = begin
    x = 1
    y = 2
    x + y
end

(x = 1;
 y = 2;
 x + y)

 # conditional evaluation
 function test(x, y)
    if x < y
        println("<")
        z = false
    elseif x > y
        println(">")
        z = false
    else
        println("=")
        z = true
    end
    z == true || println("no")  # if-blocks are leaky
end
test(1, 3)

begin
    x, y = 1, 2
    println(x < y ? "less" :
            x > y ? "greater" : "equal")
end

# short-circuit evaluation
false || (println("true"); true)
false && println("this won't be displayed")

function fact(n::Int)
    n >= 0 || error("n must be non-negative")
    n == 0 && return 1
    n * fact(n-1)
end
fact(5)

# repeated evaluation
i = 1
while i <= 5
    println(i)
    global i += 1
end

for i = 1:5
    if i == 1
        continue
    end
    println(i)
end

for i in 1:5
    println(i)
    if i > 2
        break
    end
end

for i ∈ ["foo", "bar"]
    println(i)
end

for i in 1:2, j in 3:4
    println((i,j))
end

# exceptions
try
    sqrt(-1)
catch e
    println("caught exception", e)
end

struct MyException <: Exception end
f(x) = x >= 0 ? x : throw(MyException())
f(-1)
f(1)

f(x) = try sqrt(x) catch sqrt(complex(x, 0)) end  # but exception is much slower than comparison
f(-1)

sqrt_second(x) = try
    sqrt(x[2])
catch y
    if isa(y, DomainError)
        sqrt(complex(x[2], 0))
    elseif isa(y, BoundsError)
        sqrt(x)
    end
end
sqrt_second([1 3])
sqrt_second(3)

fff = 1
try
    #
finally
    println(fff)
end

# tasks (aka coroutines)
function producer(c::Channel)
    put!(c, "start")
    for n in 1:4
        put!(c, 2n)
    end
    put!(c, "stop")
end
chnl = Channel(producer)
take!(chnl)

for x in Channel(producer)
    println(x)
end

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
abstract type distribution_hyper_params end
#Suff statistics must contain N which is the number of points associated with the cluster
abstract type sufficient_statistics end
abstract type distibution_sample end
import Base.copy

struct model_hyper_params
    global_hyper_params::distribution_hyper_params
    local_hyper_params
    α::Float64
    γ::Float64
    η::Float64
    global_weight::Float64
    local_weight::Float64
    total_dim::Int64
    local_dim::Int64
end

mutable struct cluster_parameters
    hyperparams::distribution_hyper_params
    distribution::distibution_sample
    suff_statistics::sufficient_statistics
    posterior_hyperparams::distribution_hyper_params
end

mutable struct splittable_cluster_params
    cluster_params::cluster_parameters
    cluster_params_l::cluster_parameters
    cluster_params_r::cluster_parameters
    lr_weights::AbstractArray{Float64, 1}
    splittable::Bool
    logsublikelihood_hist::AbstractArray{Float64,1}
end

mutable struct global_cluster
    cluster_params::splittable_cluster_params
    total_dim::Int64
    local_dim::Int64
    points_count::Int64
    clusters_count::Int64
    clusters_sub_counts::AbstractArray{Int64, 1}
end

mutable struct local_cluster
    cluster_params::splittable_cluster_params
    total_dim::Int64
    local_dim::Int64
    points_count::Int64
    global_weight::Float64
    local_weight::Float64
    globalCluster::Int64
    globalCluster_subcluster::Int64 #1 for left, 2 for right
    global_suff_stats::AbstractArray{Int64, 1}
end

mutable struct local_group
    model_hyperparams::model_hyper_params
    points::AbstractArray{Float64,2}
    labels::AbstractArray{Int64,2}
    labels_subcluster::AbstractArray{Int64,2}
    local_clusters::Vector{local_cluster}
    weights::Vector{Float64}
    group_num::Int64
end

mutable struct local_group_stats
    labels::AbstractArray{Int64,2}
    labels_subcluster::AbstractArray{Int64,2}
    local_clusters::Vector{local_cluster}
end

mutable struct hdp_shared_features
    model_hyperparams::model_hyper_params
    groups_dict::Dict
    global_clusters::Vector{global_cluster}
    weights::AbstractArray{Float64, 1}
end

function copy_local_cluster(c::local_cluster)
    return deepcopy(c)
end

function update_group_from_stats!(group::local_group, stats::local_group_stats)
    group.labels = stats.labels
    group.labels_subcluster = stats.labels_subcluster
    group.local_clusters = stats.local_clusters
end


text:
# Autogenerated wrapper script for Libepoxy_jll for x86_64-w64-mingw32
export libepoxy

using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("Libepoxy")
JLLWrappers.@declare_library_product(libepoxy, "libepoxy-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll)
    JLLWrappers.@init_library_product(
        libepoxy,
        "bin\\libepoxy-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TypedCodeUtils

using Base.Meta
import Core: Compiler
using .Compiler: widenconst, argextype
using Core: MethodInstance

"""
    Consumer

Downstream packages can change the the behaviour of `process_call` and `process_invoke`,
by overloading them for a different `Consumer`.

"""
abstract type Consumer end
struct DefaultConsumer <: Consumer end

##
# Reflection
##

"""
    Reflection

# Fields:
  - `CI`: CodeInfo
  - `mi`: MethodInfo
  - `slottypes`:
  - `sptypes`:
  - `world`:
"""
struct Reflection
    CI
    mi
    slottypes
    sptypes
    world
end

include("reflection.jl")

##
# Callsite processing
##
abstract type CallInfo end
include("process.jl")

##
# Reflection preprocessing
##
include("preprocess.jl")

##
# Utils
##
filter(f, code) = ((id, c) for (id, c) in enumerate(code) if f(c))

"""
    lookthrough(f, c)

Checks if `c isa Expr` and inspects the right hand side of assignments.
"""
function lookthrough(f, c)
    if c isa Expr
        # look through assignment
        if identify_assignment(c)
            return lookthrough(f, c.args[2])
        end
        return f(c)
    end
    return false
end

##
# Identify
##

identify_assignment(c::Expr) = c.head === :(=)
identify_invoke(c::Expr) = c.head === :invoke
identify_call(c::Expr) = c.head === :call
identify_foreigncall(c::Expr) = c.head === :foreigncall

end # module


text:
using ITensors
using ITensorNetworkAD
using OptimKit
using Zygote
using Test

@testset "optimization" begin
  @testset "Energy minimization" begin
    N = 6
    s = siteinds("S=1/2", N; conserve_qns=true)
    os = OpSum()
    for n in 1:(N - 1)
      os .+= 0.5, "S+", n, "S-", n + 1
      os .+= 0.5, "S-", n, "S+", n + 1
      os .+= "Sz", n, "Sz", n + 1
    end
    Hmpo = MPO(os, s)
    ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓")
    H = prod(Hmpo)
    ψ₀ = prod(ψ₀mps)
    # The Rayleigh quotient to minimize
    function E(H::ITensor, ψ::ITensor)
      ψdag = dag(ψ)
      return (ψdag' * H * ψ)[] / (ψdag * ψ)[]
    end
    E(ψ::ITensor) = E(H, ψ)
    ∇E(ψ::ITensor) = E'(ψ)
    fg(ψ::ITensor) = (E(ψ), ∇E(ψ))
    linesearch = HagerZhangLineSearch(;
      c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0, verbosity=0
    )
    algorithm = LBFGS(3; maxiter=30, gradtol=1e-8, linesearch=linesearch)
    ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)
    D, _ = eigen(H; ishermitian=true)
    @test E(H, ψ) < E(H, ψ₀)
    @test E(H, ψ) ≈ minimum(D)
  end
  @testset "Energy minimization (MPS)" begin
    N = 4
    χ = 4
    s = siteinds("S=1/2", N; conserve_qns=true)
    os = OpSum()
    for n in 1:(N - 1)
      os .+= 0.5, "S+", n, "S-", n + 1
      os .+= 0.5, "S-", n, "S+", n + 1
      os .+= "Sz", n, "Sz", n + 1
    end
    Hmpo = MPO(os, s)
    ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=χ)
    H = ITensors.data(Hmpo)
    ψ₀ = ITensors.data(ψ₀mps)
    # The Rayleigh quotient to minimize
    function E(H::Vector{ITensor}, ψ::Vector{ITensor})
      N = length(ψ)
      ψdag = dag.(addtags.(ψ, "bra"; tags="Link"))
      ψ′dag = prime.(ψdag)
      e = ITensor(1.0)
      for n in 1:N
        e = e * ψ′dag[n] * H[n] * ψ[n]
      end
      norm = ITensor(1.0)
      for n in 1:N
        norm = norm * ψdag[n] * ψ[n]
      end
      return e[] / norm[]
    end
    E(ψ) = E(H, ψ)
    ∇E(ψ) = E'(ψ)
    fg(ψ) = (E(ψ), ∇E(ψ))
    linesearch = HagerZhangLineSearch(; c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0)
    algorithm = LBFGS(5; maxiter=100, gradtol=1e-8, linesearch=linesearch, verbosity=0)
    ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)
    sweeps = Sweeps(5)
    setmaxdim!(sweeps, χ)
    fψmps, ψmps = dmrg(Hmpo, ψ₀mps, sweeps; outputlevel=0)
    time_Eψ = @elapsed E(ψ)
    time_∇Eψ = @elapsed E'(ψ)
    @test E(H, ψ) ≈ inner(ψmps, Hmpo, ψmps) / inner(ψmps, ψmps)
  end
end


text:
# naive implementationr

const Torus32 = Int32

mutable struct Polynomial {T} <: Number
    coef :: Array{T, 1}
end

degree(p::Polynomial) = length(p.coef) - 1

Int32PolynomialRing = Polynomial{Int32}
Torus32PolynomialRing = Polynomial{Torus32}

Base.:+(p1::Torus32PolynomialRing, p2::Torus32PolynomialRing) =
    Torus32PolynomialRing(p1.coef + p2.coef)

function Base.:*(p1::Int32PolynomialRing, p2::Torus32PolynomialRing)
    N_1 = degree(p1)
    coef = zeros(Torus32, degree(p1))
    for i in 0:N_1
        for j in 0:N_1
            if i + j < N_1 + 1
                coef[i+j+1] += p1.coef[i+1] * p2.coef[j+1]
            else
                coef[i+j-N_1] -= p1.coef[i+1] * p2.coef[j+1]
            end
        end
    end
    Torus32PolynomialRing(coef)
end

Base.:*(p1::Torus32PolynomialRing, p2::Int32PolynomialRing) = p2*p1

max_stars_repo_name: bbtatt/BBTTFHE

max_issues_repo_name: bbtatt/BBTTFHE

max_forks_repo_name: bbtatt/BBTTFHE


text:
# Constants (in atomic units) used in some of the potentials.

"""
Coulomb constant in atomic units.
https://en.wikipedia.org/wiki/Coulomb_constant
"""
const kₑ = 1.0

"""
Vacuum permittivity in atomic units.
https://en.wikipedia.org/wiki/Vacuum_permittivity
"""
const ϵ₀ = 1 / 4π


text:
correct_jwt = replace("""
eyJ0eXAiOiJKV1QiLCJhbGciOiJSUzI1NiJ9.
eyJleHAiOjE0NzM5NDgwNjAsImlhdCI6MTQ3Mzk0ODAwMCwiaXNzIjoxMjM0fQ.
uDZWWptpvEy4Dv1h_jHcqqROzs2HpE_XMcSfNWhpCuLNwnVPWUVoTYqiCPbNdO7
PgsuTMhtuF4PQN0bJdz7_Bnn0UqmHZ-P1ktFeMBT7IvfEBolIVLMTa2LJQC6UG7
_qCoY2kjLdriWiLFHznJvG6jfPHK-iX9VIolNjkiM9e4DG9Aq60UnZ_df40wZXd
696sRpgCakvIV3mQTmRv9IfOLVF9eRRD4yVvwTtYNGOqewpQqkPnm6K3ctYlQIX
kwKMynp6R-CgbwRedA4n0WAvy1o14TyZZ-QAChQUcS-OKb0ZM4z-fbG5ZSpWP7f
wsQxsZgWFIz6hodiw_q45bHYsLw
""",'\n' => "")
# Fix iat, to make sure the payload is reproducible
auth = GitHub.JWTAuth(1234, joinpath(dirname(@__FILE__), "not_a_real_key.pem");
    iat = DateTime("2016-9-15T14:00"))
key = MbedTLS.PKContext()
MbedTLS.parse_key!(key,
    read(joinpath(dirname(@__FILE__), "not_a_real_key.pem"), String))
auth = GitHub.JWTAuth(1234, joinpath(dirname(@__FILE__), "not_a_real_key.pem");
    iat = DateTime("2016-9-15T14:00"))
auth2 = GitHub.JWTAuth(1234, key; iat = DateTime("2016-9-15T14:00"))
# The validity of this token can be checked with jwt.io. This is just a smoke
# test to make sure things don't break.
@test auth.JWT == correct_jwt
@test auth2.JWT == correct_jwt


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Ccluster"
version = v"1.1.7"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/rimbach/Ccluster/archive/refs/tags/v1.1.7.tar.gz", "725ab22cf7e74afe5a5133ac75ee4a101d7b4ff5f0f25a6b74f5d9bfda8a18d5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Ccluster*/
if [[ ${target} == *musl* ]]; then
    export CFLAGS="-D_GNU_SOURCE=1 -std=c99 -O2 -funroll-loops -g"
elif [[ ${target} == *mingw* ]]; then
    sed -i -e "s#/lib\>#/$(basename ${libdir})#g" configure; extraflags=--build=MINGW${nbits}
fi
./configure --prefix=$prefix \
    --disable-pthread \
    --disable-static \
    --enable-shared \
    --with-gmp=$prefix \
    --with-mpfr=$prefix \
    --with-flint=$prefix \
    --with-arb=$prefix \
    ${extraflags}
make library -j${nproc}
make install LIBDIR=$(basename ${libdir})
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libccluster", :libccluster)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Arb_jll", uuid="d9960996-1013-53c9-9ba4-74a4155039c3"))
    Dependency(PackageSpec(name="FLINT_jll", uuid="e134572f-a0d5-539d-bddf-3cad8db41a82"))
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.1.2"; compat="6.1.2")
    Dependency(PackageSpec(name="MPFR_jll", uuid="3a97d323-0669-5f0c-9066-3539efd106a3"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Minimal server using the 'listen' syntax, starting with the 'inner' functions
const BAREHTML = "<meta http-equiv=\"Content-Type\" content=\"text/html; charset=UTF-8\">
 <title>Empty.html</title>"
using Sockets
using WebSockets
import WebSockets.handle
const LOCALIP = string(Sockets.getipaddr())
const PORT = 8080
const BODY =  "
Press F12

ws = new WebSocket(\"ws://$LOCALIP:$PORT\")

ws.onmessage = function(e){console.log(e.data)}

ws.send(\"Browser console says hello!\")
                "

function coroutine(ws)
    while isopen(ws)
        data, = readguarded(ws)
        s = String(data)
        if s == ""
            writeguarded(ws, "Goodbye!")
            break
        end
        println("Received: ", s)
        writeguarded(ws, "Hello! Send empty message to exit, or just leave.")
    end
end

function gatekeeper(req, ws)
    orig = WebSockets.origin(req)
    println("\nOrigin:", orig, "    Target:", target(req), "    subprotocol:", subprotocol(req))
    if occursin(LOCALIP, orig)
        coroutine(ws)
    elseif orig == ""
        @info "Non-browser clients don't send Origin. We liberally accept the update request in this case."
        coroutine(ws)
    else
        @warn "Inacceptable request"
    end
end

handler(req) = replace(BAREHTML, "" => BODY) |> WebSockets.Response
handler_wrap = WebSockets.RequestHandlerFunction(handler)

@info("Browser > $LOCALIP:$PORT , F12> console > ws = new WebSocket(\"ws://$LOCALIP:$PORT\") ")

const SERVER = Sockets.listen(Sockets.InetAddr(parse(IPAddr, LOCALIP), PORT))

task = @async try
    WebSockets.HTTP.listen(LOCALIP, PORT, server = SERVER, readtimeout = 0 ) do http
        if WebSockets.is_upgrade(http.message)
            WebSockets.upgrade(gatekeeper, http)
        else
            handle(handler_wrap, http)
        end
    end
catch err
    # Add your own error handling code; HTTP.jl sends error code to the client.
    @info err
    @info stacktrace(catch_backtrace())
end
@info("To stop serving: close(SERVER)")

nothing


text:
module UnitBlockAngular

export UnitBlockAngularMatrix, UnitBlockAngularFactor

using LinearAlgebra

# Base Interface
import Base:
    size, getindex

import LinearAlgebra.mul!

"""
    UnitBlockAngularMatrix{Tv}

Unit block-angular matrix.

A unit block-angular matrix is a matrix of the form
```
A = [ B0  B1  B2  ...  Br ]
    |  0  e'              |
    |  0      e'          |
    |  0          ...     |
    [  0               e' ]
```
where `Bk` has size `m x nk`.
"""
struct UnitBlockAngularMatrix{Tv<:Real} <: AbstractMatrix{Tv}

    m0::Int  # number of linking constraints
    n0::Int  # Number of linking columns
    n::Int   # number of non-linking columns

    M::Int  # Total number of constraints
    N::Int  # Total number of columns
    R::Int  # Number of blocks

    B0::Matrix{Tv}         # Linking block
    B::Matrix{Tv}          # Lower blocks, concatenated
    blockidx::Vector{Int}  # Block index of each column in `B`

    function UnitBlockAngularMatrix(
        B0::AbstractMatrix{Tv},
        B::AbstractMatrix{Tv},
        R::Int, blockidx::Vector{Int}
    ) where{Tv<:Real}

        # Sanity checks
        m, n = size(B)
        m0, n0 = size(B0)
        m == m0 || throw(DimensionMismatch(
            "B has $m rows, but B0 has $(m0)"
        ))
        n == length(blockidx) || throw(DimensionMismatch(
            "B has $n columns, but blockidx has $(length(blockidx))"
        ))

        if n > 0
            imin, imax = extrema(blockidx[1:n])
            (0 <= imin <= imax <= R) || throw(DimensionMismatch(
                "Specified $R blocks but indices are in range [$imin, $imax]"
            ))
        end

        A = new{Tv}(m0, n0, n, m0 + R, n0+n, R,
            Matrix{Tv}(undef, m0, n0),
            Matrix{Tv}(undef, m, n),
            Vector{Int}(undef, n)
        )

        A.B0 .= B0
        A.B  .= B
        A.blockidx .= blockidx
        return A
    end

    function UnitBlockAngularMatrix{Tv}(m0, n0, n, R) where{Tv<:Real}
        return new{Tv}(m0, n0, n, m0+R, n0+n, R, Matrix{Tv}(undef, m0, n0), Matrix{Tv}(undef, m0, n), Vector{Int}(undef, n))
    end
end

# Base matrix interface
size(A::UnitBlockAngularMatrix) = (A.M, A.N)

function getindex(A::UnitBlockAngularMatrix{Tv}, i::Integer, j::Integer) where Tv<:Real

    # Sanity check
    ((1 <= i <= A.M) && (1 <= j <= A.N)) || throw(BoundsError())

    if (i <= A.m0) && (j <= A.n0)
        # Element is in linking block
        return A.B0[i, j]

    elseif (i <= A.m0) && (j > A.n0)
        # Element is in regular block
        return A.B[i, j-A.n0]

    elseif (i > A.m0) && (j <= A.n0)
        return zero(Tv)

    elseif (i > A.m0) && (j > A.n0)
        # Check if column j belongs to block i
        if A.blockidx[j-A.n0] == i - A.m0
            return oneunit(Tv)
        else
            return zero(Tv)
        end
    end
end

# Matrix-vector products
function LinearAlgebra.mul!(
    y::AbstractVector{Tv},
    A::UnitBlockAngularMatrix{Tv},
    x::AbstractVector{Tv},
    α::Tv, β::Tv
) where{Tv<:Real}
    # Dimensions checks
    m, n = size(A)
    n == length(x) || throw(DimensionMismatch("A has size $((m, n)) but x has length $(length(x))"))
    m == length(y) || throw(DimensionMismatch("A has size $((m, n)) but y has length $(length(y))"))

    # Compute `y0 = α * A * x + β * y0`
    # `y0 = α * B0 * x0 + α B * [x1 ... xR] + β y0`
    if iszero(β)
        # We do this explicity in case `y` has not been initialized
        # Otherwise, we may have some NaN values going around
        y .= zero(Tv)
    end

    @views x0 = x[1:A.n0]
    @views x1 = x[(A.n0+1):end]
    @views y0 = y[1:A.m0]
    @views y1 = y[(A.m0+1):end]
    mul!(y0, A.B0, x0, α, β)
    mul!(y0, A.B, x1, α,  one(Tv))

    # Compute the remaining part
    y1 .*= β
    @inbounds for (j, xj) in zip(A.blockidx, x1)
        y1[j] += α * xj
    end

    return y
end

mul!(y::AbstractVector{Tv}, A::UnitBlockAngularMatrix{Tv}, x::AbstractVector{Tv}) where{Tv<:Real} = mul!(y, A, x, one(Tv), zero(Tv))

function mul!(
    x::AbstractVector{Tv},
    At::Union{
        LinearAlgebra.Transpose{Tv, UnitBlockAngularMatrix{Tv}},
        LinearAlgebra.Adjoint{Tv, UnitBlockAngularMatrix{Tv}}
    },
    y::AbstractVector{Tv},
    α::Tv,
    β::Tv
) where{Tv<:Real}
    A = At.parent

    m, n = size(A)
    n == length(x) || throw(DimensionMismatch(
        "A has size $(size(A)) but x has size $(length(x))")
    )
    m == length(y) || throw(DimensionMismatch(
        "A has size $(size(A)) but y has size $(length(y))")
    )

    # `x = B' * y0`
    @views x0 = x[1:A.n0]
    @views x1 = x[(A.n0+1):end]
    @views y0 = y[1:A.m0]
    @views y1 = y[(A.m0+1):end]
    mul!(x0, A.B0', y0, α, β)
    mul!(x1, A.B' , y0, α, β)

    # `∀i, x[i] += y[blockidx[i]]`
    @inbounds for (i, j) in enumerate(A.blockidx)
        x1[i] += α * y1[j]
    end

    return x
end

mul!(
    x::AbstractVector{Tv},
    At::Union{
        LinearAlgebra.Transpose{Tv, UnitBlockAngularMatrix{Tv}},
        LinearAlgebra.Adjoint{Tv, UnitBlockAngularMatrix{Tv}}
    },
    y::AbstractVector{Tv}
) where{Tv<:Real} = mul!(x, At, y, one(Tv), zero(Tv))

import Tulip.TLA.construct_matrix

function construct_matrix(
    ::Type{UnitBlockAngularMatrix}, M::Int, N::Int,
    aI::Vector{Int}, aJ::Vector{Int}, aV::Vector{Tv};
    m0::Int, n0::Int, n::Int, R::Int
) where{Tv<:Real}
    A = UnitBlockAngularMatrix{Tv}(m0, n0, n, R)
    A.B .= 0
    A.B0 .= 0
    A.blockidx .= 0
    # TODO: may be more efficient to first sort indices so that
    # A is accessed in column-major order.
    for(i, j, v) in zip(aI, aJ, aV)
        if 1 <= i <= m0
            # Linking constraints
            if 1 <= j <= n0
                # Update B0
                A.B0[i, j] = v
            elseif n0 < j <= n0 + n
                # Update B
                A.B[i, j - n0] = v
            else
                error("Invalid column index $j")
            end
        elseif m0 < i <= m0 + R
            # Convexity constraint
            n0 < j || error("Non-zero coefficient ($i, $j) (should be zero)")
            isone(v) || error("Invalid coefficient A[$i, $j]=$v (should be one)")

            A.blockidx[j - n0] = i - m0
        else
            error("Invalid row index $i in matrix with $m0 linking and $R convexity constraints")
        end
    end
    return A
end

include("KKTSolver.jl")

end # module


text:
# PPO loss
function ppo_loss(π, 𝒫, 𝒟; info = Dict())
    new_probs = logpdf(π, 𝒟[:s], 𝒟[:a])
    r = exp.(new_probs .- 𝒟[:logprob])

    A = 𝒟[:advantage]
    p_loss = -mean(min.(r .* A, clamp.(r, (1f0 - 𝒫[:ϵ]), (1f0 + 𝒫[:ϵ])) .* A))
    e_loss = -mean(entropy(π, 𝒟[:s]))

    # Log useful information
    ignore() do
        info[:entropy] = -e_loss
        info[:kl] = mean(𝒟[:logprob] .- new_probs)
        info[:clip_fraction] = sum((r .> 1 + 𝒫[:ϵ]) .| (r .< 1 - 𝒫[:ϵ])) / length(r)
    end
    𝒫[:λp]*p_loss + 𝒫[:λe]*e_loss
end

function PPO(;π::ActorCritic,
     ϵ::Float32 = 0.2f0,
     λp::Float32 = 1f0,
     λe::Float32 = 0.1f0,
     target_kl = 0.012f0,
     a_opt::NamedTuple=(;),
     c_opt::NamedTuple=(;),
     log::NamedTuple=(;),
     kwargs...)

     OnPolicySolver(;agent=PolicyParams(π),
                    𝒫=(ϵ=ϵ, λp=λp, λe=λe),
                    log = LoggerParams(;dir = "log/ppo", log...),
                    a_opt = TrainingParams(;loss = ppo_loss, early_stopping = (infos) -> (infos[end][:kl] > target_kl), name = "actor_", a_opt...),
                    c_opt = TrainingParams(;loss = (π, 𝒫, D; kwargs...) -> Flux.mse(value(π, D[:s]), D[:return]), name = "critic_", c_opt...),
                    post_batch_callback = (𝒟; kwargs...) -> (𝒟[:advantage] .= whiten(𝒟[:advantage])),
                    kwargs...)
end

# PPO loss with a penalty
function lagrange_ppo_loss(π, 𝒫, 𝒟; info = Dict())
    new_probs = logpdf(π, 𝒟[:s], 𝒟[:a])
    r = exp.(new_probs .- 𝒟[:logprob])

    A = 𝒟[:advantage]
    p_loss = -mean(min.(r .* A, clamp.(r, (1f0 - 𝒫[:ϵ]), (1f0 + 𝒫[:ϵ])) .* A))
    e_loss = -mean(entropy(π, 𝒟[:s]))

    penalty = ignore(() -> Flux.softplus(𝒫[:penalty_param][1]))
    cost_loss = (penalty /  (1+penalty)) * mean(r .* 𝒟[:cost_advantage])

    # Log useful information
    ignore() do
        info[:entropy] = -e_loss
        info[:kl] = mean(𝒟[:logprob] .- new_probs)
        info[:clip_fraction] = sum((r .> 1 + 𝒫[:ϵ]) .| (r .< 1 - 𝒫[:ϵ])) / length(r)
    end
    𝒫[:λp]*p_loss + 𝒫[:λe]*e_loss + 0.5f0*cost_loss
end

function lagrange_ppo_penalty_loss(π, 𝒫, 𝒟; info = Dict())
    penalty = Flux.softplus(𝒫[:penalty_param][1])
    cur_cost = mean(𝒟[:cost])

    ignore() do
        info["penalty"] = penalty
        info["cur_cost"] = cur_cost
    end

    -penalty * (cur_cost - 𝒫[:target_cost])
end

function LagrangePPO(;π::ActorCritic,
     Vc::ContinuousNetwork, # value network for estimating cost
     ϵ::Float32 = 0.2f0,
     λp::Float32 = 1f0,
     λe::Float32 = 0.1f0,
     λ_gae = 0.95f0,
     target_kl = 0.012f0,
     penalty_init = 1f0,
     target_cost = 0.025f0,
     a_opt::NamedTuple=(;),
     c_opt::NamedTuple=(;),
     penalty_opt::NamedTuple=(;),
     cost_opt::NamedTuple=(;),
     log::NamedTuple=(;),
     kwargs...)

     𝒫=(ϵ=ϵ, λp=λp, λe=λe, penalty_param=Float32[Base.log(exp(penalty_init)-1)], target_cost=target_cost)

     OnPolicySolver(;agent=PolicyParams(π),
                    𝒫=𝒫,
                    Vc=Vc,
                    log = LoggerParams(;dir = "log/lagrange_ppo", log...),
                    param_optimizers = Dict(Flux.params(𝒫[:penalty_param]) => TrainingParams(;loss=lagrange_ppo_penalty_loss, name="penalty_", penalty_opt...)),
                    a_opt = TrainingParams(;loss = lagrange_ppo_loss, early_stopping = (infos) -> (infos[end][:kl] > target_kl), name = "actor_", a_opt...),
                    c_opt = TrainingParams(;loss = (π, 𝒫, D; kwargs...) -> Flux.mse(value(π, D[:s]), D[:return]), name = "critic_", c_opt...),
                    cost_opt = TrainingParams(;loss = (π, 𝒫, D; kwargs...) -> Flux.mse(value(π, D[:s]), D[:cost_return]), name = "cost_critic_", cost_opt...),
                    required_columns = [:return, :advantage, :logprob, :cost_advantage, :cost, :cost_return],
                    post_batch_callback = (𝒟; kwargs...) -> (𝒟[:advantage] .= whiten(𝒟[:advantage])),
                    kwargs...)
end


text:
struct ThresholdExceedance <: EVA
    data::Variable
    logscale::paramfun
    shape::paramfun
end

"""
    ThresholdExceedance(exceedances::Vector{<:Real};
        logscalecov::Vector{<:DataItem} = Vector{Variable}(),
        shapecov::Vector{<:DataItem} = Vector{Variable}())::ThresholdExceedance

Creates a ThresholdExceedance structure.

"""
function ThresholdExceedance(exceedances::Variable;
    logscalecov::Vector{<:DataItem} = Vector{Variable}(),
    shapecov::Vector{<:DataItem} = Vector{Variable}())::ThresholdExceedance

    n = length(exceedances.value)
    validatelength(n, logscalecov)
    validatelength(n, shapecov)

    logscalefun = computeparamfunction(logscalecov)
    shapefun = computeparamfunction(shapecov)

    return ThresholdExceedance(exceedances, paramfun(logscalecov, logscalefun), paramfun(shapecov, shapefun))

end

function paramindex(model::ThresholdExceedance)::Dict{Symbol,Vector{<:Int}}

    sϕ = length(model.logscale.covariate) + 1
    sξ = length(model.shape.covariate) + 1

    return Dict{Symbol,Vector{<:Int}}(
        :ϕ => Int64[k for k in (1:sϕ)],
        :ξ => Int64[k for k in .+ sϕ .+ (1:sξ)]
    )

end

function getcovariatenumber(model::ThresholdExceedance)::Int

    return sum([length(model.logscale.covariate), length(model.shape.covariate)])

end

function nparameter(model::ThresholdExceedance)::Int

    return 2 + getcovariatenumber(model)

end

function getdistribution(model::ThresholdExceedance, θ::AbstractVector{<:Real})::Vector{<:Distribution}

    @assert length(θ)==nparameter(model) "The length of the parameter vector should be equal to the model number of parameters."

    pi = paramindex(model)
    ϕ = model.logscale.fun(θ[pi[:ϕ]])
    ξ = model.shape.fun(θ[pi[:ξ]])

    σ = exp.(ϕ)

    fd = GeneralizedPareto.(σ, ξ)

    return fd

end

function getinitialvalue(model::ThresholdExceedance)::Vector{<:Real}

    y = model.data.value

    # Compute stationary initial values
    ϕ₀,ξ₀ = getinitialvalue(GeneralizedPareto, y)
    # Store them in a dictionary
    θ₀ = Dict(:ϕ => ϕ₀, :ξ => ξ₀)

    initialvalues = zeros(nparameter(model))
    pi = paramindex(model)
    for param in [:ϕ, :ξ]
        ind = pi[param][1]
        initialvalues[ind] = θ₀[param]
    end

    return initialvalues

end

"""
    showEVA(io::IO, obj::ThresholdExceedance; prefix::String = "")

Displays a ThresholdExceedance with the prefix `prefix` before every line.

"""
function showEVA(io::IO, obj::ThresholdExceedance; prefix::String = "")

    println(io, prefix, "ThresholdExceedance")
    println(io, prefix, "data :\t\t",typeof(obj.data.value), "[", length(obj.data.value), "]")
    println(io, prefix, "logscale :\t", showparamfun("ϕ", obj.logscale))
    println(io, prefix, "shape :\t\t", showparamfun("ξ", obj.shape))

end


text:
# Rebase an amount of money in a given year to December 2016 dollars
# its bad form but I assume base_year and amt are singular
function opm_inflation_adjust(year, base_year, amt)
    #base <- opm_get_inflation_index(base_year)
    #div <- sapply(year, function(x) { opm_get_inflation_index(x) })
    #amt * (base / div)
    b = opm_get_inflation_index(base_year)[1]
    mult = opm_get_inflation_index(year)
    amt .* (b ./ mult)
end

# We store the index in a hash for fast access, so breaking this out into a function
function opm_get_inflation_index(year)
    #tmp <- as.character(year)
    #sapply(tmp, function(x) { ifelse(is.na(x), NA, get(x, envir = opm_inflation_hash)) })
    [(ismissing(y) ? (missing) : opm_inflation_hash[y]) for y in year]
end

# Takes a vector of years, gives percentage change between an element and the one before
function opm_percent_change_inflation(years)
    #inflation <- opm_get_inflation_index(years)
    #opm_percent_change(inflation)
    inflation = opm_get_inflation_index(years)
    opm_percent_change(inflation)
end

# december of given year inflation index numbers.
opm_inflation_data_dec = DataFrame(
#    ~year, ~ index,
    [1970  39.8;
    1971  41.1;
    1972  42.5;
    1973  46.2;
    1974  51.9;
    1975  55.5;
    1976  58.2;
    1977  62.1;
    1978  67.7;
    1979  76.7;
    1980  86.3;
    1981  94.0;
    1982  97.6;
    1983  101.3;
    1984  105.3;
    1985  109.3;
    1986  110.5;
    1987  115.4;
    1988  120.5;
    1989  126.1;
    1990  133.8;
    1991  137.9;
    1992  141.9;
    1993  145.8;
    1994  149.7;
    1995  153.5;
    1996  158.6;
    1997  161.3;
    1998  163.9;
    1999  168.3;
    2000  174.0;
    2001  176.7;
    2002  180.9;
    2003  184.3;
    2004  190.3;
    2005  196.8;
    2006  201.8;
    2007  210.036;
    2008  210.228;
    2009  215.949;
    2010  219.179;
    2011  225.672;
    2012  229.601;
    2013  233.049;
    2014  234.812;
    2015  236.525;
    2016  241.432]
)
# julia fixes for weird data frame behavior
names!(opm_inflation_data_dec, [:year, :index])
opm_inflation_data_dec[:year] = convert.(Int, opm_inflation_data_dec[:year])

# Build a hash of our inflation data becuase we do a lot of lookups into it and table filtering
# is really slow for that many iterations.
#opm_inflation_hash <- new.env(hash = TRUE)
#invisible(apply(opm_inflation_data_dec, 1, function(x) { assign(as.character(x['year']), x['index'], envir = opm_inflation_hash) }))
#opm_inflation_data_dec <- opm_inflation_data_dec %>% mutate(pct_change = opm_percent_change(index))
opm_inflation_hash  = Dict{Int, Float64}()
for row in eachrow(opm_inflation_data_dec)
    opm_inflation_hash[row[:year]] = row[:index]
end

# For lack of a better place, putting this here.
opm_presidental_transition_years = vcat(1974, 1977, 1981, 1989, 1993, 2001, 2009)
opm_presidential_names = vcat("Ford", "Carter", "Reagan", "Bush", "Clinton", "Bush", "Obama")
opm_presidents_tibble = DataFrame(years = opm_presidental_transition_years, names = opm_presidential_names)


text:
using VulkanShaders
using Documenter

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

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


text:
"""
    classical_scaling(Δ, p=2)

Given a distance matrix `Δ` and a target dimensionality `p`, find an embedding of points
in p-dimensional Euclidean space such that the between-object distances are preserved as
well as possible.

# Examples
```julia
julia> X = rand(5, 2)
5×2 Matrix{Float64}:
 0.671412  0.326635
 0.145173  0.650203
 0.677262  0.655306
 0.070194  0.4018
 0.507032  0.0532854

julia> Δ = SMACOF.dist(X)
5×5 Matrix{Float64}:
 0.0       0.617757  0.328723  0.605899  0.318969
 0.617757  0.0       0.532114  0.259472  0.698035
 0.328723  0.532114  0.0       0.657873  0.625626
 0.605899  0.259472  0.657873  0.0       0.558829
 0.318969  0.698035  0.625626  0.558829  0.0

julia> Y = classical_scaling(Δ)
5×2 Matrix{Float64}:
  0.271143  -0.0296369
 -0.34334   -0.0931504
  0.133963  -0.328368
 -0.303247   0.163206
  0.241481   0.28795

julia> SMACOF.dist(Y)
5×5 Matrix{Float64}:
 0.0       0.617757  0.328723  0.605899  0.318969
 0.617757  0.0       0.532114  0.259472  0.698035
 0.328723  0.532114  0.0       0.657873  0.625626
 0.605899  0.259472  0.657873  0.0       0.558829
 0.318969  0.698035  0.625626  0.558829  0.0
```
"""
function classical_scaling(Δ::AbstractMatrix{T}, p = 2) where T
    n = size(Δ, 1)

    # Compute A = - Δ^2 / 2
    A = Hermitian(similar(Δ), :L)
    for j in 1:n, i in j:n
        A.data[i, j] = -Δ[i, j]^2 / 2
    end

    # Doubly center A so that the row and col sums are zero
    ravg = mean(A, dims = 2)
    tavg = mean(ravg)
    for j in 1:n, i in j:n
        A.data[i, j] = A[i, j] - ravg[i] - ravg[j] + tavg
    end

    # compute low rank approximation
    Λ, V = eigen!(A, (n - p + 1):n)
    reverse!(V, dims = 2)
    reverse!(Λ)

    for j in 1:p, i in 1:n
        V[i, j] *= sqrt(abs(Λ[j]))
    end

    return V
end


text:
@testset "Transformations" begin
    function test_reverse(T, seq)
        revseq = reverse(Kmer{T,length(seq)}(seq))
        @test String(revseq) == reverse(seq)
    end

    function test_dna_complement(seq)
        comp = complement(DNAKmer{length(seq)}(seq))
        @test String(comp) == dna_complement(seq)
    end

    function test_rna_complement(seq)
        comp = complement(RNAKmer{length(seq)}(seq))
        @test String(comp) == rna_complement(seq)
    end

    function test_dna_revcomp(seq)
        revcomp = reverse_complement(DNAKmer{length(seq)}(seq))
        @test String(revcomp) == reverse(dna_complement(seq))
    end

    function test_rna_revcomp(seq)
        revcomp = reverse_complement(RNAKmer{length(seq)}(seq))
        @test String(revcomp) == reverse(rna_complement(seq))
    end

    @testset "Reverse" begin
        for len in 1:32, _ in 1:10
            test_reverse(DNA, random_dna_kmer(len))
            test_reverse(RNA, random_rna_kmer(len))
        end

        seq = dna"AAAAAAAAAAAAAAAAAAAAAAAAAAAAGATAC"
        @test reverse(seq[(length(seq)-9):length(seq)]) == dna"CATAGAAAAA"
    end

    @testset "Complement" begin
        for len in 1:32, _ in 1:10
            test_dna_complement(random_dna_kmer(len))
            test_rna_complement(random_rna_kmer(len))
        end
    end

    @testset "Reverse Complement" begin
        for len in 1:32, _ in 1:10
            test_dna_revcomp(random_dna_kmer(len))
            test_rna_revcomp(random_rna_kmer(len))
        end
    end
end