text:
module IntelGEMM

using MKL_jll

using LinearAlgebra, Libdl, Requires
using LinearAlgebra: chkstride1, checksquare, require_one_based_indexing
import LinearAlgebra.BLAS: gemm!, gemv!
import LinearAlgebra.LAPACK: chkargsok, chklapackerror, chknonsingular, chkposdef, chkfinite
import LinearAlgebra.LAPACK: geevx!, ggev!, syev!, syevr!, sygvd!

const BlasInt = Int32
const libblas = libmkl_rt
const liblapack = libmkl_rt

@enum Threading begin
    THREADING_INTEL
    THREADING_SEQUENTIAL
    THREADING_PGI
    THREADING_GNU
    THREADING_TBB
end

@enum Interface begin
    INTERFACE_LP64
    INTERFACE_ILP64
    INTERFACE_GNU
end

function set_threading_layer(layer::Threading = THREADING_INTEL)
    err = ccall((:MKL_Set_Threading_Layer, libblas), Cint, (Cint,), layer)
    if err == -1
        throw(ErrorException("return value was -1"))
    end
    return nothing
end

function set_interface_layer(interface = INTERFACE_LP64)
    err = ccall((:MKL_Set_Interface_Layer, libblas), Cint, (Cint,), interface)
    if err == -1
        throw(ErrorException("return value was -1"))
    end
    return nothing
end

function __init__()
    set_threading_layer()
    set_interface_layer()
    @require NNlib="872c559c-99b0-510c-b3b7-b6c96a88d5cd" include("nnlib.jl")
end

macro blasfunc(x)
    return Expr(:quote, x)
end

for chk in (:chkargsok, :chklapackerror, :chknonsingular, :chkposdef)
    @eval $chk(ret) = $chk(BLAS.BlasInt(ret))
end

for (gemm, elty) in
    ((:dgemm_,:Float64),
     (:sgemm_,:Float32),
     (:zgemm_,:ComplexF64),
     (:cgemm_,:ComplexF32))
@eval begin
         # SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
         # *     .. Scalar Arguments ..
         #       DOUBLE PRECISION ALPHA,BETA
         #       INTEGER K,LDA,LDB,LDC,M,N
         #       CHARACTER TRANSA,TRANSB
         # *     .. Array Arguments ..
         #       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
    function gemm!(transA::AbstractChar, transB::AbstractChar,
                   alpha::Union{($elty), Bool},
                   A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty},
                   beta::Union{($elty), Bool},
                   C::AbstractVecOrMat{$elty})
    #           if any([stride(A,1), stride(B,1), stride(C,1)] .!= 1)
    #               error("gemm!: BLAS module requires contiguous matrix columns")
    #           end  # should this be checked on every call?
            require_one_based_indexing(A, B, C)
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end
            chkstride1(A)
            chkstride1(B)
            chkstride1(C)
            ccall((@blasfunc($gemm), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                Ref{BlasInt}),
                transA, transB, m, n,
                ka, alpha, A, max(1,stride(A,2)),
                B, max(1,stride(B,2)), beta, C,
                max(1,stride(C,2)))
            C
        end
    end
end

for (fname, elty) in ((:dgemv_,:Float64),
                      (:sgemv_,:Float32),
                      (:zgemv_,:ComplexF64),
                      (:cgemv_,:ComplexF32))
    @eval begin
             #SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             #*     .. Scalar Arguments ..
             #      DOUBLE PRECISION ALPHA,BETA
             #      INTEGER INCX,INCY,LDA,M,N
             #      CHARACTER TRANS
             #*     .. Array Arguments ..
             #      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function gemv!(trans::AbstractChar, alpha::Union{($elty), Bool},
                       A::AbstractVecOrMat{$elty}, X::AbstractVector{$elty},
                       beta::Union{($elty), Bool}, Y::AbstractVector{$elty})
            require_one_based_indexing(A, X, Y)
            m,n = size(A,1),size(A,2)
            if trans == 'N' && (length(X) != n || length(Y) != m)
                throw(DimensionMismatch("A has dimensions $(size(A)), X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'C' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("the adjoint of A has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'T' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("the transpose of A has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 trans, size(A,1), size(A,2), alpha,
                 A, max(1,stride(A,2)), X, stride(X,1),
                 beta, Y, stride(Y,1))
            Y
        end
    end
end

for (geevx, ggev, elty) in
    ((:dgeevx_,:dggev_,:Float64),
     (:sgeevx_,:sggev_,:Float32))
    @eval begin
        #     SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
        #                          VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
        #                          RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
        #       INTEGER            IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
        #       DOUBLE PRECISION   ABNRM
        #       ..
        #       .. Array Arguments ..
        #       INTEGER            IWORK( * )
        #       DOUBLE PRECISION   A( LDA, * ), RCONDE( * ), RCONDV( * ),
        #      $                   SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
        #      $                   WI( * ), WORK( * ), WR( * )
        function geevx!(balanc::AbstractChar, jobvl::AbstractChar, jobvr::AbstractChar, sense::AbstractChar, A::AbstractMatrix{$elty})
            n = checksquare(A)
            chkfinite(A) # balancing routines don't support NaNs and Infs
            lda = max(1,stride(A,2))
            wr = similar(A, $elty, n)
            wi = similar(A, $elty, n)
            if balanc ∉ ['N', 'P', 'S', 'B']
                throw(ArgumentError("balanc must be 'N', 'P', 'S', or 'B', but $balanc was passed"))
            end
            ldvl = 0
            if jobvl == 'V'
                ldvl = n
            elseif jobvl == 'N'
                ldvl = 0
            else
                throw(ArgumentError("jobvl must be 'V' or 'N', but $jobvl was passed"))
            end
            VL = similar(A, $elty, ldvl, n)
            ldvr = 0
            if jobvr == 'V'
                ldvr = n
            elseif jobvr == 'N'
                ldvr = 0
            else
                throw(ArgumentError("jobvr must be 'V' or 'N', but $jobvr was passed"))
            end
            VR = similar(A, $elty, ldvr, n)
            ilo = Ref{BlasInt}()
            ihi = Ref{BlasInt}()
            scale = similar(A, $elty, n)
            abnrm = Ref{$elty}()
            rconde = similar(A, $elty, n)
            rcondv = similar(A, $elty, n)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            iworksize = 0
            if sense == 'N' || sense == 'E'
                iworksize = 0
            elseif sense == 'V' || sense == 'B'
                iworksize = 2*n - 2
            else
                throw(ArgumentError("sense must be 'N', 'E', 'V' or 'B', but $sense was passed"))
            end
            iwork = Vector{BlasInt}(undef, iworksize)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($geevx), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                       Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                       Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                       Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                       Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
                       Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                       balanc, jobvl, jobvr, sense,
                       n, A, lda, wr,
                       wi, VL, max(1,ldvl), VR,
                       max(1,ldvr), ilo, ihi, scale,
                       abnrm, rconde, rcondv, work,
                       lwork, iwork, info)
                # chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                end
            end
            A, wr, wi, VL, VR, ilo[], ihi[], scale, abnrm[], rconde, rcondv
        end

        #       SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
        #      $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBVL, JOBVR
        #       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
        # *     ..
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
        #      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
        #      $                   VR( LDVR, * ), WORK( * )
        function ggev!(jobvl::AbstractChar, jobvr::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            require_one_based_indexing(A, B)
            chkstride1(A,B)
            n, m = checksquare(A,B)
            if n != m
                throw(DimensionMismatch("A has dimensions $(size(A)), and B has dimensions $(size(B)), but A and B must have the same size"))
            end
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            alphar = similar(A, $elty, n)
            alphai = similar(A, $elty, n)
            beta = similar(A, $elty, n)
            ldvl = 0
            if jobvl == 'V'
                ldvl = n
            elseif jobvl == 'N'
                ldvl = 1
            else
                throw(ArgumentError("jobvl must be 'V' or 'N', but $jobvl was passed"))
            end
            vl = similar(A, $elty, ldvl, n)
            ldvr = 0
            if jobvr == 'V'
                ldvr = n
            elseif jobvr == 'N'
                ldvr = 1
            else
                throw(ArgumentError("jobvr must be 'V' or 'N', but $jobvr was passed"))
            end
            vr = similar(A, $elty, ldvr, n)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($ggev), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                     Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                     Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
                     Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                     Ptr{BlasInt}),
                    jobvl, jobvr, n, A,
                    lda, B, ldb, alphar,
                    alphai, beta, vl, ldvl,
                    vr, ldvr, work, lwork,
                    info)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                end
            end
            alphar, alphai, beta, vl, vr
        end
    end
end

for (geevx, ggev, elty, relty) in
    ((:zgeevx_,:zggev_,:ComplexF64,:Float64),
     (:cgeevx_,:cggev_,:ComplexF32,:Float32))
    @eval begin
        #     SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
        #                          LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
        #                          RCONDV, WORK, LWORK, RWORK, INFO )
        #
        #       .. Scalar Arguments ..
        #       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
        #       INTEGER            IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
        #       DOUBLE PRECISION   ABNRM
        #       ..
        #       .. Array Arguments ..
        #       DOUBLE PRECISION   RCONDE( * ), RCONDV( * ), RWORK( * ),
        #      $                   SCALE( * )
        #       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
        #      $                   W( * ), WORK( * )
        function geevx!(balanc::AbstractChar, jobvl::AbstractChar, jobvr::AbstractChar, sense::AbstractChar, A::AbstractMatrix{$elty})
            n = checksquare(A)
            chkfinite(A) # balancing routines don't support NaNs and Infs
            lda = max(1,stride(A,2))
            w = similar(A, $elty, n)
            if balanc ∉ ['N', 'P', 'S', 'B']
                throw(ArgumentError("balanc must be 'N', 'P', 'S', or 'B', but $balanc was passed"))
            end
            ldvl = 0
            if jobvl == 'V'
                ldvl = n
            elseif jobvl == 'N'
                ldvl = 0
            else
                throw(ArgumentError("jobvl must be 'V' or 'N', but $jobvl was passed"))
            end
            VL = similar(A, $elty, ldvl, n)
            ldvr = 0
            if jobvr == 'V'
                ldvr = n
            elseif jobvr == 'N'
                ldvr = 0
            else
                throw(ArgumentError("jobvr must be 'V' or 'N', but $jobvr was passed"))
            end
            if sense ∉ ['N','E','V','B']
                throw(ArgumentError("sense must be 'N', 'E', 'V' or 'B', but $sense was passed"))
            end
            VR = similar(A, $elty, ldvr, n)
            ilo = Ref{BlasInt}()
            ihi = Ref{BlasInt}()
            scale = similar(A, $relty, n)
            abnrm = Ref{$relty}()
            rconde = similar(A, $relty, n)
            rcondv = similar(A, $relty, n)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            rwork = Vector{$relty}(undef, 2n)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($geevx), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                       Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                       Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty},
                       Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt},
                       Ptr{$relty}, Ptr{BlasInt}),
                       balanc, jobvl, jobvr, sense,
                       n, A, lda, w,
                       VL, max(1,ldvl), VR, max(1,ldvr),
                       ilo, ihi, scale, abnrm,
                       rconde, rcondv, work, lwork,
                       rwork, info)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                end
            end
            A, w, VL, VR, ilo[], ihi[], scale, abnrm[], rconde, rcondv
        end

        # SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
        #      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBVL, JOBVR
        #       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
        # *     ..
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION   RWORK( * )
        #       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
        #      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
        #      $                   WORK( * )
        function ggev!(jobvl::AbstractChar, jobvr::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            require_one_based_indexing(A, B)
            chkstride1(A, B)
            n, m = checksquare(A, B)
            if n != m
                throw(DimensionMismatch("A has dimensions $(size(A)), and B has dimensions $(size(B)), but A and B must have the same size"))
            end
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            alpha = similar(A, $elty, n)
            beta = similar(A, $elty, n)
            ldvl = 0
            if jobvl == 'V'
                ldvl = n
            elseif jobvl == 'N'
                ldvl = 1
            else
                throw(ArgumentError("jobvl must be 'V' or 'N', but $jobvl was passed"))
            end
            vl = similar(A, $elty, ldvl, n)
            ldvr = 0
            if jobvr == 'V'
                ldvr = n
            elseif jobvr == 'N'
                ldvr = 1
            else
                throw(ArgumentError("jobvr must be 'V' or 'N', but $jobvr was passed"))
            end
            vr = similar(A, $elty, ldvr, n)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            rwork = Vector{$relty}(undef, 8n)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($ggev), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                     Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                     Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                     Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$relty},
                     Ptr{BlasInt}),
                    jobvl, jobvr, n, A,
                    lda, B, ldb, alpha,
                    beta, vl, ldvl, vr,
                    ldvr, work, lwork, rwork,
                    info)
                # chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                end
            end
            alpha, beta, vl, vr
        end
    end
