text:
module Interpolation

using MPI
import GaussQuadrature
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Geometry
using CLIMA.Mesh.Elements
using LinearAlgebra
using StaticArrays

export InterpolationBrick, interpolate_brick!, InterpolationCubedSphere, invert_trilear_mapping_hex, interpolate_cubed_sphere!

#--------------------------------------------------------
"""
    InterpolationBrick(grid::DiscontinuousSpectralElementGrid, xres, ::Type{FT}) where FT <: AbstractFloat

This interpolation structure and the corresponding functions works for a brick, where stretching/compression happens only along the x, y & z axis.
Here x1 = X1(ξ1); x2 = X2(ξ2); x3 = X3(ξ3)

# input for the inner constructor
 - `grid` DiscontinousSpectralElementGrid
 - `xres` Resolution of the interpolation grid in x1, x2 and x3 directions
"""
struct InterpolationBrick{FT <:AbstractFloat}
    realelems::UnitRange{Int64}
    poly_order::Int

    xbnd::Array{FT,2}      # domain bounds, [2(min/max),ndim]
    xres::Array{FT,1}      # resolutions in x1, x2, x3 directions for the uniform grid

    ξ1::Vector{Vector{FT}} # unique ξ1 coordinates of interpolation points within each element
    ξ2::Vector{Vector{FT}} # unique ξ2 coordinates of interpolation points within each element
    ξ3::Vector{Vector{FT}} # unique ξ3 coordinates of interpolation points within each element

    x::Vector{Array{FT,2}} # x[elno][3,npts] -> (x1,x2,x3) coordinates of portion of interpolation grid embedded within each element

    V::Vector{Vector{FT}}  # interpolated variable within each element
#--------------------------------------------------------
    function InterpolationBrick(grid::DiscontinuousSpectralElementGrid{FT}, xres) where FT <: AbstractFloat
        T = Int
        poly_order = polynomialorder(grid)
        ndim = 3
        xbnd = zeros(FT, 2, ndim) # domain bounds (min,max) in each dimension
        for dim in 1:ndim
            xbnd[1,dim], xbnd[2,dim] = extrema(grid.topology.elemtocoord[dim,:,:])
        end
        x1g = range(xbnd[1,1], xbnd[2,1], step=xres[1])
        x2g = range(xbnd[1,2], xbnd[2,2], step=xres[2])
        x3g = range(xbnd[1,3], xbnd[2,3], step=xres[3])
        #-----------------------------------------------------------------------------------
        realelems =  grid.topology.realelems # Element (numbers) on the local processor
        Nel       = length(realelems)
        n123      = zeros(T,    ndim)     # # of unique ξ1, ξ2, ξ3 points in each cell
        xsten     = zeros(T, 2, ndim)     # x1, x2, x3 start and end for each brick element
        xbndl     = zeros(T, 2, ndim)     # location of x1,x2,x3 limits (min,max) for each brick element

        ξ1 = map( i -> zeros(FT,i), zeros(T,Nel))
        ξ2 = map( i -> zeros(FT,i), zeros(T,Nel))
        ξ3 = map( i -> zeros(FT,i), zeros(T,Nel))
        V  = map( i -> zeros(FT,i), zeros(T,Nel))
        x  = map( i -> zeros(FT,ndim,i), zeros(T,Nel)) # interpolation grid points embedded in each cell
        #-----------------------------------------------------------------------------------
        for el in 1:Nel
            for (ξ,xg,dim) in zip((ξ1,ξ2,ξ3), (x1g, x2g, x3g), 1:ndim)
                xbndl[1,dim], xbndl[2,dim] = extrema(grid.topology.elemtocoord[dim,:,el])
                xsten[1,dim], xsten[2,dim] = findfirst( temp -> temp ≥ xbndl[1,dim], xg ), findlast( temp -> temp ≤ xbndl[2,dim], xg )
                n123[dim] = xsten[2,dim] - xsten[1,dim]
                ξ[el] = [ 2 * ( xg[ xsten[1,dim] + i - 1] - xbndl[1,dim] ) / (xbndl[2,dim]-xbndl[1,dim]) -  1 for i in 1:n123[dim] ]
            end

            x_el  = zeros(FT,ndim,prod(n123))
            V[el] = zeros(FT,prod(n123))

            ctr = 1

            for k in 1:n123[3], j in 1:n123[2], i in 1:n123[1]
                x_el[1,ctr] = x1g[ xsten[1,1] + i - 1 ]
                x_el[2,ctr] = x2g[ xsten[1,2] + j - 1 ]
                x_el[3,ctr] = x3g[ xsten[1,3] + k - 1 ]
                ctr += 1
            end
            x[el] = x_el
        end # el loop
        #-----------------------------------------------------------------------------------
        return new{FT}(realelems, poly_order, xbnd, xres, ξ1, ξ2, ξ3, x, V)
    end
#--------------------------------------------------------
end # struct InterpolationBrick
#--------------------------------------------------------
"""
    interpolate_brick!(intrp_brck::InterpolationBrick, sv::AbstractArray{FT}, st_idx::T, poly_order::T) where {T <: Integer, FT <: AbstractFloat}

This interpolation function works for a brick, where stretching/compression happens only along the x, y & z axis.
Here x1 = X1(ξ1); x2 = X2(ξ2); x3 = X3(ξ3)

# input
 - `intrp_brck` Initialized InterpolationBrick structure
 - `sv` State vector
 - `st_idx` # of state vector variable to be interpolated
 - `poly_order` polynomial order for the simulation
"""
function interpolate_brick!(intrp_brck::InterpolationBrick, sv::AbstractArray{FT}, st_idx::T) where {T <: Integer, FT <: AbstractFloat}
    qm1 = intrp_brck.poly_order + 1
    Nel = length(intrp_brck.realelems)
    m1_r, m1_w = GaussQuadrature.legendre(FT,qm1,GaussQuadrature.both)
    wb = Elements.baryweights(m1_r)

    vout    = FT(0)
    vout_ii = FT(0)
    vout_ij = FT(0)

    phit = Vector{FT}(undef,qm1)

    for el in 1:Nel #-----for each element elno
        if length(intrp_brck.ξ1[el]) > 0
            l1 = length(intrp_brck.ξ1[el]); l2 = length(intrp_brck.ξ2[el]); l3 = length(intrp_brck.ξ3[el])
            lag    = @view sv[:,st_idx,el]

            phir = Elements.interpolationmatrix(m1_r, intrp_brck.ξ1[el], wb)
            phis = Elements.interpolationmatrix(m1_r, intrp_brck.ξ2[el], wb)

            for k in 1:l3 # interpolating point-by-point
                interpolationvector_1pt!(m1_r,intrp_brck.ξ3[el][k],wb,phit)

                for j in 1:l2, i in 1:l1
                    vout = 0
                    for ik in 1:qm1
                        v_ij = 0
                        for ij in 1:qm1
                            v_ii = 0
                            for ii in 1:qm1
                                @inbounds v_ii += lag[ii + (ij-1)*qm1 + (ik-1)*qm1*qm1] * phir[i,ii]
                            end # ii loop

                            @inbounds v_ij += v_ii * phis[j,ij]
                        end # ij loop
                        @inbounds vout += v_ij * phit[ik]
                    end # ik loop
                    @inbounds intrp_brck.V[el][i + (j-1)*l1 + (k-1)*l1*l2] = vout
                end # j, i loop
            end # k loop
        end
  #--------------------
  end
end
#--------------------------------------------------------
function interpolate_brick_matrix!(intrp_brck::InterpolationBrick, sv::AbstractArray{FT}, st_no::T, poly_order::T) where {T <: Integer, FT <: AbstractFloat}
    qm1 = poly_order + 1
    Nel = length(intrp_brck.realelems)
    m1_r, m1_w = GaussQuadrature.legendre(FT,qm1,GaussQuadrature.both)
    wb = Elements.baryweights(m1_r)
    Nel = length(intrp_brck.realelems)
    alpha = FT(1)
    beta  = FT(0)
    #-----for each element elno
    for el in 1:Nel
        if length(intrp_brck.ξ1[el]) > 0
            lag    = @view sv[:,st_no,el]
            phir = Elements.interpolationmatrix(m1_r, intrp_brck.ξ1[el], wb)
            phis = Elements.interpolationmatrix(m1_r, intrp_brck.ξ2[el], wb)
            phit = Elements.interpolationmatrix(m1_r, intrp_brck.ξ3[el], wb)

            sr,si = size(phir); ss,sj = size(phis); st,sk = size(phit)
            #---------------------------------------------------------------
            h1 = zeros(FT, sj*sk, sr)  # this is for testing only
            h2 = zeros(FT, sk*sr, ss)  # memory allocations within the function will be removed
            #---------------------------------------------------------------
            LinearAlgebra.BLAS.gemm!('T', 'T', alpha, reshape(lag,si, sj*sk), phir, beta, h1)
            LinearAlgebra.BLAS.gemm!('T', 'T', alpha, reshape(h1, sj, sk*sr), phis, beta, h2)
            LinearAlgebra.BLAS.gemm!('T', 'T', alpha, reshape(h2, sk, sr*ss), phit, beta, reshape(intrp_brck.V[el], sr*ss, st))
            #---------------------------------------------------------------
        end
  #--------------------
  end
end
#--------------------------------------------------------
"""
    interpolationvector_1pt!(rsrc::Vector{FT}, rdst::FT,
                             wbsrc::Vector{FT}, phi::Vector{FT}) where FT <: AbstractFloat

returns the polynomial interpolation matrix for interpolating between the points
`rsrc` (with associated barycentric weights `wbsrc`) and a single point rdst. This
function writes the result to a vector since the basis functions are calculated at a single point

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
@inline function interpolationvector_1pt!(rsrc::AbstractVector{FT}, rdst::FT,
                             wbsrc::AbstractVector{FT}, phi::AbstractVector{FT}) where FT <: AbstractFloat
    qm1 = length(rsrc)
    @assert length(phi) == qm1

    for ib in 1:qm1
        if rdst==rsrc[ib]
            phi .= FT(0)
            @inbounds phi[ib] = FT(1)
            break
        else
            @inbounds phi[ib] = wbsrc[ib] / (rdst-rsrc[ib])
        end
    end
    d = sum(phi)
    phi ./= d
    return nothing
end
#--------------------------------------------------------
"""
    InterpolationCubedSphere(grid::DiscontinuousSpectralElementGrid, vert_range::AbstractArray{FT}, lat_res::FT, long_res::FT, rad_res::FT) where {FT <: AbstractFloat}

This interpolation structure and the corresponding functions works for a cubed sphere topology. The data is interpolated along a lat/long/rad grid.

# input for the inner constructor
 - `grid` DiscontinousSpectralElementGrid
 - `vert_range` vertex range along the radial coordinate
 - `lat_res` Resolution of the interpolation grid along the latitude coordinate in radians
 - `long_res` Resolution of the interpolation grid along the longitude coordinate in radians
 - `rad_res` Resolution of the interpolation grid along the radial coordinate