end

# Symmetric (real) eigensolvers
for (syev, syevr, sygvd, elty) in
    ((:dsyev_,:dsyevr_,:dsygvd_,:Float64),
     (:ssyev_,:ssyevr_,:ssygvd_,:Float32))
    @eval begin
        #       SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBZ, UPLO
        #       INTEGER            INFO, LDA, LWORK, N
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
        function syev!(jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix{$elty})
            chkstride1(A)
            n = checksquare(A)
            W     = similar(A, $elty, n)
            work  = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            info  = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($syev), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                      Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                      jobz, uplo, n, A, max(1,stride(A,2)), W, work, lwork, info)
                # chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            jobz == 'V' ? (W, A) : W
        end

        #       SUBROUTINE DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
        #      $                   ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
        #      $                   IWORK, LIWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBZ, RANGE, UPLO
        #       INTEGER            IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
        #       DOUBLE PRECISION   ABSTOL, VL, VU
        # *     ..
        # *     .. Array Arguments ..
        #       INTEGER            ISUPPZ( * ), IWORK( * )
        #       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
        function syevr!(jobz::AbstractChar, range::AbstractChar, uplo::AbstractChar, A::AbstractMatrix{$elty},
                        vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat)
            chkstride1(A)
            n = checksquare(A)
            if range == 'I' && !(1 <= il <= iu <= n)
                throw(ArgumentError("illegal choice of eigenvalue indices (il = $il, iu = $iu), which must be between 1 and n = $n"))
            end
            if range == 'V' && vl >= vu
                throw(ArgumentError("lower boundary, $vl, must be less than upper boundary, $vu"))
            end
            lda = stride(A,2)
            m = Ref{BlasInt}()
            w = similar(A, $elty, n)
            ldz = n
            if jobz == 'N'
                Z = similar(A, $elty, ldz, 0)
            elseif jobz == 'V'
                Z = similar(A, $elty, ldz, n)
            end
            isuppz = similar(A, BlasInt, 2*n)
            work   = Vector{$elty}(undef, 1)
            lwork  = BlasInt(-1)
            iwork  = Vector{BlasInt}(undef, 1)
            liwork = BlasInt(-1)
            info   = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                ccall((@blasfunc($syevr), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                        Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ref{$elty},
                        Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{BlasInt},
                        Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                        Ptr{BlasInt}),
                    jobz, range, uplo, n,
                    A, max(1,lda), vl, vu,
                    il, iu, abstol, m,
                    w, Z, max(1,ldz), isuppz,
                    work, lwork, iwork, liwork,
                    info)
                # chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                    liwork = iwork[1]
                    resize!(iwork, liwork)
                end
            end
            w[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
        end
        syevr!(jobz::AbstractChar, A::AbstractMatrix{$elty}) =
            syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0)

        # Generalized eigenproblem
        #           SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
        #      $                   LWORK, IWORK, LIWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBZ, UPLO
        #       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
        # *     ..
        # *     .. Array Arguments ..
        #       INTEGER            IWORK( * )
        #       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
        function sygvd!(itype::Integer, jobz::AbstractChar, uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            chkstride1(A, B)
            n, m = checksquare(A, B)
            if n != m
                throw(DimensionMismatch("dimensions of A, ($n,$n), and B, ($m,$m), must match"))
            end
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            w = similar(A, $elty, n)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            iwork = Vector{BlasInt}(undef, 1)
            liwork = BlasInt(-1)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                ccall((@blasfunc($sygvd), liblapack), Cvoid,
                    (Ref{BlasInt}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                     Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                     Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt},
                     Ref{BlasInt}, Ptr{BlasInt}),
                    itype, jobz, uplo, n,
                    A, lda, B, ldb,
                    w, work, lwork, iwork,
                    liwork, info)
                chkargsok(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                    liwork = iwork[1]
                    resize!(iwork, liwork)
                end
            end
            chkposdef(info[])
            w, A, B
        end
    end
end

end # module


text:
# function to generate an adequate DataFrame with capacities for each
# weekhour and district center if no external capacity plan is provided
 function capacity_heuristic(incidents::DataFrame,
                                sim_data::DataFrame,
                                total_capacity::Float64,
                                min_capacity::Int64)
    # create a DataFrame with time data
        incidents_time = DataFrame(incidentid   = Int[],
        location    = Int[],
        year        = Int[],
        yearweek    = Int[],
        weekhour    = Int[],
        workload    = Float64[])
    # fill the DataFrame with the incident data
        for i = 1:nrow(incidents)
        push!(incidents_time, (incidentid   = incidents[i, :incidentid],
            location    = sim_data[i,:location_responsible],
            year        = year(unix2datetime(incidents[i, :epoch])),
            yearweek    = week(unix2datetime(incidents[i, :epoch])),
            weekhour    = incidents[i, :weekhour],
            workload = 0))
        end
    # calculate the overall workload in minutes per incident (including
    # variability of the driving time)
        for i = 1:nrow(incidents)
            for j = 1:incidents[i,:cars]
            incidents_time[i,:workload] +=
            dispatch_drivingtime(drivingtime[sim_data[i,:location_responsible],
            incidents[i,:location]],traffic,sim_data[i,:incident_minute])
            incidents_time[i,:workload] += incidents[i,:length]
            incidents_time[i,:workload] +=
            dispatch_drivingtime(drivingtime[sim_data[i,:location_responsible],
            incidents[i,:location]],traffic,sim_data[i,:incident_minute])
            end
            incidents_time[i,:workload] += incidents[i,:backlog]
        end
    # group the data and create temporary DataFrames to include each weekhour
    # during the incident data set timeframe in each location of a district center
        incidents_time = groupby(incidents_time,[:location,:year,:yearweek,:weekhour])
        incidents_time = combine(incidents_time, :workload  => sum => :workload)
        incidents_time = groupby(incidents_time,[:location,:yearweek,:weekhour])
        incidents_time = combine(incidents_time, :workload  => mean => :workload)
        time_frame = DataFrame(Array{Int64,2}(undef,length(unique(incidents_time[!,:yearweek])),1),[:yearweek])
        time_frame[:,:yearweek] .= unique(incidents_time[!,:yearweek])
        weekhours = DataFrame(Array{Int64,2}(undef,168,1),[:weekhour])
        for i = 1:168
            weekhours[i,:weekhour] = i
        end
        time_frame = crossjoin(time_frame, weekhours, makeunique = true)
        all_locations = unique!(sim_data[:,[:location_responsible]], :location_responsible)
        all_locations = rename!(all_locations,[:location])
        time_frame = crossjoin(time_frame, all_locations, makeunique = true)
        time_frame[:,:shift] .= 0
        for i = 1:nrow(time_frame)
            time_frame[i,:shift] = shifts[time_frame[i, :weekhour]]
        end
        incidents_time = outerjoin(incidents_time,time_frame,on=[:location,:yearweek,:weekhour])
        incidents_time = coalesce.(incidents_time, 0.0)
        incidents_time = sort!(incidents_time,[:location,:yearweek,:weekhour])
    # group the data and include a safety buffer as stated in our article
        incidents_time = groupby(incidents_time,[:location,:weekhour])
        incidents_time = combine(incidents_time,
        :shift => maximum => :shift,
        :workload => mean => :workload_mean,
        :workload => std =>  :workload_std)
        incidents_time = coalesce.(incidents_time, 0.0)
        incidents_time = sort!(incidents_time, [:location,:weekhour,:shift])
        incidents_time[:,:workload] = incidents_time[:,:workload_mean] #+
            incidents_time[:,:workload_std] *
            quantile.(Normal(), capacity_service)
    # group the resulting workload to generate a mean hourly workload for each
    # district per shift
        shifts_time = groupby(incidents_time,[:location, :shift])
        shifts_time = combine(shifts_time,:workload => maximum => :workload, nrow => :shift_hours)
        shift_length = groupby(shifts_time, :shift)
        shift_length = combine(shift_length, :shift_hours => maximum => :shift_hours)
        shifts_time = select!(shifts_time, Not([:shift_hours]))
        shifts_time = unstack(shifts_time,:location,:workload)
        shifts_capacity = copy(shifts_time)
        shifts_out = copy(shifts_time)
    # assign the available ressources to the shifts during the course of a week
        shifts_capacity = shifts_capacity[:,2:end] .= min_capacity
        shifts_time = shifts_time[:,2:end] = shifts_time[:,2:end] .- 60 * min_capacity
        shifts_capacity = Array(shifts_capacity)
        shifts_time = Array(shifts_time)
        shift_length = Array(shift_length)
        all_capacity = (total_capacity - min_capacity * size(shifts_time,2)) * 168
        while all_capacity >= maximum(shift_length[:,2])
            highest_capacity_pressure = findmax(shifts_time[:,:])[2]
            shifts_time[highest_capacity_pressure] -= 60
            shifts_capacity[highest_capacity_pressure] += 1
            all_capacity -= shift_length[highest_capacity_pressure[1],2]
        end
    # export the resulting DataFrame with the planned capacity per weekhour
    # and district
        shifts_out[:,2:end] = shifts_capacity
        weekhours[:,:shift] .= shifts[:,1]
        shifts_out = innerjoin(weekhours, shifts_out,on=:shift)
        shifts_out = sort!(shifts_out, [:weekhour,:shift])
        simulation_capacity = Array{Int64,2}(shifts_out[:,3:end])
        shifts_out = sort!(stack(shifts_out, 3:size(shifts_out,2)), [:weekhour])
    return simulation_capacity, shifts_out
end

ext: jl


text:
# Helper functions for external sort

"""
Write N random 4-byte integers to file (in binary form).
"""
function create_file(n, file_name="sample.bin")
    srand(13)  # Remove in the future
    file = open(file_name, "w")
    for i in 1:n
        write(file, rand(Int32))
    end
    close(file)
end

"""
Read file and print it to screen.
"""
function read_file2screen(file_name="sample.bin", print_name=false)
    if print_name
        println(string("Begin of ", file_name, ":"))
    end

    file = open(file_name, "r")
    while ~ eof(file)
        println(read(file, Int32))
    end
    close(file)

    if print_name
        println(string("End of ", file_name, ":"))
    end
end

"""
Read file to array.
"""
function read_file2array(file_name="sample.bin")
    file = open(file_name, "r")
    array2return = Array(Int32, 0)
    while ~ eof(file)
        push!(array2return, read(file, Int32))
    end
    close(file)
    return array2return
end

"""
Check external sort.
"""
function check_sort(input_filename="sample.bin", sorted_filename="sorted-sample.bin")
    input = read_file2array(input_filename)
    sort!(input)

    sorted = read_file2array(sorted_filename)

    check_result = input == sorted

    if ~ check_result
        println(string("sorted input: ", input))
        println(string("sorted file:  ", sorted))
    end

    return check_result
end

max_stars_repo_name: rgaiacs/external-sort

max_issues_repo_name: rgaiacs/external-sort

max_forks_repo_name: rgaiacs/external-sort


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# This is an interface to FirstOrderLp for solving a single QP or LP and writing
# the solution and solve statistics to a file. Run with --help for a pretty
# description of the arguments. See the comments on solve_instance_and_output
# for a description of the output formats.

# Built-in Julia packages.
import SparseArrays

# Third-party Julia packages.
import ArgParse
import GZip
import JSON3

import FirstOrderLp

function run_with_redirected_stdio(func, stdout_path, stderr_path)
  open(stdout_path, "w") do stdout_io
    open(stderr_path, "w") do stderr_io
      redirect_stdout(stdout_io) do
        redirect_stderr(stderr_io) do
          func()
        end
      end
    end
  end
end

function write_vector_to_file(filename, vector)
  open(filename, "w") do io
    for x in vector
      println(io, x)
    end
  end
end

"""
Solves a linear or quadratic programming problem using one of the methods
implemented in FirstOrderLp. Takes a path to an instance. The instance must have
the extension .mps, .mps.gz, .qps, or .qps.gz.
Creates:
- `instance_full_log.json.gz` with a SolveLog in gzipped JSON format with all
   iteration_stats data returned by the solver
- `instance_summary.json` with a SolveLog in plain text JSON format, identical
   to `instance_full_log.json.gz` but with iteration_stats cleared
- `instance_primal.txt` with the primal solution
- `instance_dual.txt` with the dual solution

If `redirect_stdio` is `True`, then it also writes `instance_stderr.txt`
and `instance_stdout.txt`.
"""
function solve_instance_and_output(
  parameters::Union{
    FirstOrderLp.MirrorProxParameters,
    FirstOrderLp.PdhgParameters,
  },
  output_dir::String,
  instance_path::String,
  redirect_stdio::Bool,
  transform_bounds::Bool,
  fixed_format_input::Bool,
)
  if !isdir(output_dir)
    mkpath(output_dir)
  end

  instance_name =
    replace(basename(instance_path), r"\.(mps|MPS|qps|QPS)(\.gz)?$" => "")

  function inner_solve()
    lower_file_name = lowercase(basename(instance_path))
    if endswith(lower_file_name, ".mps") ||
       endswith(lower_file_name, ".mps.gz") ||
       endswith(lower_file_name, ".qps") ||
       endswith(lower_file_name, ".qps.gz")
      lp = FirstOrderLp.qps_reader_to_standard_form(
        instance_path,
        fixed_format = fixed_format_input,
      )
    else
      error(
        "Instance has unrecognized file extension: ",
        basename(instance_path),
      )
    end

    presolve_info = FirstOrderLp.presolve(
      lp;
      verbosity = parameters.verbosity,
      transform_bounds = transform_bounds,
    )

    if parameters.verbosity >= 1
      println("Instance: ", instance_name)
    end
    running_time = @elapsed begin
      output::FirstOrderLp.SaddlePointOutput =
        FirstOrderLp.optimize(parameters, lp)
    end
    println("Elapsed time: $running_time sec")

    log = FirstOrderLp.SolveLog()
    log.instance_name = instance_name
    log.command_line_invocation = join([PROGRAM_FILE; ARGS...], " ")
    log.termination_reason = output.termination_reason
    log.termination_string = output.termination_string
    log.iteration_count = output.iteration_count
    log.solve_time_sec = running_time
    # This assumes that the last iterate matches the solution returned by the
    # solver.
    log.solution_stats = output.iteration_stats[end]
    # TODO: Update this once we return more than one type.
    log.solution_type = FirstOrderLp.POINT_TYPE_AVERAGE_ITERATE

    summary_output_path = joinpath(output_dir, instance_name * "_summary.json")
    open(summary_output_path, "w") do io
      write(io, JSON3.write(log, allow_inf = true))
    end

    log.iteration_stats = output.iteration_stats
    full_log_output_path =
      joinpath(output_dir, instance_name * "_full_log.json.gz")
    GZip.open(full_log_output_path, "w") do io
      write(io, JSON3.write(log, allow_inf = true))
    end

    primal_solution, dual_solution = FirstOrderLp.undo_presolve(
      presolve_info,
      output.primal_solution,
      output.dual_solution,
    )

    primal_output_path = joinpath(output_dir, instance_name * "_primal.txt")
    write_vector_to_file(primal_output_path, primal_solution)

    dual_output_path = joinpath(output_dir, instance_name * "_dual.txt")
    write_vector_to_file(dual_output_path, dual_solution)
  end

  if redirect_stdio
    stdout_path = joinpath(output_dir, instance_name * "_stdout.txt")
    stderr_path = joinpath(output_dir, instance_name * "_stderr.txt")
    run_with_redirected_stdio(inner_solve, stdout_path, stderr_path)
  else
    inner_solve()
  end

  return
end

"""
Defines parses and args.

# Returns
A dictionary with the values of the command-line arguments.
"""
function parse_command_line()
  arg_parse = ArgParse.ArgParseSettings()

  help_method = "The optimization method to use, must be `mirror-prox` or `pdhg`."

  help_instance_path = "The path to the instance to solve in .mps.gz or .mps format."

  ArgParse.@add_arg_table! arg_parse begin
    "--method"
    help = help_method
    arg_type = String
    required = true

    "--output_dir"
    help = "The directory for output files."
    arg_type = String
    required = true

    "--instance_path"
    help = help_instance_path
    arg_type = String
    required = true

    "--l_inf_ruiz_iterations"
    help =
      "Number of l_infinity Ruiz rescaling iterations to apply to the " *
      "constraint matrix. Zero disables this rescaling pass."
    arg_type = Int
    default = 10

    "--l2_norm_rescaling"
    help = "If true, applies L2 norm rescaling after Ruiz rescaling."
    arg_type = Bool
    default = false

    "--pock_chambolle_rescaling"
    help =
      "If true, Pock-Chambolle rescaling is applied after the (optional) " *
      "L2 norm rescaling using the parameter pock_chambolle_alpha. If false, " *
      "this rescaling step is skipped."
    arg_type = Bool
    default = true

    "--pock_chambolle_alpha"
    help =
      "The exponent parameter alpha used if Pock-Chambolle rescaling is " *
      "applied. Alpha must be in the interval [0, 2]."
    arg_type = Float64
    default = 1.0

    "--primal_importance"
    help =
      "Must be positive. The default of 1 balances primal and dual " *
      "~equally. This parameter biases the initial value of the primal " *
      "weight. Note that tuning this parameter value to a particular problem " *
      "may significantly improve performance. "
    arg_type = Float64
    default = 1.0

    "--scale_invariant_initial_primal_weight"
    help =
      "If true, uses a scale-invariant choice of the initial primal weight " *
      "biased by primal_importance. If false, primal_importance is the " *
      "primal weight."
    arg_type = Bool
    default = true

    "--artificial_restart_threshold"
    help =
      "If in the past artificial_restart_threshold fraction of " *
      "iterations no restart has occurred then a restart will be " *
      "artificially triggered. The value should be strictly greater than " *
      "zero and less than or equal to one. " *
      "Smaller values will have more frequent artificial restarts than " *
      "larger values. A value of one means there is exactly one " *
      "artificial restart occurs (on iteration 1)."
    arg_type = Float64
    default = 0.5

    # TODO: Make more detailed description when restart scheme is
    # finalized.
    "--sufficient_reduction_for_restart"
    help =
      "Only supported when restart_scheme=adaptive_normalized. It is " *
      "the threshold improvement in the quality of the current/average " *
      "iterate compared with that of the last restart that will trigger a " *
      "restart. The value of this parameter should be between zero and one. " *
      "Smaller values make restarts less frequent, larger values make " *
      "restarts more frequent."
    arg_type = Float64
    default = 0.1

    "--necessary_reduction_for_restart"
    help =
      "Only supported when restart_scheme is adaptive_normalized, " *
      "adaptive_distance, or adaptive_localized. It is the threshold " *
      "improvement in the quality of the current/average iterate compared " *
      "with that of the last restart that is necessary to trigger a restart. " *
      "If this improvement threshold is met and the quality of the iterates " *
      "appear to be getting worse then a restart is triggered. " *
      "The value of this parameter should be between the value of " *
      "sufficient_reduction_for_restart and one. " *
      "Smaller values make restarts less frequent, larger values make " *
      "restarts more frequent."
    arg_type = Float64
    default = 0.9

    "--primal_weight_update_smoothing"
    help =
      "This parameter controls exponential smoothing of " *
      "log(primal_weight) when updating the primal weight. " *
      "Must be between 0.0 and 1.0 inclusive. At 0.0, the primal " *
      "weight is frozen at the initial value given by the " *
      "`primal_importance` argument. At 1.0, no smoothing is performed. It " *
      "may be preferable to freeze the primal weight if a good estimate is " *
      "available."
    arg_type = Float64
    default = 0.5

    "--verbosity"
    help =
      "The verbosity level for printing. Values between 1 and 4 " *
      "print generic information on the solve process, i.e., a table of the " *
      "iteration statistics. Values greater than 5 provide information " *
      "useful for developers."
    arg_type = Int64
    default = 2

    "--redirect_stdio"
    help = "Redirect stdout and stderr to files (for batch runs)."
    arg_type = Bool
    default = false

    "--diagonal_scaling"
    help =
      "Mirror-prox only. Supported {off, l1, l2}. Use a diagonal matrix to " *
      "define the Bregman distance or, equivalently, rescale the primal " *
      "and dual variables individually."
    arg_type = String
    default = "off"

    "--restart_scheme"
    help =
      "Supported steps: {no_restart, fixed_frequency, " *
      "adaptive_normalized, adaptive_localized, adaptive_distance}. " *
      "See FirstOrderLp.RestartScheme enum for detailed documentation."
    arg_type = String
    default = "adaptive_normalized"

    "--restart_frequency"
    help =
      "Only relevant if --restart_scheme = fixed_frequency. " *
      "Determines the number of iterations until restart."
    arg_type = Int64
    default = 1000

    "--restart_to_current_metric"
    help =
      "Options: {gap_over_distance, gap_over_distance_squared, " *
      "no_restart_to_current}. If this is value is no_restart_to_current " *
      "then always restart to the average. Otherwise, we dynamically decide " *
      "whether we reset to the average or current. There are two options for " *
      "making this decision: gap_over_distance_squared and " *
      "gap_over_distance. gap_over_distance_squared will reset to current " *
      "if the normalized duality gap divided by the distanced travelled for " *
      "the current iterate is better. Alternatively, gap_over_distance" *
      "restarts to the current if it has a better normalized duality gap than" *
      "the average iterate. " *
      "Also, note that this option has no impact if " *
      "restart_scheme=no_restart."
    arg_type = String
    default = "gap_over_distance_squared"

    "--use_approximate_localized_duality_gap"
    help =
      "Whether to use an approximate localized duality gap in the " *
      "restart scheme."
    arg_type = Bool
    default = false

    "--record_iteration_stats"
    help =
      "Whether we record iterations stats. If true then record an " *
      "IterationStats object with frequency (in iterations) " *
      "equal to termination_evaluation_frequency."
    arg_type = Bool
    default = true

    "--termination_evaluation_frequency"
    help =
      "Frequency (in iterations) that the termination criteron is " *
      "evaluated."
    arg_type = Int64
    default = 40

    # The following parameters specify termination criteria.
    # These correspond to the fields in the TerminationCriteria struct.
    # If the arguments are not provided, the default values are grabbed from
    # construct_termination_criteria in termination.jl.
    "--optimality_norm"
    help =
      "The norm for the optimality criteria. Supported options: {l2, " *
      "l_inf}"
    arg_type = String

    "--absolute_optimality_tol"
    help = "The absolute tolerance for the optimality criteria."
    arg_type = Float64

    "--relative_optimality_tol"
    help = "The relative tolerance for the optimality criteria."
    arg_type = Float64

    # The next two parameters specify infeasibility criteria.
    "--eps_primal_infeasible"
    help = "Tolerance for declaring primal infeasibility."
    arg_type = Float64

    "--eps_dual_infeasible"
    help = "Tolerance for declaring dual infeasibility."
    arg_type = Float64

    # The next two parameters are for early termination.
    "--time_sec_limit"
    help = "Time limit in seconds."
    arg_type = Float64

    "--iteration_limit"
    help = "Maximum number of iterations to run."
    arg_type = Int32

    "--kkt_matrix_pass_limit"
    help = "Terminate after this many passes through the KKT matrix."
    arg_type = Int32

    "--transform_bounds_into_linear_constraints"
    help =
      "Transform all bounds into linear constraints. This is for " *
      "ablation study only. It generally makes the problem harder to " *
      "solve."
    arg_type = Bool
    default = false

    "--fixed_format_input"
    help =
      "If true, parse the input MPS/QPS file in fixed format instead of " *
      "free (the default)."
    arg_type = Bool
    default = false

    # The following parameters define the step size policy.
    "--step_size_policy"
    help =
      "Step size policy used for PDHG. This is ignored for Mirror-prox." *
      " Supported options {constant, adaptive, malitsky-pock}. Defaults" *
      " to 'adaptive'. For the constant step size the solver computes a" *
      " provably correct step size using power iteration."
    arg_type = String
    default = "adaptive"

    "--adaptive_step_size_reduction_exponent"
    help =
      "Adaptive step size rule parameter. New step sizes are" *
      "a factor (1 - iteration^adaptive_step_size_reduction_exponent)" *
      " smaller than they could be as a margin to reduce rejected steps."
    arg_type = Float64
    default = 0.3

    "--adaptive_step_size_growth_exponent"
    help =
      "Adaptive step size rule parameter. New step sizes are at most" *
      "(1+iteration^adaptive_step_size_growth_exponent) * current_step_size."
    arg_type = Float64
    default = 0.6

    "--malitsky_pock_downscaling_factor"
    help =
      "Malitsky and Pock step size parameter. Factor by which the step size " *
      "is multiplied for in the inner loop. Corresponds to mu in the paper " *
      "(https://arxiv.org/pdf/1608.08883.pdf)."
    arg_type = Float64
    default = 0.7

    "--malitsky_pock_breaking_factor"
    help =
      "Malitsky and Pock step size parameter. The breaking factor " *
      "defines the stopping criteria of the linesearch. It should be in the " *
      "interval (0.0, 1.0]. Corresponds to delta in the paper " *
      "(https://arxiv.org/pdf/1608.08883.pdf)."
    arg_type = Float64
    default = 0.99

    "--malitsky_pock_interpolation_coefficient"
    help =
      "Malitsky and Pock step size parameter. Interpolation coefficient " *
      "to pick next step size. The next step size can be picked within an" *
      " interval [a, b] (See Step 2 of Algorithm 1 in " *
      "https://arxiv.org/pdf/1608.08883.pdf). The solver uses " *
      "a + interpolation_coefficient * (b - a)."
    arg_type = Float64
    default = 1.0

  end

  return ArgParse.parse_args(arg_parse)
end

function string_to_restart_scheme(restart_scheme::String)
  if restart_scheme == "no_restart"
    return FirstOrderLp.NO_RESTARTS
  elseif restart_scheme == "adaptive_normalized"
    return FirstOrderLp.ADAPTIVE_NORMALIZED
  elseif restart_scheme == "adaptive_distance"
    return FirstOrderLp.ADAPTIVE_DISTANCE
  elseif restart_scheme == "adaptive_localized"
    return FirstOrderLp.ADAPTIVE_LOCALIZED
  elseif restart_scheme == "fixed_frequency"
    return FirstOrderLp.FIXED_FREQUENCY
  else
    error("Unknown restart scheme $(restart_scheme)")
  end
end

function string_to_restart_to_current_metric(restart_to_current_metric::String)
  if restart_to_current_metric == "no_restart_to_current"
    return FirstOrderLp.NO_RESTART_TO_CURRENT
  elseif restart_to_current_metric == "gap_over_distance"
    return FirstOrderLp.GAP_OVER_DISTANCE
  elseif restart_to_current_metric == "gap_over_distance_squared"
    return FirstOrderLp.GAP_OVER_DISTANCE_SQUARED
  else
    error(
      "Unknown value for restart_to_current_metric $(restart_to_current_metric)",
    )
  end
end

function main()
  parsed_args = parse_command_line()

  if parsed_args["method"] == "mirror-prox" || parsed_args["method"] == "pdhg"
    restart_params = FirstOrderLp.construct_restart_parameters(
      string_to_restart_scheme(parsed_args["restart_scheme"]),
      string_to_restart_to_current_metric(
        parsed_args["restart_to_current_metric"],
      ),
      parsed_args["restart_frequency"],
      parsed_args["artificial_restart_threshold"],
      parsed_args["sufficient_reduction_for_restart"],
      parsed_args["necessary_reduction_for_restart"],
      parsed_args["primal_weight_update_smoothing"],
      parsed_args["use_approximate_localized_duality_gap"],
    )

    pock_chambolle_alpha = nothing
    if parsed_args["pock_chambolle_rescaling"]
      pock_chambolle_alpha = parsed_args["pock_chambolle_alpha"]
    end

    termination_criteria = FirstOrderLp.construct_termination_criteria()
    if parsed_args["optimality_norm"] == "l2"
      termination_criteria.optimality_norm = FirstOrderLp.OptimalityNorm.L2
    elseif parsed_args["optimality_norm"] == "l_inf"
      termination_criteria.optimality_norm = FirstOrderLp.OptimalityNorm.L_INF
    elseif parsed_args["optimality_norm"] !== nothing
      error("Unknown termination norm.")
    end
    for (field_name, arg_name) in [
      (:eps_optimal_absolute, "absolute_optimality_tol"),
      (:eps_optimal_relative, "relative_optimality_tol"),
      (:eps_primal_infeasible, "eps_primal_infeasible"),
      (:eps_dual_infeasible, "eps_dual_infeasible"),
      (:time_sec_limit, "time_sec_limit"),
      (:iteration_limit, "iteration_limit"),
      (:kkt_matrix_pass_limit, "kkt_matrix_pass_limit"),
    ]
      if parsed_args[arg_name] !== nothing
        setproperty!(termination_criteria, field_name, parsed_args[arg_name])
      end
    end

    if parsed_args["method"] == "mirror-prox"
      parameters = FirstOrderLp.MirrorProxParameters(
        parsed_args["l_inf_ruiz_iterations"],
        parsed_args["l2_norm_rescaling"],
        pock_chambolle_alpha,
        parsed_args["primal_importance"],
        parsed_args["scale_invariant_initial_primal_weight"],
        parsed_args["diagonal_scaling"],
        parsed_args["verbosity"],
        parsed_args["record_iteration_stats"],
        parsed_args["termination_evaluation_frequency"],
        termination_criteria,
        restart_params,
      )
    elseif parsed_args["method"] == "pdhg"
      if parsed_args["step_size_policy"] == "malitsky-pock"
        step_size_policy_params = FirstOrderLp.MalitskyPockStepsizeParameters(
          parsed_args["malitsky_pock_downscaling_factor"],
          parsed_args["malitsky_pock_breaking_factor"],
          parsed_args["malitsky_pock_interpolation_coefficient"],
        )
      elseif parsed_args["step_size_policy"] == "constant"
        step_size_policy_params = FirstOrderLp.ConstantStepsizeParams()
      else
        step_size_policy_params = FirstOrderLp.AdaptiveStepsizeParams(
          parsed_args["adaptive_step_size_reduction_exponent"],
          parsed_args["adaptive_step_size_growth_exponent"],
        )
      end
      parameters = FirstOrderLp.PdhgParameters(
        parsed_args["l_inf_ruiz_iterations"],
        parsed_args["l2_norm_rescaling"],
        pock_chambolle_alpha,
        parsed_args["primal_importance"],
        parsed_args["scale_invariant_initial_primal_weight"],
        parsed_args["verbosity"],
        parsed_args["record_iteration_stats"],
        parsed_args["termination_evaluation_frequency"],
        termination_criteria,
        restart_params,
        step_size_policy_params,
      )
    end
  else
    error("`method` arg must be either `mirror-prox` or `pdhg`.")
  end

  solve_instance_and_output(
    parameters,
    parsed_args["output_dir"],
    parsed_args["instance_path"],
    parsed_args["redirect_stdio"],
    parsed_args["transform_bounds_into_linear_constraints"],
    parsed_args["fixed_format_input"],
  )
end

main()


text:
import Base: size, axpy!
import Base.LinAlg: scale!
import Base.BLAS: gemm!
import ArrayViews: view, StridedView, ContiguousView
import Base: shmem_rand, shmem_randn#, acontrank
export fit!, localcols

# functions for shared arrays
function localcols(Y::SharedArray)
    idxs=localindexes(Y)
    s,t=localindexes(Y)[1],localindexes(Y)[end]
    m,n=size(Y)
    return round(floor((s-1)/m+1)):round(floor(t/m))
end

### Fitting

function fit!(glrm::ShareGLRM, params::ProxGradParams;
    ch::ConvergenceHistory=ConvergenceHistory("ProxGradShareGLRM"),verbose=true)

	### initialization (mostly name shortening)
    isa(glrm.A, SharedArray) ? A = glrm.A : A = convert(SharedArray,glrm.A)
    # make sure that we've oriented the factors correctly
    # k, m == size(X) and k, n == size(Y)
    if size(glrm.Y,1)!==size(glrm.X,1)
        glrm.X = glrm.X'
    end
    # check that we didn't initialize to zero (otherwise we will never move)
    if norm(glrm.Y) == 0
        glrm.Y = .1*shmem_randn(size(glrm.Y)...)
    end
    # at any time, glrm.X and glrm.Y will be the best model yet found, while
    # X and Y will be the working variables. All of these are shared arrays.
    X, glrm.X = glrm.X, copy(glrm.X)
    Y, glrm.Y = glrm.Y, copy(glrm.Y)

    ## a few scalars that need to be shared among all processes
    # step size (will be scaled below to ensure it never exceeds 1/\|g\|_2 or so for any subproblem)
    alpha = Base.shmem_fill(float(params.stepsize),(1,1))
    obj = Base.shmem_fill(0.0,(nprocs(),1))

    if verbose println("sending data and caching views") end
    @sync begin
        for i in procs(A)
            remotecall(i, x->(global const A=x; nothing), A)
            remotecall(i, x->(global const X=x; nothing), X)
            remotecall(i, x->(global const Y=x; nothing), Y)
            remotecall(i, x->(global const glrm=x; nothing), glrm)
            remotecall(i, x->(global const alpha=x; nothing), alpha)
            remotecall(i, x->(global const obj=x; nothing), obj)
        end
    end
    @everywhere begin
        # rename data to be easier to access on local proc
        A = LowRankModels.A # since these are shared arrays, this should only be passing the references to them
        X = LowRankModels.X
        Y = LowRankModels.Y
        glrm = LowRankModels.glrm
        k = LowRankModels.glrm.k
        losses = LowRankModels.glrm.losses
        rx = LowRankModels.glrm.rx
        ry = LowRankModels.glrm.ry # this is a normal array that is getting copied by value to all the procs
        of = LowRankModels.glrm.observed_features
        oe = LowRankModels.glrm.observed_examples
        alpha = LowRankModels.alpha
        obj = LowRankModels.obj
        axpy! = Base.BLAS.axpy!
        gemm! = Base.BLAS.gemm!
        prox! = LowRankModels.prox!
        localcols = LowRankModels.localcols
        grad = LowRankModels.grad
        evaluate = LowRankModels.evaluate

        # cache views into X and Y and the indeces of local columns of X and Y
        m,n = size(A)
        ve = LowRankModels.ContiguousView{Float64,1,Array{Float64,2}}[LowRankModels.view(X.s,:,e) for e=1:m]
        vf = LowRankModels.ContiguousView{Float64,1,Array{Float64,2}}[LowRankModels.view(Y.s,:,f) for f=1:n]
        xlcols = localcols(X)
        ylcols = localcols(Y)
        XY_x = Array(Float64,(length(xlcols), n))
        gemm!('T','N', 1.0, X[:,xlcols], Y.s, 0.0, XY_x) # initialize this for the first iteration
        XY_y = Array(Float64,(m, length(ylcols)))

        # initialize gradient
        g = zeros(k)
    end

    # stopping criterion: stop when decrease in objective < tol
    tol = params.abs_tol * mapreduce(length,+,glrm.observed_features)

    # alternating updates of X and Y
    if verbose println("fitting GLRM") end
    update_ch!(ch, 0, objective(glrm,X,Y))
    t = time()
    steps_in_a_row = 0

    for i=1:params.max_iter
        @everywhere begin
            # X update
#            XY_x = X[:,xlcols]' * Y  #(rows of the approximation matrix that this processor is in charge of)
            # this is computed before the first iteration and subsequently in the objective evaluation
            for e=xlcols
                scale!(g, 0)  # reset gradient to 0
                # compute gradient of L with respect to Xᵢ as follows:
                # ∇{Xᵢ}L = Σⱼ dLⱼ(XᵢYⱼ)/dXᵢ
                for f in of[e]
                    # but we have no function dLⱼ/dXᵢ, only dLⱼ/d(XᵢYⱼ) aka dLⱼ/du
                    # by chain rule, the result is: Σⱼ dLⱼ(XᵢYⱼ)/du * Yⱼ, where dLⱼ/du is the our grad() function
                    axpy!(grad(losses[f], XY_x[e-xlcols[1]+1,f], A[e,f]), vf[f], g)
                end
                # take a proximal gradient step
                l = length(of[e]) + 1
                scale!(g, -alpha[1]/l)
                ## gradient step: Xᵢ += -(α/l) * ∇{Xᵢ}L
                axpy!(1, g, ve[e])
                ## prox step: Xᵢ = prox_rx(Xᵢ, α/l)
                prox!(rx, ve[e], alpha[1]/l)
            end
        end
        @everywhere begin
            # Y update
            # XY_y = X'*Y[:,ylcols]  (columns of the approximation matrix that this processor is in charge of)
            gemm!('T','N', 1.0, X.s, Y[:,ylcols], 0.0, XY_y)
            for f=ylcols
                scale!(g, 0) # reset gradient to 0
                # compute gradient of L with respect to Yⱼ as follows:
                # ∇{Yⱼ}L = Σⱼ dLⱼ(XᵢYⱼ)/dYⱼ
                for e in oe[f]
                    # but we have no function dLⱼ/dYⱼ, only dLⱼ/d(XᵢYⱼ) aka dLⱼ/du
                    # by chain rule, the result is: Σⱼ dLⱼ(XᵢYⱼ)/du * Xᵢ, where dLⱼ/du is our grad() function
                    axpy!(grad(losses[f], XY_y[e,f-ylcols[1]+1], A[e,f]), ve[e], g)
                end
                # take a proximal gradient step
                l = length(oe[f]) + 1
                scale!(g, -alpha[1]/l)
                ## gradient step: Yⱼ += -(α/l) * ∇{Yⱼ}L
                axpy!(1, g, vf[f])
                ## prox step: Yⱼ = prox_ryⱼ(Yⱼ, α/l)
                prox!(ry[f], vf[f], alpha[1]/l)
            end
        end
        # evaluate objective, splitting up among processes by columns of X
        @everywhere begin
            pid = myid()
            gemm!('T','N', 1.0, X[:,xlcols], Y.s, 0.0, XY_x) # XY = X[:,xlcols]'*Y
            err = 0
            for e=xlcols
                for f in of[e]
                    err += evaluate(losses[f], XY_x[e-xlcols[1]+1,f], A[e,f])
                end
            end
            # add regularization penalty
            for e=xlcols
                err += evaluate(rx,ve[e])
            end
            for f=ylcols
                err += evaluate(ry[f],vf[f]::AbstractArray)
            end
            obj[pid] = err
        end
        # make sure parallel obj eval is the same as local (it is)
        # println("local objective = $(objective(glrm,X,Y)) while shared objective = $(obj[1])")
        # record the best X and Y yet found
        totalobj = sum(obj)
        if totalobj < ch.objective[end]
            t = time() - t
            update_ch!(ch, t, totalobj)
            #copy!(glrm.X, X); copy!(glrm.Y, Y)
            @everywhere begin
                @inbounds for i in localindexes(X)
                    glrm.X[i] = X[i]
                end
                @inbounds for i in localindexes(Y)
                    glrm.Y[i] = Y[i]
                end
            end
            alpha[1] = alpha[1] * 1.05 # sketchy constant
            steps_in_a_row = max(1, steps_in_a_row+1)
            t = time()
        else
            # if the objective went up, reduce the step size, and undo the step
            alpha[1] = alpha[1] * (1 / max(1.5, -steps_in_a_row)) # another sketchy constant
            println("objective went up to $(obj[1]); changing step size to $(alpha[1])")
            @everywhere begin
                @inbounds for i in localindexes(X)
                    X[i] = glrm.X[i]
                end
                @inbounds for i in localindexes(Y)
                    Y[i] = glrm.Y[i]
                end
            end
            @everywhere begin
                gemm!('T','N', 1.0, glrm.X[:,xlcols], glrm.Y.s, 0.0, XY_x) # reset XY to previous best
            end
            # X[:], Y[:] = copy(glrm.X), copy(glrm.Y)
            steps_in_a_row = min(0, steps_in_a_row-1)
        end
        # check stopping criterion
        if i>10 && (steps_in_a_row > 3 && ch.objective[end-1] - totalobj < tol) || alpha[1] <= params.min_stepsize
            break
        end
        if verbose && i%10==0
            println("Iteration $i: objective value = $(ch.objective[end])")
        end
    end
    t = time() - t
    update_ch!(ch, t, ch.objective[end])

    return glrm.X', glrm.Y, ch
end


text:
module Applications

using FileType

TypeWasm = FileType.Types("wasm", MIME("application/wasm"))
TypeDex  = FileType.Types("dex", MIME("application/vnd.android.dex"))
TypeDey  = FileType.Types("dey", MIME("application/vnd.android.dey"))

# Wasm detects a Web Assembly 1.0 filetype.
function Wasm(buf)::Bool
	# WASM has starts with `\0asm`, followed by the version.
	# http://webassembly.github.io/spec/core/binary/modules.html#binary-magic
	return length(buf) >= 8 &&
		buf[1] == 0x00 && buf[2] == 0x61 &&
		buf[3] == 0x73 && buf[4] == 0x6D &&
		buf[5] == 0x01 && buf[6] == 0x00 &&
		buf[7] == 0x00 && buf[8] == 0x00
end

# Dex detects dalvik executable(DEX)
function Dex(buf)::Bool
	# https://source.android.com/devices/tech/dalvik/dex-format#dex-file-magic
	return length(buf) > 36 &&
		# magic
		buf[1] == 0x64 && buf[2] == 0x65 && buf[3] == 0x78 && buf[4] == 0x0A &&
		# file size
		buf[37] == 0x70
end

# Dey Optimized Dalvik Executable(ODEX)
function Dey(buf)::Bool
	return length(buf) > 100 &&
		# dey magic
		buf[1] == 0x64 && buf[2] == 0x65 && buf[3] == 0x79 && buf[4] == 0x0A &&
		# dex
		Dex(buf[41:101])
end

Application = Dict(
	TypeWasm=>Wasm,
	TypeDex=>Dex,
	TypeDey=>Dey,
)

end


text:
using BinaryBuilder

cuda_version = v"9.2.148"   # NOTE: could be less specific

sources_linux = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/9.2_20191031/cudnn-9.2-linux-x64-v7.6.5.32.tgz",
                  "a2a2c7a8ba7b16d323b651766ee37dcfdbc2b50d920f73f8fde85005424960e4")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/9.2_20191031/cudnn-9.2-windows10-x64-v7.6.5.32.zip",
                  "ffa553df2e9af1703bb7786a784356989dac5c415bf5bca73e52b1789ddd4984")
]

include("../common.jl")


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"

function receive()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare a queue
            queue = "task_queue"
            success, queue_name, message_count, consumer_count = queue_declare(chan, queue, durable=true)
            println(" [*] Waiting for messages. To exit press CTRL+C")

            # 4. Setup function to receive message
            on_receive = (msg) -> begin
                data = String(msg.data)
                println("Received the message: $data")
                basic_ack(chan, msg.delivery_tag)
            end

            # 5. Configure Quality of Service
            basic_qos(chan, 0, 1, false)
            success, consumer_tag = basic_consume(chan, queue, on_receive)

            @assert success == true

            while true
                sleep(1)
            end
            # 5. Close the connection
        end
    end
end

# Don't exit on Ctrl-C
Base.exit_on_sigint(false)
try
    receive()
catch ex
    if ex isa InterruptException
        println("Interrupted")
    else
        println("Exception: $ex")
    end
end

max_forks_repo_name: tomashq/rabbitmq-tutorials


text:
using JuliaDB, GroupedErrors, Sputnik, Images
using Sputnik: process, SelectedData, Analysis
using StatsPlots, Statistics
using Test
school = loadtable(GroupedErrors.exampletable("school.csv"))

function compare_plots(plt1, plt2; sigma = [1,1], eps = 0.02)
    ref1 = joinpath(@__DIR__, "plots", "plot1.png")
    savefig(plt1, ref1)
    ref2 = joinpath(@__DIR__, "plots", "plot2.png")
    savefig(plt2, ref2)
    outcome = Images.test_approx_eq_sigma_eps(Images.load(ref1), Images.load(ref2), sigma, eps)
    rm(ref1)
    rm(ref2)
    outcome