"""
struct InterpolationCubedSphere{T <: Integer, FT <: AbstractFloat}

    realelems::UnitRange{Int64}
    poly_order::Int

    lat_min::FT;  long_min::FT;  rad_min::FT; # domain bounds, min
    lat_max::FT;  long_max::FT;  rad_max::FT; # domain bounds, max
    lat_res::FT;  long_res::FT;  rad_res::FT; # respective resolutions for the uniform grid

    n_lat::T; n_long::T; n_rad::T;            # # of lat, long & rad grid locations

    ξ1::Vector{Vector{FT}} # ξ1 coordinates of interpolation points within each element
    ξ2::Vector{Vector{FT}} # ξ2 coordinates of interpolation points within each element
    ξ3::Vector{Vector{FT}} # ξ3 coordinates of interpolation points within each element

    radc::Vector{Vector{FT}}  # rad coordinates of interpolation points within each element
    latc::Vector{Vector{FT}}  # lat coordinates of interpolation points within each element
    longc::Vector{Vector{FT}} # long coordinates of interpolation points within each element

    V::Vector{Vector{FT}}  # interpolated variable within each element
  #--------------------------------------------------------
    function InterpolationCubedSphere(grid::DiscontinuousSpectralElementGrid, vert_range::AbstractArray{FT}, nhor::Int, lat_res::FT, long_res::FT, rad_res::FT) where {FT <: AbstractFloat}
        poly_order = polynomialorder(grid)
        toler1 = eps(FT) * vert_range[1] * 2.0 # tolerance for unwarp function
        #toler1 = eps(FT) * 2.0 # tolerance for unwarp function
        toler2 = eps(FT) * 4.0                 # tolerance
        toler3 = eps(FT) * vert_range[1] * 4.0 # tolerance for Newton-Raphson

        lat_min,   lat_max = FT(0.0), FT(π)                 # inclination/zeinth angle range
        long_min, long_max = FT(0.0), FT(2*π)  			    # azimuthal angle range
        rad_min,   rad_max = vert_range[1], vert_range[end] # radius range

        realelems = grid.topology.realelems # Element (numbers) on the local processor
        Nel = length(realelems)

        nvert = length(vert_range) - 1              # # of elements in vertical direction
        Nel_glob = nvert * nhor * nhor * 6

        nblck = nhor * nhor * nvert
        Δh = 2 / nhor                               # horizontal grid spacing in unwarped grid

        lat_grd, long_grd, rad_grd = range(lat_min, lat_max, step=lat_res), range(long_min, long_max, step=long_res), range(rad_min, rad_max, step=rad_res)

        n_lat, n_long, n_rad = Int(length(lat_grd)), Int(length(long_grd)), Int(length(rad_grd))

        uw_grd = zeros(FT, 3, 1)
        #----------------------------------------------
        glob_ord = grid.topology.origsendorder # to account for reordering of elements after the partitioning process

        glob_elem_no = zeros(Int, nvert*length(glob_ord))

        for i in 1:length(glob_ord), j in 1:nvert
            glob_elem_no[j + (i-1)*nvert] = (glob_ord[i] - 1)*nvert + j
        end

        ξ1, ξ2, ξ3        = map( i -> zeros(FT,i), zeros(Int,Nel)), map( i -> zeros(FT,i), zeros(Int,Nel)), map( i -> zeros(FT,i), zeros(Int,Nel))
        radc, latc, longc = map( i -> zeros(FT,i), zeros(Int,Nel)), map( i -> zeros(FT,i), zeros(Int,Nel)), map( i -> zeros(FT,i), zeros(Int,Nel))

        for k in 1:n_long, j in 1:n_lat, i in 1:n_rad
            x1_grd = rad_grd[i] * sin(lat_grd[j]) * cos(long_grd[k]) # inclination -> latitude; azimuthal -> longitude.
            x2_grd = rad_grd[i] * sin(lat_grd[j]) * sin(long_grd[k]) # inclination -> latitude; azimuthal -> longitude.
            x3_grd = rad_grd[i] * cos(lat_grd[j])

            uw_grd[1], uw_grd[2], uw_grd[3] = Topologies.cubedshellunwarp(x1_grd, x2_grd, x3_grd) # unwarping from sphere to cubed shell
            rad = FT(maximum(abs.(uw_grd)))
            #--------------------------------
            x1_uw2_grd = uw_grd[1] / rad # unwrapping cubed shell on to a 2D grid (in 3D space, -1 to 1 cube)
            x2_uw2_grd = uw_grd[2] / rad
            x3_uw2_grd = uw_grd[3] / rad
            #--------------------------------
            if rad ≤ vert_range[1]       # accounting for minor rounding errors from unwarp function at boundaries
                vert_range[1] - rad < toler1 ? l_nrm = 1 :  error("fatal error, rad lower than inner radius: ", vert_range[1] - rad," $x1_grd /// $x2_grd //// $x3_grd" )
            elseif rad ≥ vert_range[end] # accounting for minor rounding errors from unwarp function at boundaries
                rad - vert_range[end] < toler1 ? l_nrm = nvert : error("fatal error, rad greater than outer radius")
            else                         # normal scenario
                l_nrm = findfirst( X -> X .- rad .> 0.0, vert_range ) - 1 # identify stack bin
            end
            #--------------------------------
            if     abs(x1_uw2_grd + 1) < toler2 # face 1 (x1 == -1 plane)
	    	    l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
    		    l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                el_glob = Int(l_nrm + (nhor-l2)*nvert + (l3-1)*nvert*nhor)
            elseif abs(x2_uw2_grd + 1) < toler2 # face 2 (x2 == -1 plane)
		        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
    		    l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                el_glob = Int(l_nrm + (l1-1)*nvert + (l3-1)*nvert*nhor + nblck*1)
            elseif abs(x1_uw2_grd - 1) < toler2 # face 3 (x1 == +1 plane)
		        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
    		    l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                el_glob = Int(l_nrm + (l2-1)*nvert + (l3-1)*nvert*nhor + nblck*2 )
            elseif abs(x3_uw2_grd - 1) < toler2 # face 4 (x3 == +1 plane)
      		    l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
	    	    l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                el_glob = Int(l_nrm + (l1-1)*nvert + (l2-1)*nvert*nhor + nblck*3)
            elseif abs(x2_uw2_grd - 1) < toler2 # face 5 (x2 == +1 plane)
	            l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
		        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                el_glob = Int(l_nrm + (l1-1)*nvert + (nhor-l3)*nvert*nhor + nblck*4 )
            elseif abs(x3_uw2_grd + 1) < toler2 # face 6 (x3 == -1 plane)
	    	    l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
		        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                el_glob = Int(l_nrm + (l1-1)*nvert + (nhor-l2)*nvert*nhor + nblck*5)
            else
                error("error: unwrapped grid does on lie on any of the 6 faces")
            end
            #--------------------------------
            el_loc = findfirst(X -> X-el_glob == 0, glob_elem_no)
            if ( el_loc ≠ nothing ) # computing inner coordinates for local elements
                ξ = invert_trilear_mapping_hex(grid.topology.elemtocoord[1,:,el_loc],
                                               grid.topology.elemtocoord[2,:,el_loc],
                                               grid.topology.elemtocoord[3,:,el_loc], uw_grd, toler3)
                push!(ξ1[el_loc],ξ[1])
                push!(ξ2[el_loc],ξ[2])
                push!(ξ3[el_loc],ξ[3])
                push!(radc[el_loc],  rad_grd[i])
                push!(latc[el_loc],  lat_grd[j])
                push!(longc[el_loc],long_grd[k])
            end
            #--------------------------------
        end

        V = [ Vector{FT}(undef, length(radc[el])) for el in 1:Nel ] # Allocating storage for interpolation variable

        return new{Int, FT}(realelems, poly_order, lat_min, long_min, rad_min, lat_max, long_max, rad_max, lat_res, long_res, rad_res,
                    n_lat, n_long, n_rad, ξ1, ξ2, ξ3, radc, latc, longc, V)
    #-----------------------------------------------------------------------------------
    end # Inner constructor function InterpolationCubedSphere
#-----------------------------------------------------------------------------------
end # structure InterpolationCubedSphere
#--------------------------------------------------------
"""
    invert_trilear_mapping_hex(X1::Array{FT}, X2::Array{FT}, X3::Array{FT}, x::Array{FT}, tol::FT) where FT <: AbstractFloat

This function computes (ξ1,ξ2,ξ3) given (x1,x2,x3) and the (8) vertex coordinates of a Hexahedron. Newton-Raphson method is used
# input
 - `X1` X1 coordinates of the (8) vertices of the hexahedron
 - `X2` X2 coordinates of the (8) vertices of the hexahedron
 - `X3` X3 coordinates of the (8) vertices of the hexahedron
 - `x` (x1,x2,x3) coordinates of the
"""
function invert_trilear_mapping_hex(X1::Array{FT}, X2::Array{FT}, X3::Array{FT}, x::Array{FT}, tol::FT) where FT <: AbstractFloat
    max_it = 10            # maximum # of iterations
    ξ      = zeros(FT,3,1) # initial guess => cell centroid

    d   = trilinear_map(ξ, X1, X2, X3) - x
    err = norm(d)
    ctr = 0
    #---Newton-Raphson iterations---------------------------
    while err > tol
        trilinear_map_IJac_x_vec!(ξ, X1, X2, X3, d)
        ξ .-= d
        d = trilinear_map(ξ, X1, X2, X3) - x
        err = norm(d)
        ctr += 1
        if ctr > max_it
            error("invert_trilinear_mapping_hex: Newton-Raphson not converging to desired tolerance after max_it = ", max_it," iterations; err = ", err,"; toler = ", tol)
        end
    end
    #-------------------------------------------------------
    return ξ
end
#--------------------------------------------------------
function trilinear_map(ξ::Array{FT}, x1v::Array{FT}, x2v::Array{FT}, x3v::Array{FT}) where FT <: AbstractFloat
    x = Array{FT}(undef,3)
    for (vert,dim) = zip((x1v,x2v,x3v),1:3)
        x[dim] = ((1 - ξ[1]) * (1 - ξ[2]) * (1 - ξ[3]) * vert[1] + (1 + ξ[1]) * (1 - ξ[2]) * (1 - ξ[3]) * vert[2] +
                  (1 - ξ[1]) * (1 + ξ[2]) * (1 - ξ[3]) * vert[3] + (1 + ξ[1]) * (1 + ξ[2]) * (1 - ξ[3]) * vert[4] +
                  (1 - ξ[1]) * (1 - ξ[2]) * (1 + ξ[3]) * vert[5] + (1 + ξ[1]) * (1 - ξ[2]) * (1 + ξ[3]) * vert[6] +
                  (1 - ξ[1]) * (1 + ξ[2]) * (1 + ξ[3]) * vert[7] + (1 + ξ[1]) * (1 + ξ[2]) * (1 + ξ[3]) * vert[8] )/ 8.0
    end
    return x
end
#--------------------------------------------------------
function trilinear_map_IJac_x_vec!(ξ::Array{FT}, x1v::Array{FT}, x2v::Array{FT}, x3v::Array{FT}, v::Array{FT}) where FT <: AbstractFloat
    Jac = MMatrix{3,3,FT,9}(undef)
    for (vert,dim) = zip((x1v,x2v,x3v),1:3)
        Jac[dim,1] = ((-1) * (1 - ξ[2]) * (1 - ξ[3]) * vert[1] + ( 1) * (1 - ξ[2]) * (1 - ξ[3]) * vert[2] +
                      (-1) * (1 + ξ[2]) * (1 - ξ[3]) * vert[3] + (+1) * (1 + ξ[2]) * (1 - ξ[3]) * vert[4] +
                      (-1) * (1 - ξ[2]) * (1 + ξ[3]) * vert[5] + (+1) * (1 - ξ[2]) * (1 + ξ[3]) * vert[6] +
                      (-1) * (1 + ξ[2]) * (1 + ξ[3]) * vert[7] + (+1) * (1 + ξ[2]) * (1 + ξ[3]) * vert[8] )/ 8.0

        Jac[dim,2] = ((1 - ξ[1]) * (-1) * (1 - ξ[3]) * vert[1] + (1 + ξ[1]) * (-1) * (1 - ξ[3]) * vert[2] +
                      (1 - ξ[1]) * (+1) * (1 - ξ[3]) * vert[3] + (1 + ξ[1]) * (+1) * (1 - ξ[3]) * vert[4] +
                      (1 - ξ[1]) * (-1) * (1 + ξ[3]) * vert[5] + (1 + ξ[1]) * (-1) * (1 + ξ[3]) * vert[6] +
                      (1 - ξ[1]) * (+1) * (1 + ξ[3]) * vert[7] + (1 + ξ[1]) * (+1) * (1 + ξ[3]) * vert[8] )/ 8.0

        Jac[dim,3] = ((1 - ξ[1]) * (1 - ξ[2]) * (-1) * vert[1] + (1 + ξ[1]) * (1 - ξ[2]) * (-1) * vert[2] +
                      (1 - ξ[1]) * (1 + ξ[2]) * (-1) * vert[3] + (1 + ξ[1]) * (1 + ξ[2]) * (-1) * vert[4] +
                      (1 - ξ[1]) * (1 - ξ[2]) * (+1) * vert[5] + (1 + ξ[1]) * (1 - ξ[2]) * (+1) * vert[6] +
                      (1 - ξ[1]) * (1 + ξ[2]) * (+1) * vert[7] + (1 + ξ[1]) * (1 + ξ[2]) * (+1) * vert[8] )/ 8.0
    end
    LinearAlgebra.LAPACK.gesv!(Jac,v)

    return nothing
end
#--------------------------------------------------------
function interpolate_cubed_sphere!(intrp_cs::InterpolationCubedSphere, sv::AbstractArray{FT}, st_no::T) where {T <: Integer, FT <: AbstractFloat}
    qm1 = intrp_cs.poly_order + 1
    Nel = length(intrp_cs.realelems)
    m1_r, m1_w = GaussQuadrature.legendre(FT,qm1,GaussQuadrature.both)
    wb = Elements.baryweights(m1_r)
    phir = MVector{qm1,FT}(undef)
    phis = MVector{qm1,FT}(undef)
    phit = MVector{qm1,FT}(undef)

    vout    = FT(0)
    vout_ii = FT(0)
    vout_ij = FT(0)

    for el in 1:Nel #-----for each element elno
        np = length(intrp_cs.ξ1[el])
        lag    = @view sv[:,st_no,el]
        for i in 1:np # interpolating point-by-point
            interpolationvector_1pt!(m1_r,intrp_cs.ξ1[el][i],wb,phir)
            interpolationvector_1pt!(m1_r,intrp_cs.ξ2[el][i],wb,phis)
            interpolationvector_1pt!(m1_r,intrp_cs.ξ3[el][i],wb,phit)
            vout = 0.0
            for ik in 1:qm1
                vout_ij = 0.0
                for ij in 1:qm1
                    vout_ii = 0.0
                    for ii in 1:qm1
                        @inbounds vout_ii += lag[ii + (ij-1)*qm1 + (ik-1)*qm1*qm1] * phir[ii]
                    end
                    @inbounds vout_ij += vout_ii * phis[ij]
                end
                @inbounds vout += vout_ij * phit[ik]
            end
            @inbounds intrp_cs.V[el][i] = vout
        end
  #--------------------
  end
  #--------------------
end
#--------------------------------------------------------
end # module interploation

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
using Documenter, Lux, Literate, Pkg

# Precompile example dependencies
Pkg.activate(joinpath(@__DIR__, "..", "examples"))
Pkg.develop(PackageSpec(; path=joinpath(@__DIR__, "..")))
Pkg.instantiate()
Pkg.precompile()
Pkg.activate(@__DIR__)

if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
    ENV["DATADEPS_ALWAYS_ACCEPT"] = true
end

deployconfig = Documenter.auto_detect_deploy_system()
Documenter.post_status(deployconfig; type="pending",
                       repo="github.com/avik-pal/Lux.jl.git")

# Tutorials
get_example_path(p) = joinpath(@__DIR__, "..", "examples", p)
OUTPUT = joinpath(@__DIR__, "src", "examples", "generated")

BEGINNER_TUTORIALS = ["Basics/main.jl", "SimpleRNN/main.jl"]
BEGINNER_TUTORIAL_NAMES = ["Julia & Lux for the Uninitiated", "Training a Simple LSTM"]
INTERMEDIATE_TUTORIALS = ["NeuralODE/main.jl", "BayesianNN/main.jl"]
INTERMEDIATE_TUTORIAL_NAMES = ["MNIST NeuralODE Classification", "Bayesian Neural Network"]
ADVANCED_TUTORIALS = []
ADVANCED_TUTORIAL_NAMES = []
MAPPING = Dict("beginner" => [], "intermediate" => [], "advanced" => [])

for (d, names, paths) in (("beginner", BEGINNER_TUTORIAL_NAMES, BEGINNER_TUTORIALS),
                          ("intermediate", INTERMEDIATE_TUTORIAL_NAMES,
                           INTERMEDIATE_TUTORIALS),
                          ("advanced", ADVANCED_TUTORIAL_NAMES, ADVANCED_TUTORIALS))
    for (n, p) in zip(names, paths)
        Literate.markdown(get_example_path(p), joinpath(OUTPUT, d, dirname(p));
                          documenter=true)
        push!(MAPPING[d],
              n => joinpath("examples/generated", d, dirname(p),
                            splitext(basename(p))[1] * ".md"))
    end
end

display(MAPPING)

makedocs(;
         sitename="Lux",
         authors="Avik Pal et al.",
         clean=true,
         doctest=false,
         modules=[Lux],
         format=Documenter.HTML(;
                                prettyurls=get(ENV, "CI", nothing) == "true",
                                assets=["assets/custom.css"], edit_branch="main",
                                analytics="G-Q8GYTEVTZ2"),
         pages=[
             "Lux: Explicitly Parameterized Neural Networks" => "index.md",
             "Introduction" => [
                 "All about Lux" => "introduction/overview.md",
                 "Ecosystem" => "introduction/ecosystem.md",
             ],
             "Examples" => [
                 "Beginner" => MAPPING["beginner"],
                 "Intermediate" => MAPPING["intermediate"],
                 "Advanced" => MAPPING["advanced"],
                 "Additional Examples" => "examples.md",
             ],
             "API" => [
                 "Layers" => "api/layers.md",
                 "Functional" => "api/functional.md",
                 "Core" => "api/core.md",
                 "Utilities" => "api/utilities.md",
             ],
             "Design Docs" => [
                 "Documentation" => "design/documentation.md",
                 "Recurrent Neural Networks" => "design/recurrent.md",
                 "Add new functionality to Lux" => "design/core.md",
             ],
         ])

deploydocs(; repo="github.com/avik-pal/Lux.jl.git", push_preview=true,
           devbranch="main")

Pkg.activate(@__DIR__)

ext: jl


text:
# Import libraries.
using Turing, Flux, Plots, Random, ReverseDiff

# Hide sampling progress.
Turing.setprogress!(false);

# Use reverse_diff due to the number of parameters in neural networks.
Turing.setadbackend(:reversediff)

# Number of points to generate.
N = 80
M = round(Int, N / 4)
Random.seed!(1234)

# Generate artificial data.
x1s = rand(M) * 4.5;
x2s = rand(M) * 4.5;
xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i in 1:M])
x1s = rand(M) * 4.5;
x2s = rand(M) * 4.5;
append!(xt1s, Array([[x1s[i] - 5; x2s[i] - 5] for i in 1:M]))

x1s = rand(M) * 4.5;
x2s = rand(M) * 4.5;
xt0s = Array([[x1s[i] + 0.5; x2s[i] - 5] for i in 1:M])
x1s = rand(M) * 4.5;
x2s = rand(M) * 4.5;
append!(xt0s, Array([[x1s[i] - 5; x2s[i] + 0.5] for i in 1:M]))

# Store all the data for later.
xs = [xt1s; xt0s]
ts = [ones(2 * M); zeros(2 * M)]

# Plot data points.
function plot_data()
    x1 = map(e -> e[1], xt1s)
    y1 = map(e -> e[2], xt1s)
    x2 = map(e -> e[1], xt0s)
    y2 = map(e -> e[2], xt0s)

    Plots.scatter(x1, y1; color="red", clim=(0, 1))
    return Plots.scatter!(x2, y2; color="blue", clim=(0, 1))
end

plot_data()

# Construct a neural network using Flux
nn_initial = Chain(Dense(2, 3, tanh), Dense(3, 2, tanh), Dense(2, 1, σ))

# Extract weights and a helper function to reconstruct NN from weights
parameters_initial, reconstruct = Flux.destructure(nn_initial)

length(parameters_initial) # number of paraemters in NN

# Create a regularization term and a Gaussian prior variance term.
alpha = 0.09
sig = sqrt(1.0 / alpha)

# Specify the probabilistic model.
@model function bayes_nn(xs, ts, nparameters, reconstruct)
    # Create the weight and bias vector.
    parameters ~ MvNormal(zeros(nparameters), sig .* ones(nparameters))

    # Construct NN from parameters
    nn = reconstruct(parameters)
    # Forward NN to make predictions
    preds = nn(xs)

    # Observe each prediction.
    for i in 1:length(ts)
        ts[i] ~ Bernoulli(preds[i])
    end
end;

# Perform inference.
N = 5000
ch = sample(
    bayes_nn(hcat(xs...), ts, length(parameters_initial), reconstruct), HMC(0.05, 4), N
);

# Extract all weight and bias parameters.
theta = MCMCChains.group(ch, :parameters).value;

# A helper to create NN from weights `theta` and run it through data `x`
nn_forward(x, theta) = reconstruct(theta)(x)

# Plot the data we have.
plot_data()

# Find the index that provided the highest log posterior in the chain.
_, i = findmax(ch[:lp])

# Extract the max row value from i.
i = i.I[1]

# Plot the posterior distribution with a contour plot
x1_range = collect(range(-6; stop=6, length=25))
x2_range = collect(range(-6; stop=6, length=25))
Z = [nn_forward([x1, x2], theta[i, :])[1] for x1 in x1_range, x2 in x2_range]
contour!(x1_range, x2_range, Z)

# Return the average predicted value across
# multiple weights.
function nn_predict(x, theta, num)
    return mean([nn_forward(x, theta[i, :])[1] for i in 1:10:num])
end;

# Plot the average prediction.
plot_data()

n_end = 1500
x1_range = collect(range(-6; stop=6, length=25))
x2_range = collect(range(-6; stop=6, length=25))
Z = [nn_predict([x1, x2], theta, n_end)[1] for x1 in x1_range, x2 in x2_range]
contour!(x1_range, x2_range, Z)

# Number of iterations to plot.
n_end = 500

anim = @gif for i in 1:n_end
    plot_data()
    Z = [nn_forward([x1, x2], theta[i, :])[1] for x1 in x1_range, x2 in x2_range]
    contour!(x1_range, x2_range, Z; title="Iteration $i", clim=(0, 1))
end every 5

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
@doc raw"""
    SpecialEuclidean(n)

Special Euclidean group $\mathrm{SE}(n)$, the group of rigid motions.

``\mathrm{SE}(n)`` is the semidirect product of the [`TranslationGroup`](@ref) on $ℝ^n$ and
[`SpecialOrthogonal(n)`](@ref)

````math
\mathrm{SE}(n) ≐ \mathrm{T}(n) ⋊_θ \mathrm{SO}(n),
````

where $θ$ is the canonical action of $\mathrm{SO}(n)$ on $\mathrm{T}(n)$ by vector rotation.

This constructor is equivalent to calling

```julia
Tn = TranslationGroup(n)
SOn = SpecialOrthogonal(n)
SemidirectProductGroup(Tn, SOn, RotationAction(Tn, SOn))
```

Points on $\mathrm{SE}(n)$ may be represented as points on the underlying product manifold
$\mathrm{T}(n) × \mathrm{SO}(n)$. For group-specific functions, they may also be
represented as affine matrices with size `(n + 1, n + 1)` (see [`affine_matrix`](@ref)), for
which the group operation is [`MultiplicationOperation`](@ref).
"""
const SpecialEuclidean{N} = SemidirectProductGroup{
    ℝ,
    TranslationGroup{Tuple{N},ℝ},
    SpecialOrthogonal{N},
    RotationAction{TranslationGroup{Tuple{N},ℝ},SpecialOrthogonal{N},LeftAction},
}

function SpecialEuclidean(n)
    Tn = TranslationGroup(n)
    SOn = SpecialOrthogonal(n)
    A = RotationAction(Tn, SOn)
    return SemidirectProductGroup(Tn, SOn, A)
end

Base.show(io::IO, ::SpecialEuclidean{n}) where {n} = print(io, "SpecialEuclidean($(n))")

Base.@propagate_inbounds function submanifold_component(
    ::SpecialEuclidean{n},
    p::AbstractMatrix,
    ::Val{1},
) where {n}
    return view(p, 1:n, n + 1)
end
Base.@propagate_inbounds function submanifold_component(
    ::SpecialEuclidean{n},
    p::AbstractMatrix,
    ::Val{2},
) where {n}
    return view(p, 1:n, 1:n)
end

function submanifold_components(G::SpecialEuclidean{n}, p::AbstractMatrix) where {n}
    @assert size(p) == (n + 1, n + 1)
    @inbounds t = submanifold_component(G, p, Val(1))
    @inbounds R = submanifold_component(G, p, Val(2))
    return (t, R)
end

Base.@propagate_inbounds function _padpoint!(
    ::SpecialEuclidean{n},
    q::AbstractMatrix,
) where {n}
    for i in 1:n
        q[n + 1, i] = 0
    end
    q[n + 1, n + 1] = 1
    return q
end

Base.@propagate_inbounds function _padvector!(
    ::SpecialEuclidean{n},
    X::AbstractMatrix,
) where {n}
    for i in 1:(n + 1)
        X[n + 1, i] = 0
    end
    return X
end

@doc raw"""
    affine_matrix(G::SpecialEuclidean, p) -> AbstractMatrix

Represent the point $p ∈ \mathrm{SE}(n)$ as an affine matrix.
For $p = (t, R) ∈ \mathrm{SE}(n)$, where $t ∈ \mathrm{T}(n), R ∈ \mathrm{SO}(n)$, the
affine representation is the $n + 1 × n + 1$ matrix

````math
\begin{pmatrix}
R & t \\
0^\mathrm{T} & 1
\end{pmatrix}.
````

See also [`screw_matrix`](@ref) for matrix representations of the Lie algebra.
"""
function affine_matrix(G::SpecialEuclidean{n}, p) where {n}
    pis = submanifold_components(G, p)
    pmat = allocate_result(G, affine_matrix, pis...)
    map(copyto!, submanifold_components(G, pmat), pis)
    @inbounds _padpoint!(G, pmat)
    return pmat
end
affine_matrix(::SpecialEuclidean{n}, p::AbstractMatrix) where {n} = p
function affine_matrix(::GT, ::Identity{GT}) where {n,GT<:SpecialEuclidean{n}}
    s = maybesize(Size(n, n))
    s isa Size && return SDiagonal{n,Float64}(I)
    return Diagonal{Float64}(I, n)
end

@doc raw"""
    screw_matrix(G::SpecialEuclidean, X) -> AbstractMatrix

Represent the Lie algebra element $X ∈ 𝔰𝔢(n) = T_e \mathrm{SE}(n)$ as a screw matrix.
For $X = (b, Ω) ∈ 𝔰𝔢(n)$, where $Ω ∈ 𝔰𝔬(n) = T_e \mathrm{SO}(n)$, the screw representation is
the $n + 1 × n + 1$ matrix

````math
\begin{pmatrix}
Ω & b \\
0^\mathrm{T} & 0
\end{pmatrix}.
````

See also [`affine_matrix`](@ref) for matrix representations of the Lie group.
"""
function screw_matrix(G::SpecialEuclidean{n}, X) where {n}
    Xis = submanifold_components(G, X)
    Xmat = allocate_result(G, screw_matrix, Xis...)
    map(copyto!, submanifold_components(G, Xmat), Xis)
    @inbounds _padvector!(G, Xmat)
    return Xmat
end
screw_matrix(::SpecialEuclidean{n}, X::AbstractMatrix) where {n} = X

function allocate_result(G::SpecialEuclidean{n}, f::typeof(affine_matrix), p...) where {n}
    return allocate(p[1], Size(n + 1, n + 1))
end
function allocate_result(G::SpecialEuclidean{n}, f::typeof(screw_matrix), X...) where {n}
    return allocate(X[1], Size(n + 1, n + 1))
end

compose(::SpecialEuclidean, p::AbstractMatrix, q::AbstractMatrix) = p * q

function compose!(
    ::SpecialEuclidean,
    x::AbstractMatrix,
    p::AbstractMatrix,
    q::AbstractMatrix,
)
    return mul!(x, p, q)
end

@doc raw"""
    group_exp(G::SpecialEuclidean{n}, X)

Compute the group exponential of $X = (b, Ω) ∈ 𝔰𝔢(n)$, where $b ∈ 𝔱(n)$ and $Ω ∈ 𝔰𝔬(n)$:

````math
\exp X = (t, R),
````

where $t ∈ \mathrm{T}(n)$ and $R = \exp Ω$ is the group exponential on $\mathrm{SO}(n)$.

In the [`screw_matrix`](@ref) representation, the group exponential is the matrix
exponential (see [`group_exp`](@ref)).
"""
group_exp(::SpecialEuclidean, ::Any)

@doc raw"""
    group_exp(G::SpecialEuclidean{2}, X)

Compute the group exponential of $X = (b, Ω) ∈ 𝔰𝔢(2)$, where $b ∈ 𝔱(2)$ and $Ω ∈ 𝔰𝔬(2)$:

````math
\exp X = (t, R) = (U(θ) b, \exp Ω),
````

where $t ∈ \mathrm{T}(2)$, $R = \exp Ω$ is the group exponential on $\mathrm{SO}(2)$,

````math
U(θ) = \frac{\sin θ}{θ} I_2 + \frac{1 - \cos θ}{θ^2} Ω,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_exp(::SpecialEuclidean{2}, ::Any)