end

@testset "statplots" begin
    data = school
    sd = SelectedData(school)

    a = Analysis(data = sd, x = :MAch, plot = density)
    plt1 = process(a)
    plt2 = @df data density(:MAch)
    @test compare_plots(plt1, plt2) < 0.005

    a = Analysis(data = sd, x = :MAch, y = :SSS, plot = scatter, plot_kwargs = [(:legend, :topleft)])
    plt1 = process(a)
    plt2 = @df data scatter(:MAch, :SSS, legend = :topleft)
    @test compare_plots(plt1, plt2) < 0.005

    sd = SelectedData(school)
    a = Analysis(data = sd, x = :MAch, y = :SSS, plot = scatter)
    plt1 = process(a)
    plt2 = @df school scatter(:MAch, :SSS, legend = false)
    @test compare_plots(plt1, plt2) < 0.005
end

@testset "groupederrors" begin
    sd = SelectedData(school)
    a = Analysis(data = sd,
                 x = :MAch,
                 y = :density,
                 plot = plot,
                 axis_type = :continuous,
                 smoother = 50.0)
    plt1 = process(a)
    bw = (51.0)*std(column(school, :MAch))/200
    plt2 = @> school begin
        @x _.MAch :continuous
        @y :density bandwidth = bw
        @plot plot()
    end
    @test compare_plots(plt1, plt2) < 0.005

    a = Analysis(data = sd,
                 compute_error = (:across, :School),
                 x = :MAch,
                 y = :density,
                 plot = plot,
                 axis_type = :continuous,
                 smoother = 50.0)
    plt1 = process(a)
    plt2 = @> school begin
        @across _.School
        @x _.MAch :continuous
        @y :density bandwidth = bw
        @plot plot()
    end
    @test compare_plots(plt1, plt2) < 0.005

    a = Analysis(data = sd,
                 compute_error = (:across, :School),
                 x = :MAch,
                 y = :SSS,
                 xfunc = mean,
                 yfunc = mean,
                 plot = scatter,
                 axis_type = :pointbypoint)
    plt1 = process(a)
    plt2 = @> school begin
        @across _.School
        @x _.MAch
        @y _.SSS
        @plot scatter()
    end
    @test compare_plots(plt1, plt2) < 0.005