@doc raw"""
    group_exp(G::SpecialEuclidean{3}, X)

Compute the group exponential of $X = (b, Ω) ∈ 𝔰𝔢(3)$, where $b ∈ 𝔱(3)$ and $Ω ∈ 𝔰𝔬(3)$:

````math
\exp X = (t, R) = (U(θ) b, \exp Ω),
````

where $t ∈ \mathrm{T}(3)$, $R = \exp Ω$ is the group exponential on $\mathrm{SO}(3)$,

````math
U(θ) = I_3 + \frac{1 - \cos θ}{θ^2} Ω + \frac{θ - \sin θ}{θ^3} Ω^2,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_exp(::SpecialEuclidean{3}, ::Any)

function group_exp!(G::SpecialEuclidean, q, X)
    Xmat = screw_matrix(G, X)
    qmat = exp(Xmat)
    map(copyto!, submanifold_components(G, q), submanifold_components(G, qmat))
    _padpoint!(G, q)
    return q
end
function group_exp!(G::SpecialEuclidean{2}, q, X)
    SO2 = submanifold(G, 2)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(R) == (2, 2)
    @assert size(t) == (2,)
    @assert size(b) == (2,)

    θ = vee(SO2, Identity(SO2, Ω), Ω)[1]
    sinθ, cosθ = sincos(θ)
    if θ ≈ 0
        α = 1 - θ^2 / 6
        β = θ / 2
    else
        α = sinθ / θ
        β = (1 - cosθ) / θ
    end

    @inbounds begin
        R[1] = cosθ
        R[2] = sinθ
        R[3] = -sinθ
        R[4] = cosθ
        t[1] = α * b[1] - β * b[2]
        t[2] = α * b[2] + β * b[1]
        _padpoint!(G, q)
    end
    return q
end
function group_exp!(G::SpecialEuclidean{3}, q, X)
    SO3 = submanifold(G, 2)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(R) == (3, 3)
    @assert size(t) == (3,)

    θ = norm(SO3, Identity(SO3, Ω), Ω) / sqrt(2)
    θ² = θ^2
    if θ ≈ 0
        α = 1 - θ² / 6
        β = θ / 2
        γ = 1 / 6 - θ² / 120
    else
        sinθ, cosθ = sincos(θ)
        α = sinθ / θ
        β = (1 - cosθ) / θ²
        γ = (1 - α) / θ²
    end

    Ω² = Ω^2
    Jₗ = I + β .* Ω .+ γ .* Ω²
    R .= I + α .* Ω .+ β .* Ω²
    copyto!(t, Jₗ * b)
    @inbounds _padpoint!(G, q)
    return q
end

@doc raw"""
    group_log(G::SpecialEuclidean{n}, p) where {n}

Compute the group logarithm of $p = (t, R) ∈ \mathrm{SE}(n)$, where $t ∈ \mathrm{T}(n)$
and $R ∈ \mathrm{SO}(n)$:

````math
\log p = (b, Ω),
````

where $b ∈ 𝔱(n)$ and $Ω = \log R ∈ 𝔰𝔬(n)$ is the group logarithm on $\mathrm{SO}(n)$.

In the [`affine_matrix`](@ref) representation, the group logarithm is the matrix logarithm
(see [`group_log`](@ref)):
"""
group_log(::SpecialEuclidean, ::Any)

@doc raw"""
    group_log(G::SpecialEuclidean{2}, p)

Compute the group logarithm of $p = (t, R) ∈ \mathrm{SE}(2)$, where $t ∈ \mathrm{T}(2)$
and $R ∈ \mathrm{SO}(2)$:

````math
\log p = (b, Ω) = (U(θ)^{-1} t, \log R),
````

where $b ∈ 𝔱(2)$, $Ω = \log R ∈ 𝔰𝔬(2)$ is the group logarithm on $\mathrm{SO}(2)$,

````math
U(θ) = \frac{\sin θ}{θ} I_2 + \frac{1 - \cos θ}{θ^2} Ω,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_log(::SpecialEuclidean{2}, ::Any)

@doc raw"""
    group_log(G::SpecialEuclidean{3}, p)

Compute the group logarithm of $p = (t, R) ∈ \mathrm{SE}(3)$, where $t ∈ \mathrm{T}(3)$
and $R ∈ \mathrm{SO}(3)$:

````math
\log p = (b, Ω) = (U(θ)^{-1} t, \log R),
````

where $b ∈ 𝔱(3)$, $Ω = \log R ∈ 𝔰𝔬(3)$ is the group logarithm on $\mathrm{SO}(3)$,

````math
U(θ) = I_3 + \frac{1 - \cos θ}{θ^2} Ω + \frac{θ - \sin θ}{θ^3} Ω^2,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_log(::SpecialEuclidean{3}, ::Any)

function group_log!(G::SpecialEuclidean, X, q)
    qmat = affine_matrix(G, q)
    Xmat = real(log_safe(qmat))
    map(copyto!, submanifold_components(G, X), submanifold_components(G, Xmat))
    _padvector!(G, X)
    return X
end
function group_log!(G::SpecialEuclidean{2}, X, q)
    SO2 = submanifold(G, 2)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(b) == (2,)

    group_log!(SO2, Ω, R)
    @inbounds θ = Ω[2]
    β = θ / 2
    α = θ ≈ 0 ? 1 - β^2 / 3 : β * cot(β)

    @inbounds begin
        b[1] = α * t[1] + β * t[2]
        b[2] = α * t[2] - β * t[1]
        _padvector!(G, X)
    end
    return X
end
function group_log!(G::SpecialEuclidean{3}, X, q)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(Ω) == (3, 3)
    @assert size(b) == (3,)

    trR = tr(R)
    cosθ = (trR - 1) / 2
    θ = acos(clamp(cosθ, -1, 1))
    θ² = θ^2
    if θ ≈ 0
        α = 1 / 2 + θ² / 12
        β = 1 / 12 + θ² / 720
    else
        sinθ = sin(θ)
        α = θ / sinθ / 2
        β = 1 / θ² - (1 + cosθ) / 2 / θ / sinθ
    end

    Ω .= (R .- transpose(R)) .* α
    Jₗ⁻¹ = I - Ω ./ 2 .+ β .* Ω^2
    mul!(b, Jₗ⁻¹, t)
    @inbounds _padvector!(G, X)
    return X
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
include("../common.jl")

gap_version = v"400.1192.000"
gap_lib_version = v"400.1192.000"
name = "NormalizInterface"
upstream_version = v"1.3.0" # when you increment this, reset offset to v"0.0.0"
offset = v"0.0.1" # increment this when rebuilding with unchanged upstream_version, e.g. gap_version changes
version = offset_version(upstream_version, offset)

# Collection of sources required to build libsingular-julia
sources = [
    ArchiveSource("https://github.com/gap-packages/$(name)/releases/download/v$(upstream_version)/$(name)-$(upstream_version).tar.bz2",
                  "5fd7a1e13bbd48cd119456745f00f8208cd09109c0f6614f0165e1965b92afdd"),
]

# Bash recipe for building across all platforms
script = raw"""
cd NormalizInterface*

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/NormalizInterface.patch

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-gaproot=${prefix}/share/gap --with-normaliz=${prefix}
make -j${nproc}