end


text:
# Additional tests of `find_registry_path` and `check_registry`. Many
# of these have the purpose to cover error cases, making them somewhat
# contrived. Another complicating factor is that some of the call
# variants have to interact with the package environment, including
# registries, of the running Julia process.

# Not a registered package.
pkg = Pkg.Types.Project(Dict("name" => "UUIDs",
                             "uuid" => "cf7118a7-6976-5b1a-9a39-7adc72f591a4"))
@test_throws ErrorException find_registry_path(nothing, pkg)

with_empty_registry() do registry_dir, packages_dir
    Pkg.Registry.add(RegistrySpec(path = registry_dir))

    # Use Multibreak as Guinea pig. The sleep is a Travis workaround. See
    # a later comment. This also tests automatically choosing the only
    # installed registry for a new package.
    sleep(1)
    register("Multibreak", push = false, gitconfig = TEST_GITCONFIG)

    # Directory already exists. Also tests code handling a trailing slash.
    create_registry("TestRegistry2", "", gitconfig = TEST_GITCONFIG,
                    push = false)
    @test_throws ErrorException create_registry("TestRegistry2/", "",
                                                gitconfig = TEST_GITCONFIG,
                                                push = false)

    # Find a registry by name.
    package_path = find_package_path("Multibreak")
    pkg = Pkg.Types.read_project(joinpath(package_path, "Project.toml"))
    @test find_registry_path("TestRegistry") == joinpath(first(DEPOT_PATH),
                                                         "registries",
                                                         "TestRegistry")

    # The named registry does not exist. In this case it is assumed to
    # be a local path if it exists, otherwise a URL, regardless how
    # much it doesn't look like a URL.
    saved_dir = pwd()
    tempdir = mktempdir()
    cd(tempdir)
    @test find_registry_path("General", pkg) == "General"
    mkdir("General")
    @test find_registry_path("General", pkg) == abspath("General")
    cd(saved_dir)
    rm(tempdir, recursive = true)

    # Find which registry contains a package.
    @test find_registry_path(nothing, pkg) == joinpath(first(DEPOT_PATH),
                                                       "registries",
                                                       "TestRegistry")

    # Workaround for bad `mtime` resolution of 1 second on MacOS workers
    # on Travis.
    #
    # The issue is that `read_registry` caches its results with respect to
    # the file `mtime`. Since `read_registry` is called from within
    # `register`, the old data will be read into the cache. If the new
    # data is written close enough to the previous registry update so that
    # `mtime` does not change, subsequent `read_registry` will keep using
    # the old data from the cache.
    sleep(1)

    # More than one registry contains the package.
    register("Multibreak", registry = "TestRegistry2",
             repo = "file://$(packages_dir)/FirstTest",
             gitconfig = TEST_GITCONFIG, push = false)
    @test_throws ErrorException find_registry_path(nothing, pkg)
end


text:
# In this file we define the interface for a number of generic functions:
# See the individual files for details on the interfaces.

include("transform.jl")
include("evaluation.jl")
include("approximation.jl")
include("differentiation.jl")

#####################################
# Operators for tensor product dictionaries
#####################################

# We make TensorProductOperator's for each generic operator, when invoked with
# TensorProductDict's.

transform_from_grid(T, s1::GridBasis, s2::TensorProductDict, grid::ProductGrid; options...) =
    tensorproduct(map( (u,v,w) -> transform_from_grid(T, u,v,w; options...), components(s1), components(s2), components(grid))...)

transform_to_grid(T, s1::TensorProductDict, s2::GridBasis, grid::ProductGrid; options...) =
    tensorproduct(map( (u,v,w) -> transform_to_grid(T, u,v,w; options...), components(s1), components(s2), components(grid))...)

for op in (:extension, :restriction, :conversion)
    @eval $op(::Type{T}, src::TensorProductDict, dest::TensorProductDict; options...) where {T} =
        tensorproduct(map( (u,v) -> $op(T, u, v; options...), components(src), components(dest))...)
end

for op in (:interpolation, :leastsquares)
    @eval $op(::Type{T}, s1::TensorProductDict, s2::TensorProductDict; options...) where {T} =
        tensorproduct(map( (u,v) -> $op(T, u, v; options...), components(s1), components(s2))...)
end

dense_evaluation(::Type{T}, s1::TensorProductDict, s2::TensorProductDict; options...) where {T} =
    tensorproduct(map( (u,v) -> dense_evaluation(T, u, v; options...), components(s1), components(s2))...)

for op in (:approximation, )
    @eval $op(s::TensorProductDict; options...) =
        tensorproduct(map( u -> $op(u; options...), components(s))...)
end


text:
import MLBase
import StatsBase

"""
"""
function avg_precision(
        allprecisions::AbstractVector{<:Real},
        allrecalls::AbstractVector{<:Real},
        allthresholds::AbstractVector{<:Real},
        )
    if length(allprecisions) != length(allrecalls)
        error("length(allprecisions) != length(allrecalls)")
    end
    if length(allprecisions) < 2
        error("length(allprecisions) < 2")
    end
    #
    permutation = sortperm(allthresholds; rev = true)
    allprecisions = allprecisions[permutation]
    allrecalls = allrecalls[permutation]
    #
    N = length(allprecisions)
    result = 0
    for k = 2:N
        result += (allrecalls[k] - allrecalls[k-1]) * allprecisions[k]
    end
    return result
end

"""
"""
function averageprecisionscore(
        ytrue::AbstractVector{<:Integer},
        yscore::AbstractVector{<:Real},
        )
    allprecisions, allrecalls, allthresholds = prcurve(
        ytrue,
        yscore,
        )
    x = allrecalls
    y = allprecisions
    avgprecision = avg_precision(
        allprecisions,
        allrecalls,
        allthresholds,
        )
    return avgprecision
end


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0,1,0]
dims = [0,0,0]
MPI.Dims_create!(nnodes, dims)
comm_cart = MPI.Cart_create(comm, dims, periods, reorder)
coords = MPI.Cart_coords!(comm_cart, ndims)
disp = 1

for i in 0:2
    neighbors = MPI.Cart_shift(comm_cart, i, disp)
    @test all( ((neighbors .>= 0) .& (neighbors .< nnodes)) .| (neighbors .== MPI.MPI_PROC_NULL) )
end

MPI.Finalize()


text:
#####
#Uploaders

export put_file, put_jld

"""
Upload a file to the swift object store.
"""
function put_file(serv::PyObject,
                  container::AbstractString,
                  name::AbstractString,
                  fname::AbstractString=name;
                  verbose::Bool=false)
    uo = service[:SwiftUploadObject](fname, object_name=name)
    async_put = serv[:upload](container, [uo])
    responses = collect(async_put)
    if verbose
        show_responses(responses)
    end
    ret = responses[end]
    ret["success"] || error(ret["error"])
    ret
end

"""Note: this reads out the `fp` IO to the end"""
function put_file(serv::PyObject,
                  container::AbstractString,
                  name::AbstractString,
                  fp::IO;
                  verbose::Bool=false)
    mktempdir() do tdir
        fname = joinpath(tdir, "lastswiftupload")
        open(fname,"w") do fp_inner
            write(fp_inner, read(fp))
        end
        put_file(serv, container, name, fname; verbose=verbose)
    end
end

"""Save data as a a JLD file, and upload to Swift"""
function put_jld(serv::PyObject,
                 container::AbstractString,
                 name::AbstractString;
                 data...)
    mktempdir() do tdir
        fname = joinpath(tdir, "lastswiftupload.jl")
        save(File(format"JLD", fname), Base.Flatten(((string(name), val) for (name,val) in data))...)
        put_file(serv, container, name, fname)
    end
end


text:
import Hiccup, Media, Images, Juno, FileIO, ModernGL, Interact

function scene2image(screen::Screen)
    GLWindow.poll_glfw()
    yield()
    render_frame(screen) # let it render
    ModernGL.glFinish()
    return Images.clamp01nan.(GLWindow.screenbuffer(screen))
end

function scene2image(scene::Scene)
    screen = getscreen(scene)
    if screen != nothing
        return scene2image(screen)
    else
        return nothing
    end
end

Base.mimewritable(::MIME"image/png", scene::Scene) = true

function Base.show(io::IO, mime::MIME"image/png", scene::Makie.Scene)
    s = to_signal(lift_node(scene, :entered_window) do value
        scene2image(scene)
    end)
    display(s)
end

function save(path::String, scene::Scene)
    img = scene2image(scene)
    if img != nothing
        FileIO.save(path, img)
    else
        error("Scene doesn't contain a plot!")
    end
end

Media.media(Scene, Media.Plot)

const use_atom_plot_pane = Ref(false)
use_plot_pane(x::Bool = true) = (use_atom_plot_pane[] = x)

function Juno.render(pane::Juno.PlotPane, plt::Scene)
    if use_atom_plot_pane[]
        img = scene2image(plt)
        if img != nothing
            Juno.render(pane, HTML("<img src=\"data:image/png;base64,$(stringmime(MIME("image/png"), img))\">"))
        end
    end
end

immutable VideoStream
    io
    buffer
    screen::GLWindow.Screen
    path::String
end

"""
    VideoStream(scene::Scene, dir = mktempdir(), name = "video")

returns a stream and a buffer that you can use to not allocate for new frames.
Use `add_frame!(stream, window, buffer)` to add new video frames to the stream.
Use `finish(stream)` to save the video to 'dir/name.mkv'. You can also call
`finish(stream, "mkv")`, `finish(stream, "mp4")`, `finish(stream, "gif")` or `finish(stream, "webm")` to convert the stream to those formats.
"""
function VideoStream(scene::Scene, dir = mktempdir(), name = "video")
    #codec = `-codec:v libvpx -quality good -cpu-used 0 -b:v 500k -qmin 10 -qmax 42 -maxrate 500k -bufsize 1000k -threads 8`
    path = joinpath(dir, "$name.mkv")
    screen = rootscreen(scene)
    if screen == nothing
        error("Scene doesn't contain any screen")
    end
    tex = GLWindow.framebuffer(screen).color
    res = size(tex)
    io, process = open(`ffmpeg -loglevel panic -f rawvideo -pixel_format rgb24 -s:v $(res[1])x$(res[2]) -i pipe:0 -vf vflip -y $path`, "w")
    VideoStream(io, Matrix{RGB{N0f8}}(res), screen, abspath(path)) # tmp buffer
end