# copy the loadable module
mkdir -p ${prefix}/lib/gap
cp bin/*/*.so ${prefix}/lib/gap/

install_license LICENSE
"""

name = gap_pkg_name(name)
platforms, dependencies = setup_gap_package(gap_version, gap_lib_version; uses_cxx = true)

push!(dependencies, Dependency("normaliz_jll", compat = "~300.900.100")) # version 3.9.1

# The products that we will ensure are always built
products = [
    FileProduct("lib/gap/NormalizInterface.so", :NormalizInterface),
]

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


text:
using LadderDGA
using SeriesAcceleration
using Test

mP_1 = ModelParameters(1.1, 1.2, 1.3, 1.4)
sP_1 = SimulationParameters(1,2,false,:nothing,:nothing,:nothing,:common,:native,false,LadderDGA.zero_χ_fill,[0,1,2,3],[0,1,2,3], 0.1)
sP_2 = SimulationParameters(1,2,false,:richardson,:coeffs,:nothing,:common,:native,false,LadderDGA.zero_χ_fill,[0,1,2,3],[0,1,2,3], 0.1)

@testset "Config" begin
    include("Config.jl")
end

@testset "Helpers" begin
    include("helpers.jl")
end

@testset "GFFit" begin
    include("GFFit.jl")
end

@testset "ladderDGATools" begin
    include("ladderDGATools.jl")
end

@testset "full run" begin
    #bubble_f, χch_f, χsp_f, trilexch_f, trilexsp_f, Σ_f = loadFortranData(dir)

    #bubble, χch, χsp, χsp_λ, usable_sp, usable_ch,trilexch, trilexsp, Σ_ladder = calculate_Σ_ladder("../config.toml")
    #@test all(bubble .≈ bubble_f)
    #@test all(χsp[usable_sp,:] .≈ χch_f[usable_sp,:])
    #@test all(χch[usable_ch,:] .≈ χch_f[usable_ch,:])
end


text:
struct Cylinder <: Shape
    id::Int64
    radius::Float64
    zmin::Float64
    zmax::Float64
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Cylinder(T::Transformation) = Cylinder(1, 1.0, -1.0, 1.0, T, inv(T))
Cylinder(id::Int64, radius::Real, zmin::Real, zmax::Real) = Cylinder(id, radius, zmin, zmax, Transformation(), Transformation())

can_intersect(C::Cylinder) = true
area(C::Cylinder) = 2π * C.radius * (C.zmax - C.zmin)

function obj_bounds(C::Cylinder)
    pMin = Point3(-C.radius, -C.radius, C.zmin)
    pMax = Point3( C.radius,  C.radius, C.zmax)
    return BoundingBox(pMin, pMax)
end

function world_bounds(C::Cylinder)
    C.obj_to_world(obj_bounds(C))
end

function (T::Transformation)(C::Cylinder)
    T2 = T * C.obj_to_world
    Cylinder(C.id, C.radius, C.zmin, C.zmax, T2, inv(T2))
end

function cylinder_coefs(ray::Ray, cyl::Cylinder)
    A = ray.direction.x^2 + ray.direction.y^2
    B = 2.0 * (ray.direction.x * ray.origin.x + ray.direction.y * ray.origin.y)
    C = ray.origin.x^2 + ray.origin.y^2 - cyl.radius^2
    return A, B, C
end

function shape_intersect(R::Ray, cyl::Cylinder)
    ray = cyl.world_to_obj(R)

    A, B, C = cylinder_coefs(ray, cyl)
    if A ≈ 0.0
        return nothing, NaN, NaN
    end
    hit, tnear, tfar = quadratic(A,B,C)

    if !hit || tnear > ray.tmax || tfar < ray.tmin || (tnear < ray.tmin && tfar > ray.tmax)
        # The ray didn't even hit the infinitely long extension of the cylinder.
        return nothing, NaN, NaN
    end

    P1 = ray(tnear)
    P2 = ray(tfar)

    if cyl.zmin <= P1.z <= cyl.zmax && tnear >= ray.tmin
        P = P1
        t = tnear
    elseif cyl.zmin <= P2.z <= cyl.zmax && tnear <= ray.tmax
        P = P2
        t = tfar
    else
        return nothing, NaN, NaN
    end

    n = normalize(Normal3(P.x, P.y, 0.0))
    s = Vector3(0.0, 0.0, 1.0)

    n = flip_normal(n, ray.direction)

    DG = DifferentialGeometry(cyl.obj_to_world(P), cyl.obj_to_world(n), cyl.obj_to_world(s))

    return DG, t, 5e-4 * t
end

function intersectP(r::Ray, cyl::Cylinder)
    ray = cyl.world_to_obj(r)
    A, B, C = cylinder_coefs(ray, cyl)
    if A ≈ 0.0
        return false
    end
    hit, tnear, tfar = quadratic(A,B,C)

    if !hit || tnear > ray.tmax || tfar < ray.tmin || (tnear < ray.tmin && tfar > ray.tmax)
        return false
    end

    P1 = ray(tnear)
    P2 = ray(tfar)

    return (cyl.zmin <= P1.z <= cyl.zmax && tnear >= ray.tmin) || (cyl.zmin <= P2.z <= cyl.zmax && tnear <= ray.tmax)

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using BinaryBuilder

# Collection of sources required to build MbedTLS
# mbedtls release: 2.16.0
sources = [
    "https://github.com/ARMmbed/mbedtls.git" =>
    "fb1972db23da39bd11d4f9c9ea6266eee665605b",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mbedtls
mkdir -p $prefix/lib
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DUSE_SHARED_MBEDTLS_LIBRARY=On
make -j${nproc} && make install
"""

# 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(prefix) = [
    LibraryProduct(prefix, "libmbedcrypto", :libmbedcrypto),
    LibraryProduct(prefix, "libmbedx509", :libmbedx509),
    LibraryProduct(prefix, "libmbedtls", :libmbedtls)
]

# 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, src_name, src_version, sources, script, platforms, products, dependencies; kwargs...)
build_tarballs(ARGS, "MbedTLS", v"2.16.0", sources, script, platforms, products, dependencies)


text:
function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
              Δt0::TimeInterval=zeroI, kwargs...) where {N}

    @unpack δ, approx_model, static, dim, recursive = alg

    if haskey(kwargs, :NSTEPS)
        NSTEPS = kwargs[:NSTEPS]
        T = NSTEPS * δ
    else
        # get time horizon from the time span imposing that it is of the form (0, T)
        T = _get_T(tspan, check_zero=true, check_positive=true)
        NSTEPS = ceil(Int, T / δ)
    end

    # normalize system to canonical form
    ivp_norm = _normalize(ivp)

    # homogeneize the initial-value problem
    if haskey(kwargs, :homogeneize) && kwargs[:homogeneize] == true
        ivp_norm = homogeneize(ivp_norm)
    end

    # discretize system
    ivp_discr = discretize(ivp_norm, δ, approx_model)
    Φ = state_matrix(ivp_discr)
    Ω0 = initial_state(ivp_discr)
    X = stateset(ivp_discr)

    # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
    got_homogeneous = !hasinput(ivp_discr)

    # this algorithm requires Ω0 to be hyperrectangle
    Ω0 = _overapproximate(Ω0, Hyperrectangle)

    # reconvert the set of initial states and state matrix, if needed
    #static = haskey(kwargs, :static) ? kwargs[:static] : alg.stati

    Ω0 = _reconvert(Ω0, static, dim)
    Φ = _reconvert(Φ, static, dim)

    # preallocate output flowpipe
    #N = eltype(Ω0)
    HT = typeof(Ω0)
    F = Vector{ReachSet{N, HT}}(undef, NSTEPS)

    if got_homogeneous
        reach_homog_BOX!(F, Ω0, Φ, NSTEPS, δ, X, recursive, Δt0)
    else
        U = inputset(ivp_discr)
        @assert isa(U, LazySet) "expcted input of type `<:LazySet`, but got $(typeof(U))"
        # TODO: can we use support function evaluations for the input set?
        U = overapproximate(U, Hyperrectangle)
        reach_inhomog_BOX!(F, Ω0, Φ, NSTEPS, δ, X, U, recursive, Δt0)
    end

    return Flowpipe(F)
end


text:
using Parameters

"""
Global configurable options for Metatheory.
"""
@with_kw mutable struct Options
    """
    Print or not information such as saturation reports
    """
    verbose::Bool = false
    """
    Print iteration numbers in equality saturation
    """
    printiter::Bool = false
end

options = Options()


text:
module BasesAndSamples

using Combinatorics

export basis_monomial, basis_laguerre, basis_jacobi, basis_chebyshev, basis_gegenbauer
export sample_points_simplex, sample_points_padua, sample_points_rescaled_laguerre, sample_points_chebyshev, sample_points_chebyshev_mod

"""
    basis_monomial(d::Int, x...)

Generate the monomial basis in variables x... up to degree d (inclusive).
"""
function basis_monomial(d::Int, x...)
    n = length(x)
    if n > 1
        q = Vector{typeof(prod(one.(x)))}(undef, binomial(n+d,d))
    elseif n == 1
        q = Vector{typeof(one(x...))}(undef, binomial(n+d,d))
    end
    idx = 1
    for k=0:d
        for exponent in multiexponents(n,k)
            q[idx] = prod([x[i]^exponent[i] for i=1:length(exponent)])
            idx+=1
        end
    end
    return q
end

"""
    basis_laguerre(d::Integer, alpha, x)

Generate the Laguerre polynomials with parameter alpha up to degree d (inclusive). (https://en.wikipedia.org/wiki/Laguerre_polynomials)
"""
function basis_laguerre(d::Integer, alpha, x)
    v = Vector{typeof(one(alpha) * one(x))}(undef, 1 + d)
    v[1] = one(x)
    d == 0 && return v
    v[2] = 1 + alpha - x
    d == 1 && return v
    for l = 2:d
        v[l+1] = 1 // big(l) * ((2l - 1 + alpha - x) * v[l] - (l + alpha - 1) * v[l-1])
    end
    return v
end

"""
    basis_jacobi(d::Integer, alpha, beta, x, normalized = true)

Generate the Jacobi polynomials with parameters alpha and beta up to degree d (inclusive). (https://en.wikipedia.org/wiki/Jacobi_polynomials)
"""
function basis_jacobi(d::Integer, alpha, beta, x, normalized = true)
    q = Vector{typeof(one(alpha) * one(x))}(undef, d + 1)
    q[1] = one(x)
    d == 0 && return q
    q[2] = x # normalized
    if !normalized
        q[2] *= (alpha + 1)
    end
    d == 1 && return q
    for k = 2:d
        # what if alpha+beta = -n for some integer 1 <= n < 2d ?
        q[k+1] =
            (2 * k + alpha + beta - 1) /
            BigFloat(2k * (k + alpha + beta) * (2k + alpha + beta - 2)) *
            ((2 * k + alpha + beta) * (2k + alpha + beta - 2) * x + beta^2 - alpha^2) *
            q[k] + -2 * (k + alpha - 1) * (k + beta - 1) * (2 * k + alpha + beta) * q[k-1]
    end
    return q
end

"""
    basis_chebyshev(d::Int,x)

Generate a basis of chebyshev polynomials up to degree d (inclusive). (https://en.wikipedia.org/wiki/Chebyshev_polynomials)
"""
function basis_chebyshev(d::Int,x)
    v = Vector{typeof(one(x))}(undef, 1 + d)
    v[1] = one(x)
    d == 0 && return v
    v[2] = x
    d == 1 && return v
    for l = 2:d
        v[l+1] = 2x*v[l] - v[l-1]
    end
    return v
end

"""
    basis_gegenbauer(d, n, x)

Basis for the Gegenbauer polynomials in dimension n up to degree d.
 This is the Gegenbauer polynomial with parameter lambda = n/2-1,
 or the Jacobi polynomial with alpha = beta = (n-3)/2.
 Normalized to evaluate to 1 at 1.
 Taken from arxiv/2001.00256, ancillary files, SemidefiniteProgramming.jl.
"""
function basis_gegenbauer(d, n, x)
    v = Vector{typeof(one(x))}(undef, 1 + d)
    v[1] = one(x)
    d == 0 && return v
    v[2] = x
    d == 1 && return v
    for l = 2:d
        v[l+1] = (2l + n - 4) // (l + n - 3) * x * v[l] - (l - 1) // (l + n - 3) * v[l-1]
    end
    return v
end

"""
    sample_points_simplex(n, d)

Generate the rational sample points in the unit simplex with denominator d.
"""
function sample_points_simplex(n, d)
    #rational points in the unit simplex with denominator d
    #probably not very efficient, but I dont know how to do it better for general n.
    x = [zeros(BigFloat, n) for i = 1:binomial(n + d, d)] #need n+d choose d points for a unisolvent set, if symmetry is not used.
    idx = 1
    for I in CartesianIndices(ntuple(k -> 0:d, n)) #all tuples with elements in 0:d of length n
        if sum(Tuple(I)) <= d #in unit simplex
            x[idx] = [i / BigFloat(d) for i in Tuple(I)]
            idx += 1
        end
    end
    return x
end

"""
    sample_points_padua(d)

Generate the Padua points for degree d. (https://en.wikipedia.org/wiki/Padua_points)
"""
function sample_points_padua(d)
    #padua points:
    z = [Array{BigFloat}(undef, 2) for i = 1:binomial(2 + d, d)]
    z_idx = 1
    for j = 0:d
        delta_j = j % 2 == d % 2 == 1 ? 1 : 0
        mu_j = cospi(j / d)
        for k = 1:(div(d, 2)+1+delta_j)
            eta_k = j % 2 == 1 ? cospi((2 * k - 2) / (d + 1)) : cospi((2 * k - 1) / (d + 1))
            z[z_idx] = [mu_j, eta_k]
            z_idx += 1
        end
    end
    return z
end

"""
    sample_points_rescaled_laguerre(d)

Generate 'rescaled laguerre' points, as in SDPB.
"""
function sample_points_rescaled_laguerre(d)
    #as done in sdpb: ('rescaled Laguerre')
    # x[k] = sqrt(pi)/(64*(d+1)*log( 3- 2*sqrt(2))) * (-1+4*k)^2, with k=0:d
    constant = -sqrt(BigFloat(pi)) / (64 * (d + 1) * log(3 - 2 * sqrt(BigFloat(2))))
    x = zeros(BigFloat, d + 1)
    for k = 0:d
        x[k+1] = constant * (-1 + 4 * k)^2
    end
    return [[i] for i in x]
end

"""
    sample_points_chebyshev(d, a = -1, b = 1)

Generate the d+1 chebyshev points in [a,b]. (https://en.wikipedia.org/wiki/Chebyshev_nodes)
"""
function sample_points_chebyshev(d, a = -1, b = 1)
    #roots of chebyshev polynomials of the first kind, unisolvent for polynomials up to degree d
    points = [
        (a + b) / BigFloat(2) +
        (b - a) / BigFloat(2) * cos((2k - 1) / BigFloat(2(d + 1)) * BigFloat(pi)) for
        k = 1:d+1
    ]
    return [[i] for i in points]
end

"""
    sample_points_chebyshev_mod(d, a = -1, b = 1)

Generate the d+1 modified chebyshev points in [a,b], the chebyshev points divided by cos(pi/(2(d+1))).
"""
function sample_points_chebyshev_mod(d, a = -1, b = 1)
    #roots of chebyshev polynomials of the first kind, divided by cos(pi/2(d+1)) to get a lower lebesgue constant
    points = [
        (a + b) / BigFloat(2) +
        (b - a) / BigFloat(2) * cos((2k - 1) / BigFloat(2(d + 1)) * BigFloat(pi)) /
        cos(BigFloat(pi) / 2(d + 1)) for k = 1:d+1
    ]
    return [[i] for i in points]
end

end # end module


text:
struct Bar
    baz
    qux::Float64
end


text:
using EllZeroTrendFiltering

N = 10
t = [0.0:N-1...]
I_g = rand(N)
I_g2 = rand(N)
I_tg = rand(N)

ℓ = EllZeroTrendFiltering.TransitionCostContinuous{Float64}(t, I_g, I_g2, I_tg)

#ℓ = compute_transition_costs(x -> exp(2^x), range(0, stop=2, length=N))

ℓ = compute_transition_costs(x -> exp(x^2), range(0, stop=1, length=20))

V_N = QuadraticPolynomial(0.0, 0.0, 0.0)

M = 10

@time Λ =  pwq_dp_constrained(ℓ, V_N, M)

unique_polys(Λ) = length(unique([λ.p for λ in Λ]))

sizes = [(isassigned(Λ,i,j) ? unique_polys(Λ[i,j]) : 0) for i = 1:size(Λ,1), j = 1:size(Λ,2)]
sizes = [(isassigned(Λ,i,j) ? length(Λ[i,j]) : 0) for i = 1:size(Λ,1), j = 1:size(Λ,2)]

# Check compared to bound R <= N
nom_sizes = fill(1, M, 1)*(20-1:-1:0)'

println(maximum(sizes - nom_sizes))


text:
using Documenter, Statistics

# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

makedocs(
    modules = [Statistics],
    sitename = "Statistics",
    pages = Any[
        "Statistics" => "index.md"
        ]
    )

deploydocs(repo = "github.com/JuliaStats/Statistics.jl.git")


text:
using PQEcolaPoint

# init some functions -
function _fitness(world::CRRJITContractPremiumLatticeModel, test_point::PQContractPremiumLatticePoint)::Float64

    # how would we impl constraints?
    p_value = premium(world, test_point)

    # fitness is premium -
    return p_value
end

function _neighborhood(point::PQContractPremiumLatticePoint)::Array{PQContractPremiumLatticePoint,1}

    # initialize -
    N_array = Array{PQContractPremiumLatticePoint,1}()
    D_max = 301
    D_min = 2
    I_max = 16
    I_min = 2

    if (point.d < (D_max - 1) && point.d >= D_min && point.i < (I_max - 1) && point.i >= I_min)

        push!(N_array, build(PQContractPremiumLatticePoint; i = (point.i + 1), j = point.j, s = point.s, d = point.d))
        push!(N_array, build(PQContractPremiumLatticePoint; i = (point.i - 1), j = point.j, s = point.s, d = point.d))
        push!(N_array, build(PQContractPremiumLatticePoint; i = point.i, j = point.j, s = point.s, d = (point.d + 1)))
        push!(N_array, build(PQContractPremiumLatticePoint; i = point.i, j = point.j, s = point.s, d = (point.d - 1)))
        push!(N_array, build(PQContractPremiumLatticePoint; i = (point.i + 1), j = point.j, s = point.s, d = (point.d + 1)))
        push!(N_array, build(PQContractPremiumLatticePoint; i = (point.i + 1), j = point.j, s = point.s, d = (point.d - 1)))
        push!(N_array, build(PQContractPremiumLatticePoint; i = (point.i - 1), j = point.j, s = point.s, d = (point.d + 1)))
        push!(N_array, build(PQContractPremiumLatticePoint; i = (point.i - 1), j = point.j, s = point.s, d = (point.d - 1)))
    else
        push!(N_array, build(PQContractPremiumLatticePoint; i = I_min, j = point.j, s = point.s, d = D_min))
    end

    @show N_array

    # return -
    return N_array
end

# setup DTE -
T = 300.0
DTE = range(T, stop=0.0,step=-1.0) |> collect

# setup underlying range -
Sₒ = 76.04
IV = 0.4456
σ = Sₒ*(IV)*(√(T/365))
L = round(Sₒ - 1.96*σ)  # 95%
U = round(Sₒ + 1.96*σ)  # 95%
S = range(max.(L, 0.0),stop=U, step=0.5) |> collect

# setup IV array -
IV_array = range(0.1, stop = 1.2, step = 0.05) |> collect

# setup strike coordinates -
K = range(50.0, stop = 125.0, step = 5.0) |> collect

# build JIT lattice model -
JIT_lattice = build(CRRJITContractPremiumLatticeModel, PutContractModel, S, K, IV_array, DTE; number_of_levels = 80)

# build a point -
start_point = encode(JIT_lattice; S = 73.0, DTE = 35.0, IV = 0.4498, K = 75.0) # we are underwater here -

# try the tabu -
(F,B) = tabu(JIT_lattice, start_point; fitness = _fitness, neighborhood = _neighborhood)


text:
reciprocal_distretization(x, L) = reciprocal_distretization(lattice(x), L)
function reciprocal_distretization(lattice::AbstractLattice, L)
    r1, r2 = reciprocal_vectors(lattice, L)
    reciprocal_distretization(r1, r2, L)
end
function reciprocal_distretization(v1, v2, L)
    pgram = map(((i, j) for i in 0:L-1 for j in 0:L-1)) do (i, j)
        out = [0., 0.]
        temp = [0., 0.]
        min = Inf
        for n in -1:1, m in -1:1
            temp .= v1 * (n + i/L) .+ v2 * (m + j/L)
            if norm(temp) < min
                min = norm(temp)
                out .= temp
            end
        end
        out
    end
    return pgram
end

structure_factor(args...) = fourier_transform(args...)
function fourier_transform(mc::DQMC, key::Symbol, args...)
    fourier_transform(mc, mc[key], args...)
end
function fourier_transform(mc::DQMC, m::DQMCMeasurement, L::Integer, args...)
    qs = reciprocal_distretization(mc, L)
    return qs, fourier_transform(mc, qs, m, args...)
end
function fourier_transform(mc::DQMC, qs::Vector, m::DQMCMeasurement, args...)
    dirs = directions(mc)
    values = mean(m)
    fourier_transform(qs, dirs, values, args...)
end

function fourier_transform(qs::Vector, dirs::Vector, values::Vector)
    @boundscheck length(dirs) == length(qs) == length(values)
    map(qs) do q
        sum(cis(dot(q, v)) * o for (v, o) in zip(dirs, values))
    end
end

function fourier_transform(qs::Vector, dirs::Vector, values::Matrix, weights::Vector)
    @boundscheck begin
        length(weights) ≤ size(values, 2) &&
        length(dirs) == length(qs) &&
        length(dirs) == size(values, 1)
    end
    map(qs) do q
        out = zero(ComplexF64)
        for i in eachindex(dirs)
            temp = zero(ComplexF64)
            for j in eachindex(weights)
                temp += weights[j] * values[i, j]
            end
            out += cis(dot(q, dirs[i])) * temp
        end
        out
    end
end

function fourier_transform(qs::Vector, dirs::Vector, values::Array{T, 3}, weights::Vector) where {T}
    @boundscheck begin
        length(weights) ≤ size(values, 2) &&
        length(dirs) == length(qs) &&
        length(dirs) == size(values, 1)
    end
    map(qs) do q
        out = zero(ComplexF64)
        @turbo for i in eachindex(dirs)
            temp = zero(ComplexF64)
            for j in eachindex(weights), k in eachindex(weights)
                temp += weights[j] * weights[k] * values[i, j, k]
            end
            out += cis(dot(q, dirs[i])) * temp
        end
        out
    end
end

uniform_fourier(mc::DQMC, key::Symbol) = uniform_fourier(mc[key])
uniform_fourier(m::DQMCMeasurement) = sum(mean(m))

apply_symmetry(mc::DQMC, key::Symbol, weights=(1)) = apply_symmetry(mc[key], weights)
apply_symmetry(m::DQMCMeasurement, weights=(1)) = apply_symmetry(mean(m), weights)
function apply_symmetry(data::AbstractArray{T, 3}, weights=(1)) where {T}
    out = zeros(T, size(data, 1))
    @turbo for i in eachindex(out)
        for j in eachindex(weights), k in eachindex(weights)
            out[i] += weights[j] * weights[k] * data[i,j,k]
        end
    end
    out
end
function apply_symmetry(data::AbstractArray{T, 1}, weights=(1)) where {T}
    out = zero(T)
    @turbo for i in eachindex(weights)
        out += weights[i] * data[i]
    end
    out
end
function apply_symmetry(data::AbstractArray{T, 2}, weights=(1)) where {T}
    out = zeros(T, size(data, 1))
    @turbo for i in eachindex(out)
        for j in eachindex(weights)
            out[i] += weights[j] * data[i,j]
        end
    end
    out
end

################################################################################
### Superfluid Stiffness
################################################################################

function superfluid_stiffness(
        mc::DQMC, G::DQMCMeasurement, ccs::DQMCMeasurement;
        shift_dir = [1., 0.], calculate_error = false
    )
    # find all hopping directions (skipping on-site)
    valid_directions = hopping_directions(mc, mc.model)

    # reduce to directions that have a positive component in shift_dir
    dirs = directions(lattice(mc))
    filter!(idx -> dot(dirs[idx], shift_dir) > 0, valid_directions)

    # Reduce to valid (ccs_idx, dir_idx) pairs (we may have more or less in
    # the measurement)
    ind_dir = Pair{Int, Int}[]
    for (i, j) in enumerate(ccs.lattice_iterator.directions)
        if j in valid_directions
            push!(ind_dir, i => j)
        end
    end
    if length(ind_dir) != length(valid_directions)
        @warn("Missing directions in Superfluid Stiffness")
    end

    Kx  = dia_K_x(mc, mean(G), last.(ind_dir))
    Λxx = para_ccc(mc, mean(ccs), ind_dir)

    if calculate_error
        @warn "Errors not implemented for sfs yet"
        return 0.25 * (-Kx - Λxx), 0.0
    else
        return 0.25 * (-Kx - Λxx)
    end
end

function dia_K_x(mc, G, idxs)
    T = Matrix(MonteCarlo.hopping_matrix(mc, mc.model))

    push!(mc.lattice_iterator_cache, MonteCarlo.Dir2SrcTrg(), lattice(mc))
    dir2srctrg = mc.lattice_iterator_cache[MonteCarlo.Dir2SrcTrg()]
    N = length(lattice(mc))

    Kx = ComplexF64(0)
    f = 1.0 + (model(mc).U > 0)
    for shift in 0:N:size(G, 1)-1
        for i in idxs
            for (src, trg) in dir2srctrg[i]
                # c_j^† c_i = δ_ij - G[i, j], but δ always 0 cause no onsite
                Kx -= f * T[trg+shift, src+shift] * G[src+shift, trg+shift]
                # reverse directions cause we filter those out beforehand
                Kx -= f * T[src+shift, trg+shift] * G[trg+shift, src+shift]
            end
        end
    end
    Kx /= N
end

# uses directions with filter > 0
function para_ccc(mc, ccs, ind_dir)
    #dq = [2pi / size(lattice(mc))[2], 0.0]
    dq = [0.0, 0.0]
    dirs = directions(lattice(mc))
    Λxx = ComplexF64(0)

    for (i, dir) in enumerate(dirs)
        for (j, jdir) in ind_dir, (k, kdir) in ind_dir
            # result += ccs[i, j] * (cis(dot(dir, long)) - cis(dot(dir, trans)))
            Λxx += ccs[i, j, k] * cis(-dot(dir + 0.5(dirs[jdir] .- dirs[kdir]), dq))
        end
    end

    Λxx
end


text:
### Congestion management with renewable energy sources and storages
function CongestionManagement(hours,
							  nodes::Nodes,
							  lines::Lines,
							  powerplants::PowerPlants,
							  renewables::Renewables,
							  storages::Storages,
							  price,
							  P_opt,
							  P_R_opt,
							  P_S_opt,
							  D_S_opt,
						      )

	h,b = calc_h_b(lines, nodes)

	T = hours
	G = powerplants.unit
	R = renewables.unit
	S = storages.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	L = lines.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, ΔP_R_up[T, R] == 0)
	@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_flow[T, L])
	@variable(CM_mod, P_load_lost[T, J] >= 0)
	@variable(CM_mod, P_gen_lost[T, J] >= 0)

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	# Constraints
	# Market clearing/power balance
	# @constraint(CM_mod, MarketClearing[t = T],
	#     sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) for g in G) +
	# 	sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) for r in R) +
	# 	sum(P_S_opt[t, s]/Sbase for s in S) -
	# 	sum(D_S_opt[t, s]/Sbase for s in S)
	# 		== nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g][t]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g][t]/Sbase);

	# Upper generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitUp[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
			renewables.infeed[r][t]/Sbase);

	# Lower generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitDown[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
			0);

	# # Ramp-up limit
	# @constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
	# 	(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
	# 		powerplants.rup[g]/Sbase)
	#
	# # Ramp-down limit
	# @constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
	# 	-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
	# 		powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) +
		sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r])
			for r in getKeyVector(renewables.node, j)) +
		sum((P_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) -
		sum((D_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) -
		P_gen_lost[t, j] +
		sum(b[j,jj] * Θ[t, jj] for jj in J) ==
		nodes.load[j][t]/Sbase - nodes.exchange[j][t]/Sbase - P_load_lost[t, j]);

		@constraint(CM_mod, LinePowerFlowMax[l = L, t = T],
			P_flow[t, l] <= lines.pmax[l]/Sbase);

		@constraint(CM_mod, LinePowerFlowMin[l = L, t = T],
			P_flow[t, l] >= -lines.pmax[l]/Sbase);

		@constraint(CM_mod, LineFlow[l = L, t = T],
			P_flow[t, l] == sum(h[l,j] * Θ[t,j] for j in J));

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T) +
		sum(0 * (ΔP_R_up[t, r]) * Sbase +
			(price[t] - 0) * ΔP_R_dn[t, r] *
				Sbase for r in R, t in T) +
		1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
		1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
	ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	price = JuMP.dual.(PowerInjectionBal)
	P_gen_lost_opt = JuMP.value.(P_gen_lost) * Sbase
	P_load_lost_opt = JuMP.value.(P_load_lost) * Sbase

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
		   ΔP_R_up_opt,
		   ΔP_R_dn_opt,
	       Θ_opt,
	       P_flow_opt,
		   price,
		   P_gen_lost_opt,
		   P_load_lost_opt
	       )
end

function CongestionManagementPHS(hours,
								 nodes::Nodes,
								 lines::Lines,
								 powerplants::PowerPlants,
								 renewables::Renewables,
								 storages::Storages,
								 price,
								 P_opt,
								 P_R_opt,
								 P_S_opt,
								 D_S_opt,
								 L_S_opt,
							     )

	h,b = calc_h_b(lines, nodes)

	T = hours
	G = powerplants.unit
	R = renewables.unit
	S = storages.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	L = lines.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, ΔP_R_up[T, R] == 0)
	@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_flow[T, L])
	@variable(CM_mod, P_load_lost[T, J] >= 0)
	@variable(CM_mod, P_gen_lost[T, J] >= 0)

	# Storage voriables
	@variable(CM_mod, P_S_up[T, S] >= 0)	# Power generation from storage unit
	# @variable(ED_mod, D_S[T, S] >= 0)	# Electricity demand from storage unit
	@variable(CM_mod, L_S[T, S] >= 0)

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	# Fix storage level in period 0 to zero
	if hours[end] <= 24
		for s in S
			fix(L_S[T[1], s], 0; force = true)
		end
	else
		for s in S
			fix(L_S[T[1], s], L_S_opt[(T[1]), s]; force = true)
		end
	end

	# # Fix storage parameters in first hour to the last hour value of the
	# # previous day slice.
	# for s in S
	# 	fix(L_S[T[1], s], L_S_opt[(T[end]), s]/Sbase; force = true)
	# end

	# Constraints
	# Market clearing/power balance
	# @constraint(CM_mod, MarketClearing[t = T],
	#     sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) for g in G) +
	# 	sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) for r in R) +
	# 	sum(P_S_opt[t, s]/Sbase for s in S) -
	# 	sum(D_S_opt[t, s]/Sbase for s in S)
	# 		== nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g][t]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g][t]/Sbase);

	# Upper generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitUp[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
			renewables.infeed[r][t]/Sbase);

	# Lower generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitDown[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
			0);

	# # Ramp-up limit
	# @constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
	# 	(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
	# 		powerplants.rup[g]/Sbase)
	#
	# # Ramp-down limit
	# @constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
	# 	-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
	# 		powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) +
		sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r])
			for r in getKeyVector(renewables.node, j)) +
		sum(((P_S_opt[t, s]/Sbase) + P_S_up[t, s]) for s in getKeyVector(storages.node, j)) -
		sum((D_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) -
		P_gen_lost[t, j] +
		sum(b[j,jj] * Θ[t, jj] for jj in J) ==
		nodes.load[j][t]/Sbase - nodes.exchange[j][t]/Sbase - P_load_lost[t, j]);

	@constraint(CM_mod, LinePowerFlowMax[l = L, t = T],
		P_flow[t, l] <= lines.pmax[l]/Sbase);

	@constraint(CM_mod, LinePowerFlowMin[l = L, t = T],
		P_flow[t, l] >= -lines.pmax[l]/Sbase);

	@constraint(CM_mod, LineFlow[l = L, t = T],
		P_flow[t, l] == sum(h[l,j] * Θ[t,j] for j in J));

	@constraint(CM_mod, LostLoadBound[j = J, t = T],
		P_load_lost[t, j] <= nodes.load[j][t]/Sbase);

	### PHS constraints
	@constraint(CM_mod, StoragePowerOutput[s = S, t = T],
 		P_S_opt[t, s]/Sbase + P_S_up[t, s] <= storages.power[s]/Sbase);

	# Maximum storage power input
	# @constraint(CM_mod, StoragePowerInput[s = S, t = T],
	# 	D_S[t, s] <= storages.power[s]);

	# Generation CM_mod determined by remaining storage level
	@constraint(CM_mod, StorageLevelGen[s = S, t = T],
		P_S_opt[t, s]/Sbase + P_S_up[t, s] <= L_S[t, s]);

	# # Storace capacity limit
	@constraint(CM_mod, StorageLevelCap[s = S, t = T],
		L_S[t, s] <= storages.storage[s]/Sbase);

	# Storage
	@constraint(CM_mod, Storage[s = S, t = T; t > T[1]],
		L_S[t-1, s] - (P_S_opt[t-1, s]/Sbase + P_S_up[t-1, s]) + storages.efficiency[s] * D_S_opt[t-1, s]/Sbase == L_S[t, s]);

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T) +
		sum(0 * (ΔP_R_up[t, r]) * Sbase +
			(price[t] - 0) * ΔP_R_dn[t, r] *
				Sbase for r in R, t in T) +
		sum(5 * P_S_up[t, s] * Sbase for s in S, t in T) +
		1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
		1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
	ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	price = JuMP.dual.(PowerInjectionBal)
	P_gen_lost_opt = JuMP.value.(P_gen_lost) * Sbase
	P_load_lost_opt = JuMP.value.(P_load_lost) * Sbase

	P_S_up_opt = JuMP.value.(P_S_up) * Sbase
	L_S_opt_after_CM = JuMP.value.(L_S) * Sbase

	P_S_opt_after_CM = JuMP.Containers.DenseAxisArray(P_S_opt.data .+ P_S_up_opt.data, P_S_up_opt.axes[1], P_S_up_opt.axes[2])

	# if length(P_S_opt.axes[1]) <= 24
	# 	P_S_opt_after_CM = JuMP.Containers.DenseAxisArray(P_S_opt.data .+ P_S_up_opt.data, P_S_up_opt.axes[1], P_S_up_opt.axes[2])
	# 	# L_S_opt_after_CM = JuMP.Containers.DenseAxisArray(L_S_opt.data .- P_S_up_opt.data, P_S_up_opt.axes[1], P_S_up_opt.axes[2])
	# else
	# 	t_observed = 1:24
	# 	P_S_opt_after_CM = JuMP.Containers.DenseAxisArray(P_S_opt.data[t_observed, :] .+ P_S_up_opt.data, P_S_up_opt.axes[1], P_S_up_opt.axes[2])
	# 	# L_S_opt_after_CM = JuMP.Containers.DenseAxisArray(L_S_opt.data[t_observed, :] .- P_S_up_opt.data, P_S_up_opt.axes[1], P_S_up_opt.axes[2])
	# end

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
		   ΔP_R_up_opt,
		   ΔP_R_dn_opt,
	       Θ_opt,
	       P_flow_opt,
		   price,
		   P_gen_lost_opt,
		   P_load_lost_opt,
		   P_S_opt_after_CM,
		   D_S_opt,
		   L_S_opt_after_CM
	       )
end

### Congestion management with renewable energy sources and storages
function CongestionManagementPHS80(hours,
							  nodes::Nodes,
							  lines::Lines,
							  powerplants::PowerPlants,
							  renewables::Renewables,
							  storages::Storages,
							  price,
							  P_opt,
							  P_R_opt,
							  P_S_opt,
							  D_S_opt,
							  L_S_opt,
						      )

	h,b = calc_h_b(lines, nodes)

	T = hours
	G = powerplants.unit
	R = renewables.unit
	S = storages.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	L = lines.id

	# Initialise JuMP model: Congestion management
	CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(CM_mod, ΔP_up[T, G] >= 0)
	@variable(CM_mod, ΔP_dn[T, G] >= 0)
	@variable(CM_mod, ΔP_R_up[T, R] == 0)
	@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
	@variable(CM_mod, Θ[T, J])
	@variable(CM_mod, P_flow[T, L])
	@variable(CM_mod, P_load_lost[T, J] >= 0)
	@variable(CM_mod, P_gen_lost[T, J] >= 0)

	# Storage voriables
	@variable(CM_mod, P_S_up[T, S] >= 0)	# Power generation from storage unit
	# @variable(ED_mod, D_S[T, S] >= 0)	# Electricity demand from storage unit
	@variable(CM_mod, L_S[T, S] >= 0)

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	for s in S
		fix(L_S[T[1], s], 0.2*storages.storage[s]/Sbase; force = true)
	end

	# Constraints
	# Market clearing/power balance
	# @constraint(CM_mod, MarketClearing[t = T],
	#     sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) for g in G) +
	# 	sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) for r in R) +
	# 	sum(P_S_opt[t, s]/Sbase for s in S) -
	# 	sum(D_S_opt[t, s]/Sbase for s in S)
	# 		== nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(CM_mod, GenerationLimitUp[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
			powerplants.pmax[g][t]/Sbase);

	# Lower generation limit
	@constraint(CM_mod, GenerationLimitDown[g = G, t = T],
	    (P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
			powerplants.pmin[g][t]/Sbase);

	# Upper generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitUp[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
			renewables.infeed[r][t]/Sbase);

	# Lower generation limit for renewables
	@constraint(CM_mod, ResGenerationLimitDown[r = R, t = T],
	    (P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
			0);

	# # Ramp-up limit
	# @constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
	# 	(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
	# 		powerplants.rup[g]/Sbase)
	#
	# # Ramp-down limit
	# @constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
	# 	-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
	# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
	# 		powerplants.rdn[g]/Sbase)

	### DC power flow constraints
	# Power injection balance
	@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
		sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g])
			for g in getKeyVector(powerplants.node, j)) +
		sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r])
			for r in getKeyVector(renewables.node, j)) +
		sum((P_S_opt[t, s]/Sbase + P_S_up[t, s]) for s in getKeyVector(storages.node, j)) -
		sum((D_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) -
		P_gen_lost[t, j] +
		sum(b[j,jj] * Θ[t, jj] for jj in J) ==
		nodes.load[j][t]/Sbase - nodes.exchange[j][t]/Sbase - P_load_lost[t, j]);

	@constraint(CM_mod, LostLoadBound[j = J, t = T],
		P_load_lost[t, j] <= nodes.load[j][t]/Sbase);

	@constraint(CM_mod, LinePowerFlowMax[l = L, t = T],
		P_flow[t, l] <= lines.pmax[l]/Sbase);

	@constraint(CM_mod, LinePowerFlowMin[l = L, t = T],
		P_flow[t, l] >= -lines.pmax[l]/Sbase);

	@constraint(CM_mod, LineFlow[l = L, t = T],
		P_flow[t, l] == sum(h[l,j] * Θ[t,j] for j in J));

	### PHS constraints
	@constraint(CM_mod, StoragePowerOutput[s = S, t = T],
		P_S_opt[t, s]/Sbase + P_S_up[t, s] <= storages.power[s]/Sbase);

	# Maximum storage power input
	# @constraint(CM_mod, StoragePowerInput[s = S, t = T],
	# 	D_S[t, s] <= storages.power[s]);

	# Generation CM_mod determined by remaining storage level
	@constraint(CM_mod, StorageLevelGen[s = S, t = T],
		P_S_up[t, s] <= L_S[t, s]);

	# Storage
	@constraint(CM_mod, Storage[s = S, t = T; t > T[1]],
		L_S[t-1, s] - P_S_up[t-1, s] == L_S[t, s]);

	# Objective function
	@objective(CM_mod, Min,
	    sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
			(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
				Sbase for g in G, t in T) +
		sum(0 * (ΔP_R_up[t, r]) * Sbase +
			(price[t] - 0) * ΔP_R_dn[t, r] *
				Sbase for r in R, t in T) +
		sum((price[t]/storages.efficiency[s]) * P_S_up[t, s] * Sbase for s in S, t in T) +
		1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
		1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T));

	# Initiate optimisation process
	JuMP.optimize!(CM_mod)

	# Export results
	ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
	ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
	ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
	ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	price = JuMP.dual.(PowerInjectionBal)
	P_gen_lost_opt = JuMP.value.(P_gen_lost) * Sbase
	P_load_lost_opt = JuMP.value.(P_load_lost) * Sbase
	P_S_up_opt = JuMP.value.(P_S_up) * Sbase

	return(CM_mod,
	       ΔP_up_opt,
		   ΔP_dn_opt,
		   ΔP_R_up_opt,
		   ΔP_R_dn_opt,
	       Θ_opt,
	       P_flow_opt,
		   price,
		   P_gen_lost_opt,
		   P_load_lost_opt,
		   P_S_up_opt
	       )
end

function CongestionManagementPHS80_LoadShift(hours,
	nodes::Nodes,
	lines::Lines,
	powerplants::PowerPlants,
	renewables::Renewables,
	storages::Storages,
	price,
	P_opt,
	P_R_opt,
	P_S_opt,
	D_S_opt,
	L_S_opt,
	x_ag,
	x_re,
	x_ind,
	s_ag,
	s_re,
	s_ind,
	isworkhour

	)
""" parameters: distribution of ind/ag/re per node """

h,b = calc_h_b(lines, nodes)

T = hours
G = powerplants.unit
R = renewables.unit
S = storages.unit

slack = nodes.slack
Sbase = nodes.bmva # MVA

###################
# MODEL ###########
###################

# calculate average load per node
loadmean = Dict(i => sum(nodes.load[i][hours])/(length(hours)) for i in keys(nodes.load))

# Create a subset of buses I and J
J = nodes.id
L = lines.id

# Initialise JuMP model: Congestion management
CM_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

# Variables
@variable(CM_mod, ΔP_up[T, G] >= 0)
@variable(CM_mod, ΔP_dn[T, G] >= 0)
@variable(CM_mod, ΔP_R_up[T, R] == 0)
@variable(CM_mod, ΔP_R_dn[T, R] >= 0)
@variable(CM_mod, Θ[T, J])
@variable(CM_mod, P_flow[T, L])
@variable(CM_mod, P_load_lost[T, J] >= 0)
@variable(CM_mod, P_gen_lost[T, J] >= 0)

@variable(CM_mod, ΔD_up_ag[T, J] >= 0)
@variable(CM_mod, ΔD_up_re[T, J] >= 0)
@variable(CM_mod, ΔD_up_ind[T, J] >= 0)
@variable(CM_mod, ΔD_dn_ag[T, J] >= 0)
@variable(CM_mod, ΔD_dn_re[T, J] >= 0)
@variable(CM_mod, ΔD_dn_ind[T, J] >= 0)

@variable(CM_mod, ΔD_up[T, J] >= 0)
@variable(CM_mod, ΔD_dn[T, J] >= 0)

@variable(CM_mod, Load[T, J] >= 0)

@constraint(CM_mod, ShiftUpTotal[j = J, t = T], ΔD_up[t, j] == ΔD_up_ag[t, j] + ΔD_up_re[t, j] + ΔD_up_ind[t, j]);
@constraint(CM_mod, ShiftDnTotal[j = J, t = T], ΔD_dn[t, j] == ΔD_dn_ag[t, j] + ΔD_dn_re[t, j] + ΔD_dn_ind[t, j]);

@constraint(CM_mod, ShiftAgUpMax[j = J, t = T], ΔD_up_ag[t, j] <= x_ag*s_ag*loadmean[j]/Sbase);
@constraint(CM_mod, ShiftAgDnMax[j = J, t = T], ΔD_dn_ag[t, j] <= x_ag*s_ag*loadmean[j]/Sbase);
@constraint(CM_mod, ShiftReUpMax[j = J, t = T], ΔD_up_re[t, j] <= x_re*s_re*loadmean[j]/Sbase);
@constraint(CM_mod, ShiftReDnMax[j = J, t = T], ΔD_dn_re[t, j] <= x_re*s_re*loadmean[j]/Sbase);

# if mod(hours[24],168) < 120
@constraint(CM_mod, ShiftIndUpMax[j = J, t = T], ΔD_up_ind[t, j] <= x_ind*s_ind*loadmean[j]/Sbase*isworkhour[mod(t,24)]);
@constraint(CM_mod, ShiftIndDnMax[j = J, t = T], ΔD_dn_ind[t, j] <= x_ind*s_ind*loadmean[j]/Sbase*isworkhour[mod(t,24)]);
#	else
#@constraint(CM_mod, ShiftIndUpMax[j = J, t = T], ΔD_up_ind[t, j] == 0);
#@constraint(CM_mod, ShiftIndDnMax[j = J, t = T], ΔD_dn_ind[t, j] == 0);
# end

@constraint(CM_mod, ShiftBalAg[j = J],
sum(ΔD_up_ag[t, j] - ΔD_dn_ag[t, j] for t in hours) == 0);
@constraint(CM_mod, ShiftBalRe[j = J],
sum(ΔD_up_re[t, j] - ΔD_dn_re[t, j] for t in hours) == 0);
@constraint(CM_mod, ShiftBalInd[j = J],
sum(ΔD_up_ind[t, j] - ΔD_dn_ind[t, j] for t in hours) == 0);

@constraint(CM_mod, ShiftedLoad[j = J, t = T],
Load[t, j] == nodes.load[j][t]/Sbase + ΔD_up[t, j] - ΔD_dn[t, j]
);

# Storage voriables
@variable(CM_mod, P_S_up[T, S] >= 0)	# Power generation from storage unit
# @variable(ED_mod, D_S[T, S] >= 0)	# Electricity demand from storage unit
@variable(CM_mod, L_S[T, S] >= 0)

# Fix the voltage angle of the slack bus
for t in T
fix(Θ[t, slack], 0)
end

for s in S
fix(L_S[T[1], s], 0.2*storages.storage[s]/Sbase; force = true)
end

# Constraints
# Market clearing/power balance
# @constraint(CM_mod, MarketClearing[t = T],
#     sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) for g in G) +
# 	sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) for r in R) +
# 	sum(P_S_opt[t, s]/Sbase for s in S) -
# 	sum(D_S_opt[t, s]/Sbase for s in S)
# 		== nodes.systemload[t]/Sbase);

# Upper generation limit
@constraint(CM_mod, GenerationLimitUp[g = G, t = T],
(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) <=
powerplants.pmax[g][t]/Sbase);

# Lower generation limit
@constraint(CM_mod, GenerationLimitDown[g = G, t = T],
(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) >=
powerplants.pmin[g][t]/Sbase);

# Upper generation limit for renewables
@constraint(CM_mod, ResGenerationLimitUp[r = R, t = T],
(P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) <=
renewables.infeed[r][t]/Sbase);

# Lower generation limit for renewables
@constraint(CM_mod, ResGenerationLimitDown[r = R, t = T],
(P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) >=
0);

# # Ramp-up limit
# @constraint(CM_mod, RampUp[g = G, t = T; t > T[1]],
# 	(P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g]) <=
# 		powerplants.rup[g]/Sbase)
#
# # Ramp-down limit
# @constraint(CM_mod, RampDown[g = G, t = T; t > T[1]],
# 	-((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) -
# 		(P_opt[t-1, g]/Sbase + ΔP_up[t-1, g] - ΔP_dn[t-1, g])) <=
# 		powerplants.rdn[g]/Sbase)

### DC power flow constraints
# Power injection balance
@constraint(CM_mod, PowerInjectionBal[j = J, t = T],
sum((P_opt[t, g]/Sbase + ΔP_up[t, g] - ΔP_dn[t, g]) for g in getKeyVector(powerplants.node, j)) #Sum of power plant generation
+ sum((P_R_opt[t, r]/Sbase + ΔP_R_up[t, r] - ΔP_R_dn[t, r]) for r in getKeyVector(renewables.node, j)) #Sum of renewable generation
+ sum((P_S_opt[t, s]/Sbase + P_S_up[t, s]) for s in getKeyVector(storages.node, j)) #Sum of storage generation
- sum((D_S_opt[t, s]/Sbase) for s in getKeyVector(storages.node, j)) #Sum of storage demand
- P_gen_lost[t, j] #Lost generation
+ sum(b[j,jj] * Θ[t, jj] for jj in J) #Sum of line power flow?
==
# nodes.load[j][t]/Sbase +
Load[t, j]
- nodes.exchange[j][t]/Sbase # Always 0
- P_load_lost[t, j]); #Lost Load

@constraint(CM_mod, LostLoadBound[j = J, t = T],
P_load_lost[t, j] <= nodes.load[j][t]/Sbase);

@constraint(CM_mod, LinePowerFlowMax[l = L, t = T],
P_flow[t, l] <= lines.pmax[l]/Sbase);

@constraint(CM_mod, LinePowerFlowMin[l = L, t = T],
P_flow[t, l] >= -lines.pmax[l]/Sbase);

@constraint(CM_mod, LineFlow[l = L, t = T],
P_flow[t, l] == sum(h[l,j] * Θ[t,j] for j in J));

### PHS constraints
@constraint(CM_mod, StoragePowerOutput[s = S, t = T],
P_S_opt[t, s]/Sbase + P_S_up[t, s] <= storages.power[s]/Sbase);

# Maximum storage power input
# @constraint(CM_mod, StoragePowerInput[s = S, t = T],
# 	D_S[t, s] <= storages.power[s]);

# Generation CM_mod determined by remaining storage level
@constraint(CM_mod, StorageLevelGen[s = S, t = T],
P_S_up[t, s] <= L_S[t, s]);

# Storage
@constraint(CM_mod, Storage[s = S, t = T; t > T[1]],
L_S[t-1, s] - P_S_up[t-1, s] == L_S[t, s]);

# Objective function
@objective(CM_mod, Min,
sum(powerplants.mc[g][1] * (ΔP_up[t, g]) * Sbase +
(price[t] - powerplants.mc[g][1]) * ΔP_dn[t, g] *
Sbase for g in G, t in T) +
sum(0 * (ΔP_R_up[t, r]) * Sbase +
(price[t] - 0) * ΔP_R_dn[t, r] *
Sbase for r in R, t in T) +
sum((price[t]/storages.efficiency[s]) * P_S_up[t, s] * Sbase for s in S, t in T) +
1000 * sum(P_load_lost[t, j] * Sbase for j in J, t in T) +
1000 * sum(P_gen_lost[t, j] * Sbase for j in J, t in T));

# Initiate optimisation process
JuMP.optimize!(CM_mod)

# Export results
ΔP_up_opt = JuMP.value.(ΔP_up) * Sbase
ΔP_dn_opt = JuMP.value.(ΔP_dn) * Sbase
ΔP_R_up_opt = JuMP.value.(ΔP_R_up) * Sbase
ΔP_R_dn_opt = JuMP.value.(ΔP_R_dn) * Sbase
Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
P_flow_opt = JuMP.value.(P_flow)*Sbase
price = JuMP.dual.(PowerInjectionBal)
P_gen_lost_opt = JuMP.value.(P_gen_lost) * Sbase
P_load_lost_opt = JuMP.value.(P_load_lost) * Sbase
P_S_up_opt = JuMP.value.(P_S_up) * Sbase
Load_shifted_opt = JuMP.value.(Load) * Sbase
Loadshift_up_opt = JuMP.value.(ΔD_up) * Sbase
Loadshift_dn_opt = JuMP.value.(ΔD_dn) * Sbase

return(CM_mod,
ΔP_up_opt,
ΔP_dn_opt,
ΔP_R_up_opt,
ΔP_R_dn_opt,
Θ_opt,
P_flow_opt,
price,
P_gen_lost_opt,
P_load_lost_opt,
P_S_up_opt,
Load_shifted_opt,
Loadshift_up_opt,
Loadshift_dn_opt,
)
end


text:
module WorkerUtilities

export Lockable

const WORK_QUEUE = Channel{Task}(0)
const WORKER_TASKS = Task[]

"""
  WorkerUtilities.@spawn expr
  WorkerUtilities.@spawn passthroughstorage::Bool expr

Similar to `Threads.@spawn`, schedule and execute a task (given by `expr`)
that will be run on a "background worker" (see [`WorkerUtilities.init`]((@ref))).

In the 2-argument invocation, `passthroughstorage` controls whether the task-local storage of the
`current_task()` should be "passed through" to the spawned task.
"""
macro spawn(thunk)
    esc(quote
        tsk = @task $thunk
        tsk.storage = current_task().storage
        put!(WorkerUtilities.WORK_QUEUE, tsk)
        tsk
    end)
end

"""
  WorkerUtilities.@spawn expr
  WorkerUtilities.@spawn passthroughstorage::Bool expr

Similar to `Threads.@spawn`, schedule and execute a task (given by `expr`)
that will be run on a "background worker" (see [`WorkerUtilities.init`]((@ref))).

In the 2-argument invocation, `passthroughstorage` controls whether the task-local storage of the
`current_task()` should be "passed through" to the spawned task.
"""
macro spawn(passthroughstorage, thunk)
    esc(quote
        tsk = @task $thunk
        if $passthroughstorage
            tsk.storage = current_task().storage
        end
        put!(WorkerUtilities.WORK_QUEUE, tsk)
        tsk
    end)
end

"""
  WorkerUtilities.init(nworkers=Threads.nthreads() - 1)

Initialize background workers that will execute tasks spawned via
[`WorkerUtilities.@spawn`](@ref). If `nworkers == 1`, a single worker
will be started on thread 1 where tasks will be executed in contention
with other thread 1 work. Background worker tasks can be inspected by
looking at `WorkerUtilities.WORKER_TASKS`.
"""
function init(nworkers=Threads.nthreads()-1)
    maxthreadid = nworkers + 1
    tids = Threads.nthreads() == 1 ? (1:1) : 2:maxthreadid
    resize!(WORKER_TASKS, nworkers)
    Threads.@threads for tid in 1:maxthreadid
        if tid in tids
            WORKER_TASKS[tid == 1 ? 1 : (tid - 1)] = Base.@async begin
                for task in WORK_QUEUE
                    schedule(task)
                    wait(task)
                end
            end
        end
    end
    return
end

"""
  Lockable(value, lock = ReentrantLock())

Creates a `Lockable` object that wraps `value` and
associates it with the provided `lock`.
"""
struct Lockable{T, L <: Base.AbstractLock}
    value::T
    lock::L
end

Lockable(value) = Lockable(value, ReentrantLock())

"""
  lock(f::Function, l::Lockable)

Acquire the lock associated with `l`, execute `f` with the lock held,
and release the lock when `f` returns. `f` will receive one positional
argument: the value wrapped by `l`. If the lock is already locked by a
different task/thread, wait for it to become available.
When this function returns, the `lock` has been released, so the caller should
not attempt to `unlock` it.
"""
function Base.lock(f, l::Lockable)
    lock(l.lock) do
        f(l.value)
    end
end

# implement the rest of the Lock interface on Lockable
Base.islocked(l::Lockable) = islocked(l.lock)
Base.lock(l::Lockable) = lock(l.lock)
Base.trylock(l::Lockable) = trylock(l.lock)
Base.unlock(l::Lockable) = unlock(l.lock)

end # module


text:
include("fisher_yates!.jl")

"""
    simplexpt(unif)

Takes a vector of length `d-1` of numbers between `0` and `1` and converts it a point on the standard `d-1` dimensional simplex.
"""
function simplexpt(unif)
    d = length(unif) + 1
    T = eltype(unif)
    w = zeros(T, d + 1)
    w[2:d] .= sort(unif)
    w[d+1] = one(T)
    diff(w)
end

"""
    randprobvec([rng], d::Int)

Generates a random probability vector of length `d` using an algorithm by [Smith and Tromble](http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf).
"""
randprobvec(rng::AbstractRNG, ::Type{N}, d::Int) where {N<:Number} =
    simplexpt(rand(rng, N, d - 1))

randprobvec(rng::AbstractRNG, d::Int) = simplexpt(rand(rng, Float64, d - 1))

randprobvec(::Type{N}, d::Int) where {N<:Number} = randprobvec(Random.GLOBAL_RNG, N, d)

randprobvec(d::Int) = randprobvec(Random.GLOBAL_RNG, Float64, d)

"""
    randprobvec([rng], d::Int, N::T) where {T <: Integer}

Generates a random probability vector of length `d` whose entries are rational numbers with denominator at most `N` using an algorithm by [Smith and Tromble](http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf).
"""
function randprobvec(rng::AbstractRNG, d::Int, N::T) where {T<:Integer}
    unif = zeros(T, d - 1)
    fisher_yates_sample!(rng, 1:N, unif)
    unif = unif .// N
    simplexpt(unif)
end

randprobvec(d::Int, N::T) where {T<:Integer} = randprobvec(Random.GLOBAL_RNG, d, N)

"""
    randunitary([rng], d::Int)

Generates a unitary matrix of dimension `d` at random according to the Haar measure, using an algorithm described by Maris Ozols in ["How to generate a random unitary matrix"](http://home.lu.lv/~sd20008/papers/essays/Random%20unitary%20%5Bpaper%5D.pdf).
"""
function randunitary(rng::AbstractRNG, d::Int)
    RG = randn(rng, d, d) + im * randn(rng, d, d)
    Q, R = qr!(RG)
    r = diag(R)
    L = Diagonal(r ./ abs.(r))
    return Q * L
end

randunitary(d::Int) = randunitary(Random.GLOBAL_RNG, d)

"""
    randdm([rng], d::Int)

Generates a density matrix of dimension `d` at random.
"""
function randdm(rng, d::Int)
    eigs = Diagonal(randprobvec(rng, d))
    U = randunitary(rng, d)
    ρ = U * eigs * (U')
    return Hermitian((ρ + ρ') / 2)
end

randdm(d::Int) = randdm(Random.GLOBAL_RNG, d)

"""
    randmajpair([rng], [T], d::Int, [N])

Generates a pair of probability vectors `p`, `q` of length `d`
such that `p ≺ q` by rejection sampling. Supply a type `T` to
choose the numeric type of the elements of the vectors, and an
integer `N` to use rational numbers with maximum denominator `N`.
"""
function randmajpair(args...)
    while true
        p = randprobvec(args...)
        q = randprobvec(args...)
        p ≺ q && return p, q
        p ≻ q && return q, p
    end
end

"""
    randincomppair([rng], [T], d::Int, [N]; sortby = entropy)

Generates a pair of probability vectors `p`, `q` of length `d`
such that `!(p ≺ q)` and `!(q ≺ p)` by rejection sampling. Supply
a type `T` to choose the numeric type of the elements of the vectors,
and an integer `N` to use rational numbers with maximum denominator `N`.

Optionally, pass a function (from probability vectors to an
ordered set) to the keyword argument `sortby` to sort the results.
"""
function randincomppair(args...; sortby = entropy)
    while true
        p = randprobvec(args...)
        q = randprobvec(args...)
        p ≺ q && continue
        p ≻ q && continue
        if sortby(p) >= sortby(q)
            return p, q
        else
            return q, p
        end
    end
end


text:
module AdvancedTopics

import Distributions
import KernelDensity
import LinearAlgebra
import MLDatasets
import Plots
import Random

include("activation_functions.jl")
include("common.jl")
include("elm.jl")
include("sieve.jl")

include("mnist_data.jl")
include("synthetic_data.jl")
include("error.jl")

export ELM, Sieve, DataPassingSieve, ProjectionPassingSieve, predict,
       Cube, Linear, ReLU, Sigmoid, Square, Tanh,
       KDEComparator, train_kde_comparators,
       calculate_error,
       mnist_training_data, mnist_testing_data,
       synthetic_circular_data, synthetic_spiral_data

end # module


text:
using DataFrames
using PySide

## Test DataFrameModel.

## test roles
d = DataFrame(x=[randstring(10) for i in 1:10],
              x__ToolTipRole = [randstring(10) for i in 1:10],
              x__BackgroundRole= rep("yellow", 10),
              y = 1:10
              )

w = Widget()
lyt = VBoxLayout(w)
setLayout(w, lyt)

view = TableView(w)
push!(lyt, view)

m = DataFrameModel(d, view)
setModel(view, m)

## test scaling
n = 10^5
d1 = DataFrame(x = randn(n),
              y = 1:n)

view = TableView(w)
push!(lyt, view)

m = DataFrameModel(d1, view, editable=true)
setModel(view, m)

raise(w)


text:
module TestMetadata
    using Base.Test
    using Languages
    using TextAnalysis

    sample_text1 = "This is a string"
    sample_text2 = "This is also a string"
    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    sd = StringDocument(sample_text1)

    @assert isequal(name(sd), "Unnamed Document")
    @assert isequal(language(sd), TextAnalysis.Languages.EnglishLanguage)
    @assert isequal(author(sd), "Unknown Author")
    @assert isequal(timestamp(sd), "Unknown Time")
end