"""
Adds a video frame to the VideoStream
"""
function recordframe!(io::VideoStream)
    #codec = `-codec:v libvpx -quality good -cpu-used 0 -b:v 500k -qmin 10 -qmax 42 -maxrate 500k -bufsize 1000k -threads 8`
    GLWindow.poll_glfw()
    yield()
    render_frame(io.screen) # let it render
    ModernGL.glFinish()
    tex = GLWindow.framebuffer(io.screen).color
    write(io.io, map(RGB{N0f8}, gpu_data(tex)))
    return
end

"""
    finish(io::VideoStream, typ = "mkv"; remove_mkv = true)

Flushes the video stream and optionally converts the file to `typ` which can
be (`mkv` is default and doesn't need convert) gif, mp4 and webm.
If you want to convert the original mkv to multiple formats you should choose
`remove_mkv = false`, and remove it manually after you're done (with `rm(videostream.path)`)
webm yields the smallest file size, mp4 and mk4 are marginally bigger and gifs are up to
6 times bigger!
"""
function finish(io::VideoStream, typ = "mkv"; remove_mkv = true)
    close(io.io)
    path = dirname(io.path)
    typ == "mkv" && return io.path
    name, ext = splitext(basename(io.path))
    path = dirname(io.path)
    out = joinpath(path, name * ".$typ")

    if typ == "mp4"
        run(`ffmpeg -loglevel quiet -i $(io.path) -vcodec copy -acodec copy -y $out`)
        remove_mkv && rm(io.path)
        return out
    elseif typ == "webm"
        run(`ffmpeg -loglevel quiet -i $(io.path) -c:v libvpx-vp9 -threads 16 -b:v 2000k -c:a libvorbis -threads 16 -vf scale=iw:ih -y $out`)
        remove_mkv && rm(io.path)
        return out
    elseif typ == "gif"
        palette_path = mktempdir()
        pname = joinpath(palette_path, "palette.png")
        filters = "fps=15,scale=iw:ih:flags=lanczos"
        run(`ffmpeg -loglevel quiet -v warning -i $(io.path) -vf "$filters,palettegen" -y $pname`)
        run(`ffmpeg -loglevel quiet -v warning -i $(io.path) -i $pname -lavfi "$filters [x]; [x][1:v] paletteuse" -y $out`)
        remove_mkv && rm(io.path)
        return out
    else
        error("Video type $typ not known")
    end
end
Base.mimewritable(::MIME"text/html", scene::VideoStream) = true

function Base.show(io::IO, mime::MIME"text/html", vs::VideoStream)
    path = finish(vs, "mp4")
    html = """
        <video width="100%" autoplay controls="false">
            <source src="$path" type="video/mp4">
            Your browser does not support the video tag. Please use a modern browser like Chrome or Firefox.
        </video>
    """
    display(
        Main.IJulia.InlineDisplay(),
        mime,
        string(
            """<video autoplay controls><source src="data:video/x-m4v;base64,""",
            base64encode(open(read, path)),"""" type="video/mp4"></video>"""
        )
    )
end

# function Base.show(io::IO, mime::MIME"text/html", vs::VideoStream)
#     path = "file/" * finish(vs, "mp4")
#     html = """
#         <video width="100%" autoplay controls="false">
#             <source src="%path" type="video/mp4">
#             Your browser does not support the video tag. Please use a modern browser like Chrome or Firefox.
#         </video>
#     """
#     display(
#         Main.IJulia.InlineDisplay(),
#         mime,
#         html
#     )
# end

export VideoStream, recordframe!, finish

#
# function to_gif()
#     ffmpeg -y -i ijulia.ogv -vf fps=20, scale=320:-1:flags=lanczos,palettegen palette.png
#     ffmpeg -v warning -i ijulia.ogv -vf scale=300:-1 -gifflags +transdiff -y bbb-trans.gif
#


text:
""" Rewriting for GAT expressions.

The current content of this module is just a stopgap until I can implement
a generic term rewriting system.
"""
module Rewrite
export associate, associate_unit, distribute_unary, anti_involute

using ..Syntax

""" Simplify associative binary operation.

Maintains the normal form `op(e1,e2,...)` where `e1`,`e2`,... are expressions
that are *not* applications of `op()`
"""
function associate(expr::E)::E where E <: GATExpr
  op, e1, e2 = head(expr), first(expr), last(expr)
  args1 = head(e1) == op ? args(e1) : [e1]
  args2 = head(e2) == op ? args(e2) : [e2]
  E([args1; args2], type_args(expr))
end

""" Simplify associative binary operation with unit.

Reduces a freely generated (typed) monoid to normal form.
"""
function associate_unit(expr::GATExpr, unit::Function)::GATExpr
  e1, e2 = first(expr), last(expr)
  if (head(e1) == nameof(unit)) e2
  elseif (head(e2) == nameof(unit)) e1
  else associate(expr) end
end

""" Distribute unary operation over a binary operation.
"""
function distribute_unary(raw_expr::GATExpr, un_op::Function,
                          bin_op::Function)::GATExpr
  if head(raw_expr) != nameof(un_op)
    return raw_expr
  end
  expr = first(raw_expr)
  if head(expr) == nameof(bin_op)
    bin_op([un_op(A) for A in args(expr)]...)
  else
    raw_expr
  end
end

""" Simplify unary operation that is an anti-involution on a (typed) monoid.
"""
function anti_involute(raw_expr::GATExpr, inv::Function, op::Function,
                       unit::Function)::GATExpr
  if head(raw_expr) != nameof(inv)
    return raw_expr
  end
  expr = first(raw_expr)
  if head(expr) == nameof(inv)
    first(expr)
  elseif head(expr) == nameof(op)
    op([inv(A) for A in reverse(args(expr))]...)
  elseif head(expr) == nameof(unit)
    expr
  else raw_expr end
end

end


text:
export pred_ared_computation

function pred_ared_computation(gₖ, fₖ, dersec, d, ftestTR, gtestTR;
                               seuil = -1e-10)
  pred = gₖ .* d + 0.5 .* dersec .* d .^ 2

  if true in (pred .> seuil)
    ared = (gₖ .+ gtestTR) .* d ./ 2.0
  else
    ared = ftestTR .- fₖ
  end

  ratio = ared ./ pred

  return (pred[1], ared[1], ratio[1])
end


text:
using Test
using OmegaCore
using Distributions
using OmegaTest
using OmegaCore.Interventions

function test_mergetags()
  x = 1 ~ Normal(0.0, 1.0)
  y(ω) = x(ω) + 10.0
  c1 = (ω -> 100.0)
  c2 = (ω -> 200.0)
  yi = intervene(y, x => c1)
  yii = intervene(yi, x => c2)
  yiii = intervene(yi, (x => c1, x => c2))
  yiii2 = intervene(yi, (x => c2, x => c1))

  nt1 = (intervene = yi.i,)
  nt2 = (intervene = yii.i,)
  ntmerged = mergetags(nt1, nt2)
  @test isinferred(mergetags, nt1, nt2)

  # Only one of these should be true
  # @test ntmerged.intervene == yiii.i
  @test ntmerged.intervene == yiii2.i
end

function test_changed_rettype_merge()
  xx = 1 ~ Categorical([0.5, 0.5])
  y(ω) = xx(ω) + 1 # int
  yi = intervene(y, xx => (ω -> 200.0)) # float
  @test randsample(yi) == 201
  ω = defω()
  @test Base.return_types(yi, Base.typesof(ω))[1] == Union{Int64, Float64}
end

function test_merge_1()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  yi = intervene(y, xx => (ω -> 200.0))
  yi2 = intervene(yi, xx => (ω -> 300.0))
  @test randsample(yi2) == 210
  # yi2(ω)
  @test isinferred(randsample, yi2)
end

function test_merge_2()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  xr = 2 ~ Normal(30, 1)
  yi = intervene(y, xx => xr)
  yi2 = intervene(yi, xr => (ω -> 300.0))
  @test randsample(yi2) == 310
  @test isinferred(randsample, yi2)
end

function test_merge_3()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  yi = intervene(y, (xx => (ω -> 200.0), xx => (ω -> 300.0)))
  yi3 = intervene(yi, xx => (ω -> 400.0))
  @test randsample(yi3) == 210
  @test isinferred(randsample, yi3)
end

function test_merge_4()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  xr = 2 ~ Normal(30, 1)
  xrr(ω) = xr(ω) * xr(ω)
  yi = intervene(y, xx => xrr)
  yi2 = intervene(yi, xr => (ω -> 30.0))
  @test randsample(yi2) == 910.0
  @test isinferred(randsample, yi2)
end

function test_merge_more_than_5_interventions()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  yi = intervene(y, xx => (ω -> 200.0))
  yi2 = intervene(yi, xx => (ω -> 300.0))
  yi3 = intervene(yi2, xx => (ω -> 400.0))
  yi4 = intervene(yi3, xx => (ω -> 500.0))
  yi5 = intervene(yi4, xx => (ω -> 600.0))
  yi6 = intervene(yi5, xx => (ω -> 700.0))
  @test randsample(yi6) == 210
  # isinferred fails
end

function minimal_example()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  yi = intervene(y, xx => (ω -> 200.0))
  yi2 = intervene(yi, xx => (ω -> 300.0))
  @test randsample(yi2) == 200
end

function minimal_example()
  xx = 1 ~ Normal(0, 1)
  y(ω) = xx(ω) + 10
  xr = 2 ~ Normal(30, 1)
  yi = intervene(y, xx => xr)
  yi2 = intervene(yi, xr => (ω -> 300.0))
  @test randsample(yi2) == 310
end

function test_model()
  # Normally distributed random variable with id 1
  x = 1 ~ Normal(0, 1)

    # Normally distributed random variable with id 2 and x as mean
    function y(ω)
      xt = x(ω)
      u = Uniform(xt, xt + 1)
      u((2,), ω)
    end
    x_ = 0.1
    y_ = 0.3

    # An omega object -- like a trace
    ω = SimpleΩ(Dict((1,) => x_, (2,) => y_))

    # model -- tuple-valued random variable ω -> (x(ω), y(ω)), becase we want joint pdf
    m(ω) = (x(ω), y(ω))
    (x, y, m)
end

function samplemodel()
  x,y,z = test_model()
  ω = defω()
  y(ω)
end

function test_intervention()
  x, y, m = test_model()
  yⁱ = y |ᵈ (x => (ω -> 100.0))
  @test 100.0 <= randsample(yⁱ) <= 101.0
  @test isinferred(randsample, yⁱ)
end

function test_intervene_diff_parents()
  x = 1 ~ Normal(0, 1)
  function y(ω)
    x_ = x(ω)
    (2 ~ Normal(x_, 1))(ω)
  end
  x2 = 3 ~ Normal(0, 1)
  yi = y |ᵈ (x => ω -> 100.0)
  yi2 = y |ᵈ (x2 => ω -> 100.0)
  yi_, yi2_ = randsample((yi, yi2))
  @test yi_ != yi2_
end

function test_two_interventions()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Uniform(10.0, 20.0)
  z(ω) = Normal(x(ω), y(ω))((3,), ω)
  (x, y, z)
  zi = z |ᵈ (x => (ω -> 100.0), y => (ω -> 0.1))
  @test 99 <= randsample(zi) <= 101
  @test isinferred(randsample, zi)
end

function test_three_interventions()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Uniform(10.0, 20.0)
  c = 3 ~ Uniform(2.0, 3.0)
  z(ω) = Normal(x(ω)*c(ω), y(ω))((3,), ω)
  (x, y, z)
  zi = z |ᵈ (x => (ω -> 100.0), y => (ω -> 0.1), c => (w -> 1.0))
  @test 99 <= randsample(zi) <= 101
end

function test_intervention_logpdf()
  # Log density of model wrt ω
  l = logpdf(m, ω)

  # Check it is what it should be
  @test l == logpdf(Normal(0, 1), x_) + logpdf(Normal(x_, 1), y_)

  # Intervened model
  v_ = 100.0

  # y had x been v_
  yⁱ = y | had(x => v_)

  # new model with Intervened variables
  mⁱ = rt(x, yⁱ)

  # log pdf of interved model on same ω
  lⁱ = logpdf(mⁱ, ω)
  logpdf(x, ω)

  @test lⁱ == logpdf(Normal(0, 1), x_) + logpdf(Normal(v_, 1), y_)
  @test lⁱ < l
end

@testset "intervene" begin
  #test_intervention()
  #test_intervene_diff_parents()
  test_two_interventions()
  test_three_interventions()
  # test_intervention_logpdf()
  test_mergetags()
  test_merge_1()
  test_merge_2()
  test_merge_3()
  test_merge_4()
  test_changed_rettype_merge()
  test_merge_more_than_5_interventions()
end


text:
function test_memory_output_layer(backend::Backend, T, eps)
  println("-- Testing Memory Output Layer on $(typeof(backend)){$T}...")

  tensor_dim = abs(rand(Int)) % 4 + 2
  dims = tuple((abs(rand(Int, tensor_dim)) % 8 + 1)...)
  println("    > $dims")

  inputs = Array[rand(T, dims), rand(T, dims)]
  input_blobs = Blob[make_blob(backend, x) for x in inputs]

  layer = MemoryOutputLayer(bottoms=[:input1, :input2])
  state = setup(backend, layer, input_blobs, Array(Blob, length(inputs)))

  # repeat 2 times
  forward(backend, state, input_blobs)
  forward(backend, state, input_blobs)

  @test length(state.outputs) == length(inputs)
  for i = 1:length(inputs)
    @test length(state.outputs[i]) == 2
    @test all(abs(state.outputs[i][1]-inputs[i]) .< eps)
    @test all(abs(state.outputs[i][2]-inputs[i]) .< eps)
  end

  reset_outputs(state)
  forward(backend, state, input_blobs)
  @test length(state.outputs) == length(inputs)
  for i = 1:length(inputs)
    @test length(state.outputs[i]) == 1
    @test all(abs(state.outputs[i][1]-inputs[i]) .< eps)
  end

  shutdown(backend, state)
end

function test_memory_output_layer(backend::Backend)
  test_memory_output_layer(backend, Float32, 1e-5)
  test_memory_output_layer(backend, Float64, 1e-10)
end

if test_cpu
  test_memory_output_layer(backend_cpu)
end
if test_gpu
  test_memory_output_layer(backend_gpu)
end


text:
export drawLine

# drawLine using Bresenham's Line algorithm (3-D implementation)
"""
    drawLine(p1::Tuple{Real,Real,Real}, p2::Tuple{Real,Real,Real}, block::PiCraft.Block = PiCraft.GOLD_BLOCK, width::Real = 1)

Draw a line from point `p1` to point `p2` using `block`.
`block` is gold block by default.

"""
function drawLine(p1::Tuple{Real,Real,Real}, p2::Tuple{Real,Real,Real}, block::PiCraft.Block = PiCraft.GOLD_BLOCK, width::Real = 1.0)
    p = p1 = round.(Int, p1)
    p2 = round.(Int, p2)

    if(width < 1)
        error("Width must be greater than or equal to 1")
    end
    if width > 1.0
       r = (-width/2):(width/2)
        for a = r, b = r, c = r
            if sum((a, b, c).^2) <= (width/2)^2
                setBlock(p .+ (a, b, c), block)
            end
        end
    else
        setBlock(p, block)
    end

    Δx, Δy, Δz = p2 .- p1
    δx, δy, δz = ifelse.((Δx, Δy, Δz) .< 0, -1, 1)
    l, m, n = abs.((Δx, Δy, Δz))
    Δx2, Δy2, Δz2 = abs.(2 .* (Δx, Δy, Δz))
    if (l >= m && l >= n)
        ϵ1, ϵ2 = (Δy2, Δz2) .- (l, l)
        for i in 1:l
            if width > 1.0
               r = (round(-width/2 - 1)):(round(width/2 + 1))
                for a = r, b = r, c = r
                    if sum((a, b, c).^2) <= (width/2)^2
                        setBlock(p .+ (a, b, c), block)
                    end
                end
            else
                setBlock(p, block)
            end
            if ϵ1 > 0
                p = p .+ (0, δy, 0)
                ϵ1 -= Δx2
            end
            if ϵ2 > 0
                p = p .+ (0, 0, δz)
                ϵ2 -= Δx2
            end
            ϵ1 += Δy2
            ϵ2 += Δz2
            p = p .+ (δx, 0, 0)
        end
    elseif m >= l && m >= n
        ϵ1, ϵ2 = (Δx2, Δz2) .- (m, m)
        for i in 1:m
            if width > 1.0
               r = (-width/2):(width/2)
                for a = r, b = r, c = r
                    if sum((a, b, c).^2) <= (width/2)^2
                        setBlock(p .+ (a, b, c), block)
                    end
                end
            else
                setBlock(p, block)
            end
            if ϵ1 > 0
                p = p .+ (δx, 0, 0)
                ϵ1 -= Δy2
            end
            if ϵ2 > 0
                p = p .+ (0, 0, δz)
                ϵ2 -= Δy2
            end
            ϵ1 += Δx2
            ϵ2 += Δz2
            p = p .+ (0, δy, 0)
        end
    else
        ϵ1 , ϵ2 = (Δy2, Δx2) .- (n, n)
        for i in 1:n
            if width > 1.0
               r = (-width/2):(width/2)
                for a = r, b = r, c = r
                    if sum((a, b, c).^2) <= (width/2)^2
                        setBlock(p .+ (a, b, c), block)
                    end
                end
            else
                setBlock(p, block)
            end
            if ϵ1 > 0
                p = p .+ (0, δy, 0)
                ϵ1 -= Δz2
            end
            if ϵ2 > 0
                p = p .+ (δx, 0, 0)
                ϵ2 -= Δz2
            end
            ϵ1 += Δy2
            ϵ2 += Δx2
            p = p .+ (0, 0, δz)
        end
    end
    if width > 1.0
       r = (-width/2):(width/2)
        for a = r, b = r, c = r
            if sum((a, b, c).^2) <= (width/2)^2
                setBlock(p2 .+ (a, b, c), block)
            end
        end
    else
        setBlock(p2, block)
    end
    return
end


text:
using Neighborhood: Theiler, KDTree, BruteForce, bulkisearch, searchstructure
using Distances: Metric, Euclidean
export NaiveKernel, KDTree, BruteForce

"""
    NaiveKernel(ϵ::Real, ss = KDTree; w = 0, metric = Euclidean()) <: ProbabilitiesEstimator

Estimate probabilities/entropy using a "naive" kernel density estimation approach (KDE), as
discussed in Prichard and Theiler (1995) [^PrichardTheiler1995].

Probabilities ``P(\\mathbf{x}, \\epsilon)`` are assigned to every point ``\\mathbf{x}`` by
counting how many other points occupy the space spanned by
a hypersphere of radius `ϵ` around ``\\mathbf{x}``, according to:

```math
P_i( X, \\epsilon) \\approx \\dfrac{1}{N} \\sum_{s} B(||X_i - X_j|| < \\epsilon),
```

where ``B`` gives 1 if the argument is `true`. Probabilities are then normalized.

The search structure `ss` is any search structure supported by Neighborhood.jl.
Specifically, use `KDTree` to use a tree-based neighbor search, or `BruteForce` for
the direct distances between all points. KDTrees heavily outperform direct distances
when the dimensionality of the data is much smaller than the data length.

The keyword `w` stands for the [Theiler window](@ref), and excludes indices ``s``
that are within ``|i - s| ≤ w`` from the given point ``X_i``.

[^PrichardTheiler1995]: Prichard, D., & Theiler, J. (1995). Generalized redundancies for time series analysis. Physica D: Nonlinear Phenomena, 84(3-4), 476-493.
"""
struct NaiveKernel{KM, M <: Metric} <: ProbabilitiesEstimator
    ϵ::Float64
    method::KM
    w::Int
    metric::M
end
function NaiveKernel(ϵ::Real, method = KDTree; w = 0, metric = Euclidean())
    ϵ > 0 || error("Radius ϵ must be larger than zero!")
    return NaiveKernel(ϵ, method, w, metric)
end

function probabilities(x::DelayEmbeddings.AbstractDataset, est::NaiveKernel)
    theiler = Theiler(est.w)
    ss = searchstructure(est.method, x.data, est.metric)
    idxs = bulkisearch(ss, x.data, WithinRange(est.ϵ), theiler)
    p = Float64.(length.(idxs))
    return Probabilities(p)
end


text:
"""
	stationarydist(c::MarkovChain)

Phân phối dừng của xích markov
"""
function stationarydist(c::MarkovChain)
	P = m.P
	n = size(P, 1)
	A = [(transpose(p) - I)[1:end-1, :]; ones(1, n)]
	B = [zeros(n - 1); 1]
	inv(A) * B
end

function Base.show(io::IO, m::MarkovChain)
	display("Xích Markov với ma trận chuyển")
	display(m.P)
end


text:
using HTTP

@testset "HTTP.download" begin
    @testset "Content-Disposition" begin
        invalid_content_disposition_fn = HTTP.download(
            "http://test.greenbytes.de/tech/tc2231/attonlyquoted.asis")
        @test isfile(invalid_content_disposition_fn)
        @test basename(invalid_content_disposition_fn) == "attonlyquoted.asis" # just last part  of name

        content_disposition_fn = HTTP.download(
            "http://test.greenbytes.de/tech/tc2231/inlwithasciifilenamepdf.asis")
        @test isfile(content_disposition_fn)
        @test basename(content_disposition_fn) == "foo.pdf"

        if Sys.isunix() # Don't try this on windows, quotes are not allowed in windows filenames.
            escaped_content_disposition_fn = HTTP.download(
                "http://test.greenbytes.de/tech/tc2231/attwithasciifnescapedquote.asis")
            @test isfile(escaped_content_disposition_fn)
            @test basename(escaped_content_disposition_fn) == "\"quoting\" tested.html"
        end
    end

    @testset "Provided Filename" begin
        provided_filename = tempname()
        returned_filename = HTTP.download(
            "http://test.greenbytes.de/tech/tc2231/inlwithasciifilenamepdf.asis",
            provided_filename
        )
        @test provided_filename == returned_filename
        @test isfile(provided_filename)

    end

    @testset "Content-Encoding" begin
        gzip_content_encoding_fn = HTTP.download("https://httpbin.org/gzip")
        @test isfile(gzip_content_encoding_fn)
        @test last(splitext(gzip_content_encoding_fn)) == ".gz"
    end
end


text:
using Plotly
using JSON
using Test

#test signin only one endpoint specified
@test_throws ErrorException Plotly.signin("fake_username", "fake_api_key", Dict("plotly_domain"=>"test"))

#test signin and get_credentials
Plotly.signin("test_username", "test_api_key")
creds = Plotly.get_credentials()
@test creds.username == "test_username"
@test creds.api_key == "test_api_key"

#test signin + endpoints and get_config
endpoints = Dict("plotly_domain"=>"my_plotly_domain", "plotly_api_domain"=>"my_plotly_api_domain")
Plotly.signin("test_username", "test_api_key", endpoints)
config = Plotly.get_config()
@test config.plotly_domain == "my_plotly_domain"
@test config.plotly_api_domain== "my_plotly_api_domain"