text:
libcublas=dlopen("libcublas")

# Matrix multiply
for (fname, elty) in ((:cublasSgemm, :Float32),
                      (:cublasDgemm, :Float64),
                      (:cublasCgemm, :Complex64),
                      (:cublasZgemm, :Complex128))
    @eval begin
        function cuda_gemm(transA::Char, transB::Char,
                           m::Int32, n::Int32, k::Int32,
                           alpha::($elty), A::Ptr{$elty}, lda::Int32,
                           B::Ptr{$elty}, ldb::Int32, beta::($elty),
                           C::Ptr{$elty}, ldc::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Int32, Int32, Int32,
                         $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32,
                         $elty, Ptr{$elty}, Int32),
                  transA, transB, m, n, k, alpha, A, lda, B, ldb,
                  beta, C, ldc)
        end
    end
end

# Find the index of the absolute maximum value
for (fname, elty) in ((:cublasIdamax, :Float64),
                      (:cublasIsamax, :Float32),
                      (:cublasIcamax, :Complex64),
                      (:cublasIzamax, :Complex128))
    @eval begin
        function cuda_amax(num::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Int32, (Int32, Ptr{$elty}, Int32),
                  num, x, 1)
        end
    end
end

# Find the index of the absolute minimum value
for (fname, elty) in ((:cublasIdamin, :Float64),
                      (:cublasIsamin, :Float32),
                      (:cublasIcamin, :Complex64),
                      (:cublasIzamin, :Complex128))
    @eval begin
        function cuda_amin(num::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Int32, (Int32, Ptr{$elty}, Int32),
                  num, x, 1)
        end
    end
end

# Find the sum of absolute values
for (fname, elty) in ((:cublasDasum, :Float64),
                      (:cublasSasum, :Float32),
                      (:cublasScasum, :Complex64),
                      (:cublasDzasum, :Complex128))
    @eval begin
        function cuda_asum(num::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  $elty, (Int32, Ptr{$elty}, Int32),
                  num, x, 1)
        end
    end
end

# Scale the matrix
for (fname, elty1, elty2) in ((:cublasSscal, :Float32, :Float32),
                              (:cublasDscal, :Float64, :Float64),
                              (:cublasCscal, :Complex64, :Complex64),
                              (:cublasCsscal, :Float32, :Complex64),
                              (:cublasZscal, :Complex128, :Complex128),
                              (:cublasZdscal, :Float64, :Complex128))
    @eval begin
        function cuda_scal(num::Int32, x::Ptr{$elty2}, alpha::($elty1))
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, $elty1, Ptr{$elty2}, Int32),
                  num, alpha, x, 1)
        end
    end
end

# Dot product
for (fname, elty) in ((:cublasSdot, :Float32),
                      (:cublasDdot, :Float64),
                      (:cublasCdotu, :Complex64),
                      (:cublasCdotc, :Complex64),
                      (:cublasZdotu, :Complex128),
                      (:cublasZdotc, :Complex128))
    @eval begin
        function cuda_dot(num::Int32, x::Ptr{$elty}, y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  ($elty) , (Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  num, x, 1, y, 1)
        end
    end
end

for (fname, elty) in ((:cublasSaxpy, :Float32),
                      (:cublasDaxpy, :Float64),
                      (:cublasCaxpy, :Complex64),
                      (:cublasZaxpy, :Complex128))
    @eval begin
        function cuda_axpy(n::Int32, alpha::($elty), x::Ptr{$elty}, y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, $elty, Ptr{$elty}, Int32,
                         Ptr{$elty}, Int32),
                  n, alpha, x, 1, y, 1)
        end
    end
end

# copy matrix
for (fname, elty) in ((:cublasScopy, :Float32),
                      (:cublasDcopy, :Float64),
                      (:cublasCcopy, :Complex64),
                      (:cublasZcopy, :Complex128))
    @eval begin
        function cuda_copy(n::Int32, x::Ptr{$elty}, y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Ptr{$elty}, Int32,
                         Ptr{$elty}, Int32),
                  n, x, 1, y, 1)
        end
    end
end

# Norm
for (fname, elty) in ((:cublasSnrm2, :Float32),
                      (:cublasDnrm2, :Float64),
                      (:cublasScnrm2, :Complex64),
                      (:cublasDznrm2, :Complex128))
    @eval begin
        function cuda_nrm2(n::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  ($elty), (Int32, Ptr{$elty}, Int32),
                  n, x, 1)
        end
    end
end

# Swap
for (fname, elty) in ((:cublasSswap, :Float32),
                      (:cublasDswap, :Float64),
                      (:cublasCswap, :Complex64),
                      (:cublasZswap, :Complex128))
    @eval begin
        function cuda_swap(n::Int32, x::Ptr{$elty}, y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  n, x, 1, y, 1)
        end
    end
end

# Rotate
for (fname, elty1, elty2, elty3) in ((:cublasSrot, :Float32, :Float32, :Float32),
                                     (:cublasDrot, :Float64, :Float64, :Float64),
                                     (:cublasCrot, :Complex64, :Float32, :Complex64),
                                     (:cublasZrot, :Complex128, :Float64, :Complex128),
                                     (:cublasCsrot, :Complex64, :Float32, :Float32),
                                     (:cublasZdrot, :Complex128, :Float64, :Float64))
    @eval begin
        function cuda_rot(n::Int32, x::Ptr{$elty1}, y::Ptr{$elty1}, sc::($elty2), ss::($elty3))
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Ptr{$elty1}, Int32, Ptr{$elty1}, Int32, $elty2, $elty3),
                  n, x, 1, y, 1, sc, ss)
        end
    end
end

# rotg
for (fname, elty1, elty2) in ((:cublasSrotg, :Float32, :Float32),
                              (:cublasDrotg, :Float64, :Float64),
                              (:cublasCrotg, :Complex64, :Float32),
                              (:cublasZrotg, :Complex128, :Float64))
    @eval begin
        function cuda_rotg(sa::Ptr{$elty1}, sb::Ptr{$elty1}, sc::Ptr{$elty2}, ss::Ptr{$elty1})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Ptr{$elty1}, Ptr{$elty1}, Ptr{$elty2}, Ptr{$elty1}),
                  sa, sb, sc, ss)
        end
    end
end

# rotm
for (fname, elty) in ((:cublasSrotm, :Float32),
                      (:cublasDrotm, :Float64))
    @eval begin
        function cuda_rotm(n::Int32, x::Ptr{$elty}, y::Ptr{$elty}, sparam::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32, Ptr{$elty}),
                  n, x, 1, y, 1, sparam)
        end
    end
end

# rotmg
for (fname, elty) in ((:cublasSrotmg, :Float32),
                      (:cublasDrotmg, :Float64))
    @eval begin
        function cuda_rotmg(sd1::Ptr{$elty}, sd2::Ptr{$elty}, sx1::Ptr{$elty},
                            sy1::Ptr{$elty}, sparam::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}),
                  sd1, sd2, sx1, sy1, sparam)
        end
    end
end

# gemv
for (fname, elty) in ((:cublasSgemv, :Float32),
                      (:cublasDgemv, :Float64),
                      (:cublasCgemv, :Complex64),
                      (:cublasZgemv, :Complex128))
    @eval begin
        function cuda_gemv(trans::Char, m::Int32, n::Int32, alpha::($elty), A::Ptr{$elty},
                           lda::Int32, x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, Int32, $elty, Ptr{$elty}, Int32,
                         Ptr{$elty}, Int32, $elty, Ptr{$elty}, Int32),
                  trans, m, n, alpha, A, lda, x, 1, beta, y, 1)
        end
    end
end

# gbmv
for (fname, elty) in ((:cublasSgbmv, :Float32),
                      (:cublasDgbmv, :Float64),
                      (:cublasCgbmv, :Complex64),
                      (:cublasZgbmv, :Complex128))
    @eval begin
        function cuda_gbmv(trans::Char, m::Int32, n::Int32, kl::Int32, ku::Int32, alpha::($elty),
                           A::Ptr{$elty}, lda::Int32, x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, Int32, Int32, Int32, $elty, Ptr{$elty}, Int32,
                         Ptr{$elty}, Int32, $elty, Ptr{$elty}, Int32),
                  trans, m, n, kl, ku, alpha, A, lda, x, 1, beta, y, 1)
        end
    end
end

# trmv
for (fname, elty) in ((:cublasStrmv, :Float32),
                      (:cublasDtrmv, :Float64),
                      (:cublasCtrmv, :Complex64),
                      (:cublasZtrmv, :Complex128))
    @eval begin
        function cuda_trmv(uplo::Char, trans::Char, diag::Char, n::Int32,
                           A::Ptr{$elty}, lda::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Char, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, trans, diag, n, A, lda, x, 1)
        end
    end
end

# tbmv
for (fname, elty) in ((:cublasStbmv, :Float32),
                      (:cublasDtbmv, :Float64),
                      (:cublasCtbmv, :Complex64),
                      (:cublasZtbmv, :Complex128))
    @eval begin
        function cuda_tbmv(uplo::Char, trans::Char, diag::Char, n::Int32, k::Int32,
                           A::Ptr{$elty}, lda::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Char, Int32, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, trans, diag, n, k, A, lda, x, 1)
        end
    end
end

# tpmv
for (fname, elty) in ((:cublasStpmv, :Float32),
                      (:cublasDtpmv, :Float64),
                      (:cublasCtpmv, :Complex64),
                      (:cublasZtpmv, :Complex128))
    @eval begin
        function cuda_tpmv(uplo::Char, trans::Char, diag::Char, n::Int32,
                           A::Ptr{$elty}, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Char, Int32, Ptr{$elty}, Ptr{$elty}, Int32),
                  uplo, trans, diag, n, A, x, 1)
        end
    end
end

# trsv
for (fname, elty) in ((:cublasStrsv, :Float32),
                      (:cublasDtrsv, :Float64),
                      (:cublasCtrsv, :Complex64),
                      (:cublasZtrsv, :Complex128))
    @eval begin
        function cuda_trsv(uplo::Char, trans::Char, diag::Char, n::Int32,
                           A::Ptr{$elty}, lda::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Char, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, trans, diag, n, A, lda, x, 1)
        end
    end
end

# tpsv
for (fname, elty) in ((:cublasStpsv, :Float32),
                      (:cublasDtpsv, :Float64),
                      (:cublasCtpsv, :Complex64),
                      (:cublasZtpsv, :Complex128))
    @eval begin
        function cuda_tpsv(uplo::Char, trans::Char, diag::Char, n::Int32,
                           A::Ptr{$elty}, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Char, Int32, Ptr{$elty}, Ptr{$elty}, Int32),
                  uplo, trans, diag, n, A, x, 1)
        end
    end
end

# tbsv
for (fname, elty) in ((:cublasStbsv, :Float32),
                      (:cublasDtbsv, :Float64),
                      (:cublasCtbsv, :Complex64),
                      (:cublasZtbsv, :Complex128))
    @eval begin
        function cuda_tbsv(uplo::Char, trans::Char, diag::Char, n::Int32, k::Int32,
                           A::Ptr{$elty}, lda::Int32, x::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Char, Char, Int32, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, trans, diag, n, k, A, lda, x, 1)
        end
    end
end

# symv
for (fname, elty) in ((:cublasSsymv, :Float32),
                      (:cublasDsymv, :Float64))
    @eval begin
        function cuda_symv(uplo::Char, n::Int32, alpha::($elty), A::Ptr{$elty},
                           lda::Int32, x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, $elty, Ptr{$elty}, Int32),
                  uplo, n, alpha, A, lda, x, 1, beta, y, 1)
        end
    end
end

# hemv
for (fname, elty) in ((:cublasChemv, :Complex64),
                      (:cublasZhemv, :Complex128))
    @eval begin
        function cuda_hemv(uplo::Char, n::Int32, alpha::($elty), A::Ptr{$elty},
                           lda::Int32, x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, $elty, Ptr{$elty}, Int32),
                  uplo, n, alpha, A, lda, x, 1, beta, y, 1)
        end
    end
end

# sbmv
for (fname, elty) in ((:cublasSsbmv, :Float32),
                      (:cublasDsbmv, :Float64))
    @eval begin
        function cuda_sbmv(uplo::Char, n::Int32, k::Int32, alpha::($elty), A::Ptr{$elty},
                           lda::Int32, x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, $elty, Ptr{$elty}, Int32),
                  uplo, n, k, alpha, A, lda, x, 1, beta, y, 1)
        end
    end
end

# hbmv
for (fname, elty) in ((:cublasChemv, :Complex64),
                      (:cublasZhemv, :Complex128))
    @eval begin
        function cuda_hbmv(uplo::Char, n::Int32, k::Int32, alpha::($elty), A::Ptr{$elty},
                           lda::Int32, x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, $elty, Ptr{$elty}, Int32),
                  uplo, n, k, alpha, A, lda, x, 1, beta, y, 1)
        end
    end
end

# spmv
for (fname, elty) in ((:cublasSspmv, :Float32),
                      (:cublasDspmv, :Float64))
    @eval begin
        function cuda_spmv(uplo::Char, n::Int32, alpha::($elty), A::Ptr{$elty},
                           x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Ptr{$elty},
                         Int32, $elty, Ptr{$elty}, Int32),
                  uplo, n, alpha, A, x, 1, beta, y, 1)
        end
    end
end

# hpmv
for (fname, elty) in ((:cublasChpmv, :Complex64),
                      (:cublasZhpmv, :Complex128))
    @eval begin
        function cuda_hpmv(uplo::Char, n::Int32, alpha::($elty), A::Ptr{$elty},
                           x::Ptr{$elty}, beta::($elty), y::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Ptr{$elty},
                         Int32, $elty, Ptr{$elty}, Int32),
                  uplo, n, alpha, A, x, 1, beta, y, 1)
        end
    end
end

# ger
for (fname, elty) in ((:cublasSger, :Float32),
                      (:cublasDger, :Float64))
    @eval begin
        function cuda_ger(m::Int32, n::Int32, alpha::($elty), x::Ptr{$elty},
                          y::Ptr{$elty}, A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, Ptr{$elty}, Int32),
                  m, n, alpha, x, 1, y, 1,  A, lda)
        end
    end
end

# geru
for (fname, elty) in ((:cublasSgeru, :Complex64),
                      (:cublasDgeru, :Complex128))
    @eval begin
        function cuda_geru(m::Int32, n::Int32, alpha::($elty), x::Ptr{$elty},
                          y::Ptr{$elty}, A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, Ptr{$elty}, Int32),
                  m, n, alpha, x, 1, y, 1,  A, lda)
        end
    end
end

# gerc
for (fname, elty) in ((:cublasSgerc, :Complex64),
                      (:cublasDgerc, :Complex128))
    @eval begin
        function cuda_gerc(m::Int32, n::Int32, alpha::($elty), x::Ptr{$elty},
                          y::Ptr{$elty}, A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Int32, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty},
                         Int32, Ptr{$elty}, Int32),
                  m, n, alpha, x, 1, y, 1,  A, lda)
        end
    end
end

# syr
for (fname, elty) in ((:cublasSsyr, :Float32),
                      (:cublasDsyr, :Float64))
    @eval begin
        function cuda_syr(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty},
                          A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, A, lda)
        end
    end
end

# her
for (fname, elty) in ((:cublasSher, :Complex64),
                      (:cublasDher, :Complex128))
    @eval begin
        function cuda_her(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty},
                          A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, A, lda)
        end
    end
end

# spr
for (fname, elty) in ((:cublasSspr, :Float32),
                      (:cublasDspr, :Float64))
    @eval begin
        function cuda_syr(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty},
                          A::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, A)
        end
    end
end

# hpr
for (fname, elty) in ((:cublasShpr, :Complex64),
                      (:cublasDhpr, :Complex128))
    @eval begin
        function cuda_hpr(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty},
                          A::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, A)
        end
    end
end

# syr2
for (fname, elty) in ((:cublasSsyr2, :Float32),
                      (:cublasDsyr2, :Float64))
    @eval begin
        function cuda_syr2(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty}, y::Ptr{$elty},
                           A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, y, 1, A, lda)
        end
    end
end

# her2
for (fname, elty) in ((:cublasSher2, :Complex64),
                      (:cublasDher2, :Complex128))
    @eval begin
        function cuda_her2(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty}, y::Ptr{$elty},
                           A::Ptr{$elty}, lda::Int32)
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, y, 1, A, lda)
        end
    end
end

# spr2
for (fname, elty) in ((:cublasSspr2, :Float32),
                      (:cublasDspr2, :Float64))
    @eval begin
        function cuda_spr2(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty},
                           y::Ptr{$elty}, A::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, y::Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, y, 1, A)
        end
    end
end

# hpr2
for (fname, elty) in ((:cublasShpr2, :Complex64),
                      (:cublasDhpr2, :Complex128))
    @eval begin
        function cuda_hpr2(uplo::Char, n::Int32, alpha::($elty), x::Ptr{$elty},
                           y::Ptr{$elty}, A::Ptr{$elty})
            ccall(dlsym(libcublas, $string(fname)),
                  Void, (Char, Int32, $elty, Ptr{$elty}, Int32, y::Ptr{$elty}, Int32, Ptr{$elty}, Int32),
                  uplo, n, alpha, x, 1, y, 1, A)
        end
    end
end

max_stars_repo_name: stefan-k/julia-CuMatrix

max_issues_repo_name: stefan-k/julia-CuMatrix

max_forks_repo_name: stefan-k/julia-CuMatrix


text:
# Plotting functions
function smoothLivePlotGeneral(plotFunction, plotArgs)
    data_obs = Observable{Any}(plotArgs)
    plt_obs = Observable{Any}(plotFunction(plotArgs...))
    map!(mapArg -> plotFunction(mapArg...), plt_obs, data_obs)

    # Create figure
    ui = dom"div"(plt_obs);
    display(ui)
    sleep(0.4)
    return data_obs
end
function modifyPlotObject!(mutablePlotArray; args...)
    # Modify plot objects by passing "argX = newValue" pairs
    Nargs = length(mutablePlotArray[])
    NargsSupplied = length(args)

    # Test arguments
    if .!all("arg" .== map(x -> string(x)[1:3], args.itr))
        error("""Mutable arguments must begin with \"arg\"""")
    end
    argStringNums = map(x -> parse(Int, string(x)[4:end]), args.itr)
    if .!all(argStringNums .<= Nargs)
        error("Indices must be less than or equal to length of mutable array")
    end

    mutablePlotArray[] = map(x -> x in argStringNums ? args.data[findfirst(x .== argStringNums)] : mutablePlotArray[][x], 1:Nargs)

end

ext: jl


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractDynamicModelImpl

export AbstractDynamicModel

abstract AbstractDynamicModel

end

module AbstractDynamicModelInterfaces

export initializeDynamicModel,
       simulateDynamicModel

function initializeDynamicModel()
end

function simulateDynamicModel()
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
# Example by @walra356 (modified)
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(13)
    n= 5
    x = rand(n)                        # specify x steplength
    y = rand(n)                        # specify y steplength
    σ = rand(length(x), length(y))      # beware of dims
    # custom ticks
    function steps(x::Vector{T} where T<:Real)
    sum(x .< 0) == 0 || error("Error: $x - negative step length not allowed")
    (s = append!(eltype(x)[0],x); [Base.sum(s[1:i]) for i ∈ Base.eachindex(s)])
    end
    function stepcenters(x::Vector{T} where T<:Real)
    δ = x .* 0.5
    s = append!(eltype(x)[0],x)
    [Base.sum(s[1:i]) for i ∈ Base.eachindex(x)] .+ δ
    end
    stepedges(x::Vector{T} where T<:Real) = steps(x)

    # the actual plot
    theme = Theme(fontsize = 16, colormap = :gist_earth, resolution = (750,450))
    set_theme!(theme)

    attr1 = (xticks = (stepcenters(x),string.(1:n)),
        yticks = (stepcenters(y), string.(1:n)), xlabel = "cat", )
    attr2 = (xticks = (stepedges(x),string.(0:n)),
        yticks = (stepedges(y), string.(0:n)), xlabel="cat", )

    fig = Figure()
    ax1 = Axis(fig; attr1...)
    ax2 = Axis(fig; attr2...)
    heatmap!(ax1, steps(x), steps(y), σ)
    heatmap!(ax2, steps(x), steps(y), σ)
    fig[1,1] = ax1
    fig[1,2] = ax2
    fig
    save(joinpath(@__DIR__, "output", "heatmapIrregularCategories.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using ManifoldLearning
import ManifoldLearning: BruteForce, knn, swiss_roll
using Test
import Random

Random.seed!(3483743871)

@testset "Nearest Neighbors" begin

    # setup parameters
    k = 12
    X, L = swiss_roll(100)
    DD, EE = knn(X,k)
    @test_throws AssertionError knn(rand(3,10), k)

    NN = fit(BruteForce, X, k)
    D, E = knn(NN, X)
    @test size(X,2) == size(D,2) && size(D, 1) == k
    @test size(X,2) == size(E,2) && size(E, 1) == k
    @test E == EE
    @test D ≈ DD

    D, E = knn(NN, X[:,1:k+1])
    @test k+1 == size(D,2) && size(D, 1) == k
    @test k+1 == size(E,2) && size(E, 1) == k

    D, E = knn(NN, X[:,1:k+1], self=true)
    @test D[1,:] == zeros(k+1)
    @test E[1,:] == collect(1:k+1)

    @test_throws AssertionError knn(NN, X[:,1:k])
end

@testset "Manifold Learning" begin

    # setup parameters
    k = 12
    d = 2
    X, L = swiss_roll()

    # test algorithms
    #@testset for algorithm in [Isomap, LEM, LLE, HLLE, LTSA, DiffMap]
    @testset for algorithm in [DiffMap]
        for (k, T) in zip([5, 12], [Float32, Float64])
            # construct KW parameters
            kwargs = [:maxoutdim=>d]
            if algorithm == DiffMap
                push!(kwargs, :t => k)
            else
                push!(kwargs, :k => k)
            end

            # call transformation
            Y = fit(algorithm, convert(Array{T}, X); kwargs...)

            # test results
            @test outdim(Y) == d
            @test size(transform(Y), 2) == size(X, 2)
            @test length(split(sprint(show, Y), '\n')) > 1
            @test length(eigvals(Y)) == d
            if algorithm !== DiffMap
                @test neighbors(Y) == k
                @test length(vertices(Y)) > 1
            end

            # test if we provide pre-computed Gram matrix
            if algorithm == DiffMap
                kernel = (x, y) -> exp(-sum((x .- y) .^ 2)) # default kernel
                n_obs = size(X)[2]
                custom_K = zeros(T, n_obs, n_obs)
                for i = 1:n_obs
                    for j = i:n_obs
                        custom_K[i, j] = kernel(convert(Array{T}, X[:, i]), convert(Array{T}, X[:, j]))
                        custom_K[j, i] = custom_K[i, j]
                    end
                end
                Y_custom_K = fit(algorithm, custom_K; kwargs..., kernel=nothing)
                @test Y_custom_K.proj ≈ Y.proj
            end
        end
    end
end

@testset "OOS" begin
    k = 10
    d = 2
    X, _ = swiss_roll()
    M = fit(Isomap, X; k=k, maxoutdim=d)

    @test all(sum(abs2, transform(M) .- transform(M,X), dims=1) .< eps())
end


text:
using Base: @deprecate, @deprecate_binding, depwarn

if !isdefined(Base, :stderr)
    @deprecate stderr(obj::StatisticalModel) stderror(obj)
else
    function (io::typeof(stderr))(obj::StatisticalModel)
        Base.depwarn("stderr(obj::StatisticalModel) is deprecated, use stderror(obj) instead", :stderr)
        io === stderr ? stderror(obj) : throw(MethodErrror(io, (obj,)))
    end
end

@deprecate model_response(obj::StatisticalModel) response(obj)

@deprecate norepeats(a::AbstractArray) allunique(a)

@deprecate(mad!(v::AbstractArray{<:Real}, center;
                constant::Real = BigFloat("1.482602218505601860547076529360423431326703202590312896536266275245674447622701")),
           mad!(v, center=center, constant=constant))


text:
function block(i)
    color = [:dark_gray, :red, :light_red, :yellow, :green, :cyan, :blue, :magenta][i+1]
    return string(Crayon(foreground = color), i==0 ? " □ " : " ■ ")
end

function update_board!(old::Board, b::Board)
    buf = IOBuffer()
    for i in findall(old.data .⊻ b.data .!= 0)
        y,x = Tuple(i)
        put(buf, [(3*x)-2,y], string(block.(b.data[y,x])...))
    end
    if (old.level != b.level) || (old.score != b.score)
        put(buf, [3,21], :dark_gray, "Level: $(b.level)\tScore:$(b.score)")
    end
    print(String(take!(buf)))
end

update_board!(b::Board) = update_board!(Board(2), b)

function print_preview(b::Board)
    buf = IOBuffer()
    put(buf, [33, 2], :dark_gray, string("Next Tiles:"))
    for (nr, tile) in enumerate(b.nexttiles)
        _print_tile(buf, tile, 2+(nr-1)*4)
    end
    print(String(take!(buf)))
end

function print_hold(b::Board)
    buf = IOBuffer()
    put(buf, [33, 16], :dark_gray, string("Hold Tile:"))
    _print_tile(buf, b.holdtile, 16)
    print(String(take!(buf)))
end

function _print_tile(buf, tile, x)
    dt = data(tile)'
    _, dy = size(dt)
    for i in 1:2
        put(buf, [35, x+i+1], :dark_gray, string(" □ "^4))
        i <= dy && put(buf, [35, x+i+1], string(block.(dt[:, i])...))
    end
end

max_stars_repo_name: ChristianKurz/REPLTetris

max_issues_repo_name: ChristianKurz/REPLTetris

max_forks_repo_name: ChristianKurz/REPLTetris


text:
################################################################################
### AVX-powered matrix multiplications for real matrices
################################################################################

function vmul!(C::Matrix{T}, A::Matrix{T}, B::Matrix{T}) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
end
function vmul!(C::Matrix{T}, A::Matrix{T}, B::Diagonal{T}) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(A, 2)
        C[m,n] = A[m,n] * B.diag[n]
    end
end
function vmul!(C::Matrix{T}, A::Diagonal{T}, B::Matrix{T}) where {T <: Real}
    @turbo for m in 1:size(C, 1), n in 1:size(C, 2)
        C[m,n] = A.diag[m] * B[m,n]
    end
end
function vmul!(C::Matrix{T}, A::Matrix{T}, X::Adjoint{T}) where {T <: Real}
    B = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * conj(B[n, k])
        end
        C[m,n] = Cmn
    end
end
function vmul!(C::Matrix{T}, X::Adjoint{T}, B::Matrix{T}) where {T <: Real}
    A = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += conj(A[k,m]) * B[k,n]
        end
        C[m,n] = Cmn
    end
end
function vmul!(C::Matrix{T}, X1::Adjoint{T}, X2::Adjoint{T}) where {T <: Real}
    A = X1.parent
    B = X2.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[k,m] * B[n,k]
        end
        C[m,n] = Cmn
    end
end
function rvmul!(A::Matrix{T}, B::Diagonal{T}) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(A, 2)
        A[m,n] = A[m,n] * B.diag[n]
    end
end
function lvmul!(A::Diagonal{T}, B::Matrix{T}) where {T <: Real}
    @turbo for m in 1:size(B, 1), n in 1:size(B, 2)
        B[m,n] = A.diag[m] * B[m, n]
    end
end
function rvadd!(A::Matrix{T}, D::Diagonal{T}) where {T <: Real}
    @turbo for i in axes(A, 1)
        A[i, i] = A[i, i] + D.diag[i]
    end
end
function rvadd!(A::Matrix{T}, B::Matrix{T}) where {T <: Real}
    @turbo for i in axes(A, 1), j in axes(A, 2)
        A[i, j] = A[i, j] + B[i, j]
    end
end
function vsub!(O::Matrix{T}, A::Matrix{T}, ::UniformScaling) where {T <: Real}
    T1 = one(T)
    @turbo for i in axes(O, 1), j in axes(O, 2)
        O[i, j] = A[i, j]
    end
    @turbo for i in axes(O, 1)
        O[i, i] -= T1
    end
end

# NOTE
# These should only be called with real Vectors, no need to implement extra
# methods
function vmin!(v::Vector{T}, w::Vector{T}) where {T<:Real}
    T1 = one(T)
    @turbo for i in eachindex(w)
        v[i] = min(T1, w[i])
    end
    v
end
function vmininv!(v::Vector{T}, w::Vector{T}) where {T<:Real}
    T1 = one(T)
    @turbo for i in eachindex(w)
        v[i] = T1 / min(T1, w[i])
    end
    v
end
function vmax!(v::Vector{T}, w::Vector{T}) where {T<:Real}
    T1 = one(T)
    @turbo for i in eachindex(w)
        v[i] = max(T1, w[i])
    end
    v
end
function vmaxinv!(v::Vector{T}, w::Vector{T}) where {T<:Real}
    T1 = one(T)
    @turbo for i in eachindex(w)
        v[i] = T1 / max(T1, w[i])
    end
    v
end
function vinv!(v::Vector{T}) where {T<:Real}
    T1 = one(T)
    @turbo for i in eachindex(v)
        v[i] = T1 / v[i]
    end
    v
end

"""
    rdivp!(A, T, O, pivot)

Computes `A * T^-1` where `T` is an upper triangular matrix which should be
pivoted according to a pivoting vector `pivot`. `A`, `T` and `O` should all be
square matrices of the same size. The result will be written to `A` without
changing `T`.

This function is written to work with (@ref)[`udt_AVX_pivot!`].
"""
function rdivp!(A::Matrix, T, O, pivot)
    # assume Diagonal is ±1!
    @inbounds begin
        N = size(A, 1)

        # Apply pivot
        for (j, p) in enumerate(pivot)
            @turbo for i in 1:N
                O[i, j] = A[i, p]
            end
        end

        # do the rdiv
        # @turbo will segfault on `k in 1:0`, so pull out first loop
        @turbo for i in 1:N
            A[i, 1] = O[i, 1] / T[1, 1]
        end
        for j in 2:N
            @turbo for i in 1:N
                x = O[i, j]
                for k in 1:j-1
                    x -= A[i, k] * T[k, j]
                end
                A[i, j] = x / T[j, j]
            end
        end
    end
    A
end

################################################################################
### Fallbacks for complex matrices
################################################################################

vmul!(C, A, B) = mul!(C, A, B)
rvmul!(A, B) = rmul!(A, B)
lvmul!(A, B) = lmul!(A, B)
rvadd!(A, B) = A .+= B
vsub!(A, B, C) = A .= B .- C
function vsub!(A, B, ::UniformScaling)
    T1 = one(eltype(A))
    T0 = zero(eltype(A))
    @inbounds for i in axes(A, 1), j in axes(A, 2)
        A[i, j] = B[i, j] - ifelse(i==j, T1, T0)
    end
    A
end

function rdivp!(A::Matrix{<: Complex}, T::Matrix{<: Complex}, O::Matrix{<: Complex}, pivot)
    # assume Diagonal is ±1!
    @inbounds begin
        N = size(A, 1)

        # Apply pivot
        for (j, p) in enumerate(pivot)
            for i in 1:N
                O[i, j] = A[i, p]
            end
        end

        # do the rdiv
        # @turbo will segfault on `k in 1:0`, so pull out first loop
        for j in 1:N
            for i in 1:N
                x = O[i, j]
                @simd for k in 1:j-1
                    x -= A[i, k] * T[k, j]
                end
                A[i, j] = x / T[j, j]
            end
        end
    end
    A
end


text:
@enum Shape CIRCLE RECTANGLE ROUNDED_RECTANGLE DISTANCEFIELD TRIANGLE
@enum CubeSides TOP BOTTOM FRONT BACK RIGHT LEFT

struct Grid{N, T <: AbstractRange}
    dims::NTuple{N, T}
end
Base.ndims(::Grid{N,T}) where {N,T} = N

Grid(ranges::AbstractRange...) = Grid(ranges)
function Grid(a::Array{T, N}) where {N, T}
    s = Vec{N, Float32}(size(a))
    smax = maximum(s)
    s = s./smax
    Grid(ntuple(Val{N}) do i
        range(0, stop=s[i], length=size(a, i))
    end)
end

Grid(a::AbstractArray, ranges...) = Grid(a, ranges)

"""
This constructor constructs a grid from ranges given as a tuple.
Due to the approach, the tuple `ranges` can consist of NTuple(2, T)
and all kind of range types. The constructor will make sure that all ranges match
the size of the dimension of the array `a`.
"""
function Grid(a::AbstractArray{T, N}, ranges::Tuple) where {T, N}
    length(ranges) =! N && throw(ArgumentError(
        "You need to supply a range for every dimension of the array. Given: $ranges
        given Array: $(typeof(a))"
    ))
    Grid(ntuple(Val(N)) do i
        range(first(ranges[i]), stop=last(ranges[i]), length=size(a, i))
    end)
end

Base.length(p::Grid) = prod(size(p))
Base.size(p::Grid) = map(length, p.dims)
function Base.getindex(p::Grid{N,T}, i) where {N,T}
    inds = ind2sub(size(p), i)
    return Point{N, eltype(T)}(ntuple(Val(N)) do i
        p.dims[i][inds[i]]
    end)
end

Base.iterate(g::Grid, i = 1) = i <= length(g) ? (g[i], i + 1) : nothing

GLAbstraction.isa_gl_struct(x::Grid) = true
GLAbstraction.toglsltype_string(t::Grid{N,T}) where {N,T} = "uniform Grid$(N)D"
function GLAbstraction.gl_convert_struct(g::Grid{N, T}, uniform_name::Symbol) where {N,T}
    return Dict{Symbol, Any}(
        Symbol("$uniform_name.start") => Vec{N, Float32}(minimum.(g.dims)),
        Symbol("$uniform_name.stop") => Vec{N, Float32}(maximum.(g.dims)),
        Symbol("$uniform_name.lendiv") => Vec{N, Cint}(length.(g.dims) .- 1),
        Symbol("$uniform_name.dims") => Vec{N, Cint}(map(length, g.dims))
    )
end
function GLAbstraction.gl_convert_struct(g::Grid{1, T}, uniform_name::Symbol) where T
    x = g.dims[1]
    return Dict{Symbol, Any}(
        Symbol("$uniform_name.start") => Float32(minimum(x)),
        Symbol("$uniform_name.stop") => Float32(maximum(x)),
        Symbol("$uniform_name.lendiv") => Cint(length(x) - 1),
        Symbol("$uniform_name.dims") => Cint(length(x))
    )
end
import Base: getindex, length, iterate, ndims, setindex!, eltype

struct GridZRepeat{G, T, N} <: AbstractArray{Point{3, T}, N}
    grid::G
    z::Array{T, N}
end
Base.size(g::GridZRepeat) = size(g.z)
Base.size(g::GridZRepeat, i) = size(g.z, i)
Base.IndexStyle(::Type{<:GridZRepeat}) = Base.IndexLinear()

function Base.getindex(g::GridZRepeat{G, T}, i) where {G,T}
    pxy = g.grid[i]
    Point{3, T}(pxy[1], pxy[2], g.z[i])
end

struct GLVisualizeShader <: AbstractLazyShader
    paths::Tuple
    kw_args::Dict{Symbol, Any}
    function GLVisualizeShader(paths::String...; view = Dict{String, String}(), kw_args...)
        # TODO properly check what extensions are available
        @static if !Sys.isapple()
            view["GLSL_EXTENSIONS"] = "#extension GL_ARB_conservative_depth: enable"
            view["SUPPORTED_EXTENSIONS"] = "#define DETPH_LAYOUT"
        end
        args = Dict{Symbol, Any}(kw_args)
        args[:view] = view
        args[:fragdatalocation] = [(0, "fragment_color"), (1, "fragment_groupid")]
        new(map(x-> assetpath("shader", x), paths), args)
    end
end


text:
using Oceananigans.Operators: assumed_field_location

"""
    FieldBoundaryConditions

An alias for `NamedTuple{(:x, :y, :z)}` that represents a set of three `CoordinateBoundaryCondition`s
applied to a field along x, y, and z.
"""
const FieldBoundaryConditions = NamedTuple{(:x, :y, :z)}

"""
    FieldBoundaryConditions(x, y, z)

Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for each of the
`x`, `y`, and `z` coordinates.
"""
FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))

"""
    DefaultBoundaryCondition(::Type{Periodic}, loc)

Returns [`PeriodicBoundaryCondition`](@ref).
"""
DefaultBoundaryCondition(::Type{Grids.Periodic}, loc) = PeriodicBoundaryCondition()

"""
    DefaultBoundaryCondition(::Type{Flat}, loc)

Returns `nothing`.
"""
DefaultBoundaryCondition(::Type{Flat}, loc) = nothing

"""
    DefaultBoundaryCondition(::Type{Bounded}, ::Type{Cell})

Returns [`NoFluxBoundaryCondition`](@ref).
"""
DefaultBoundaryCondition(::Type{Bounded}, ::Type{Cell}) = NoFluxBoundaryCondition()

"""
    DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face})

Returns [`ImpenetrableBoundaryCondition`](@ref).
"""
DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face}) = ImpenetrableBoundaryCondition()

function validate_bcs(topology, left_bc, right_bc, default_bc, left_name, right_name, dir)

    if topology isa Periodic && (left_bc != default_bc || right_bc != default_bc)
        e = "Cannot specify $left_name or $right_name boundary conditions with " *
            "$topology topology in $dir-direction."
        throw(ArgumentError(e))
    end

    return true
end

"""
    FieldBoundaryConditions(grid, loc;   east = DefaultBoundaryCondition(topology(grid, 1), loc[1]),
                                         west = DefaultBoundaryCondition(topology(grid, 1], loc[1]),
                                        south = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                        north = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                       bottom = DefaultBoundaryCondition(topology(grid, 3), loc[3]),
                                          top = DefaultBoundaryCondition(topology(grid, 3), loc[3]))

Construct `FieldBoundaryConditions` for a field with location `loc` (a 3-tuple of `Face` or `Cell`)
defined on `grid`.

Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)` and `loc`.
"""
function FieldBoundaryConditions(grid, loc;   east = DefaultBoundaryCondition(topology(grid, 1), loc[1]),
                                              west = DefaultBoundaryCondition(topology(grid, 1), loc[1]),
                                             south = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                             north = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                            bottom = DefaultBoundaryCondition(topology(grid, 3), loc[3]),
                                               top = DefaultBoundaryCondition(topology(grid, 3), loc[3]))

    TX, TY, TZ = topology(grid)
    x_default_bc = DefaultBoundaryCondition(topology(grid, 1), loc[1])
    y_default_bc = DefaultBoundaryCondition(topology(grid, 2), loc[2])
    z_default_bc = DefaultBoundaryCondition(topology(grid, 3), loc[3])

    validate_bcs(TX, west,   east, x_default_bc, :west,   :east, :x)
    validate_bcs(TY, south, north, y_default_bc, :south, :north, :y)
    validate_bcs(TZ, bottom,  top, z_default_bc, :bottom,  :top, :z)

    x = CoordinateBoundaryConditions(west, east)
    y = CoordinateBoundaryConditions(south, north)
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

"""
    UVelocityBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Face),
                                        west = DefaultBoundaryCondition(topology(grid, 1), Face),
                                       south = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                       north = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                      bottom = DefaultBoundaryCondition(topology(grid, 3), Cell),
                                         top = DefaultBoundaryCondition(topology(grid, 3), Cell))

Construct `FieldBoundaryConditions` for the `u`-velocity field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face, Cell, Cell); user_defined_bcs...)

"""
    VVelocityBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                        west = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                       south = DefaultBoundaryCondition(topology(grid, 2), Face),
                                       north = DefaultBoundaryCondition(topology(grid, 2), Face),
                                      bottom = DefaultBoundaryCondition(topology(grid, 3), Cell),
                                         top = DefaultBoundaryCondition(topology(grid, 3), Cell))

Construct `FieldBoundaryConditions` for the `v`-velocity field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Face, Cell); user_defined_bcs...)

"""
    WVelocityBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                        west = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                       south = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                       north = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                      bottom = DefaultBoundaryCondition(topology(grid, 3), Face),
                                         top = DefaultBoundaryCondition(topology(grid, 3), Face))

Construct `FieldBoundaryConditions` for the `w`-velocity field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Face); user_defined_bcs...)

"""
    TracerBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                     west = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                    south = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                    north = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                   bottom = DefaultBoundaryCondition(topology(grid, 3), Cell),
                                      top = DefaultBoundaryCondition(topology(grid, 3), Cell))

Construct `FieldBoundaryConditions` for a tracer field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)

const PressureBoundaryConditions = TracerBoundaryConditions
const DiffusivityBoundaryConditions = TracerBoundaryConditions

# Here we overload setproperty! and getproperty to permit users to call
# the 'left' and 'right' bcs in the z-direction 'bottom' and 'top'
# and the 'left' and 'right' bcs in the y-direction 'south' and 'north'.
Base.setproperty!(bcs::FieldBoundaryConditions, side::Symbol, bc) = setbc!(bcs, Val(side), bc)

setbc!(bcs::FieldBoundaryConditions, ::Val{S}, bc) where S = setfield!(bcs, S, bc)
setbc!(bcs::FieldBoundaryConditions, ::Val{:west},   bc) = setfield!(bcs.x, :left,  bc)
setbc!(bcs::FieldBoundaryConditions, ::Val{:east},   bc) = setfield!(bcs.x, :right, bc)
setbc!(bcs::FieldBoundaryConditions, ::Val{:south},  bc) = setfield!(bcs.y, :left,  bc)
setbc!(bcs::FieldBoundaryConditions, ::Val{:north},  bc) = setfield!(bcs.y, :right, bc)
setbc!(bcs::FieldBoundaryConditions, ::Val{:bottom}, bc) = setfield!(bcs.z, :left,  bc)
setbc!(bcs::FieldBoundaryConditions, ::Val{:top},    bc) = setfield!(bcs.z, :right, bc)

@inline Base.getproperty(bcs::FieldBoundaryConditions, side::Symbol) = getbc(bcs, Val(side))

@inline getbc(bcs::FieldBoundaryConditions, ::Val{S}) where S = getfield(bcs, S)
@inline getbc(bcs::FieldBoundaryConditions, ::Val{:west})   = getfield(bcs.x, :left)
@inline getbc(bcs::FieldBoundaryConditions, ::Val{:east})   = getfield(bcs.x, :right)
@inline getbc(bcs::FieldBoundaryConditions, ::Val{:south})  = getfield(bcs.y, :left)
@inline getbc(bcs::FieldBoundaryConditions, ::Val{:north})  = getfield(bcs.y, :right)
@inline getbc(bcs::FieldBoundaryConditions, ::Val{:bottom}) = getfield(bcs.z, :left)
@inline getbc(bcs::FieldBoundaryConditions, ::Val{:top})    = getfield(bcs.z, :right)

function regularize_field_boundary_conditions(bcs::FieldBoundaryConditions, grid, tracer_names, field_name)

    model_field_names = tuple(:u, :v, :w, tracer_names...)

    X, Y, Z = assumed_field_location(field_name)

    east   = regularize_boundary_condition(bcs.east,   Nothing, Y, Z, grid.Nx, model_field_names)
    west   = regularize_boundary_condition(bcs.west,   Nothing, Y, Z, 1,       model_field_names)
    north  = regularize_boundary_condition(bcs.north,  X, Nothing, Z, grid.Ny, model_field_names)
    south  = regularize_boundary_condition(bcs.south,  X, Nothing, Z, 1,       model_field_names)
    top    = regularize_boundary_condition(bcs.top,    X, Y, Nothing, grid.Nz, model_field_names)
    bottom = regularize_boundary_condition(bcs.bottom, X, Y, Nothing, 1,       model_field_names)

    x = CoordinateBoundaryConditions(west, east)
    y = CoordinateBoundaryConditions(south, north)
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

function regularize_field_boundary_conditions(boundary_conditions::NamedTuple, grid, tracer_names, field_name=nothing)
    boundary_conditions_names = propertynames(boundary_conditions)
    boundary_conditions_tuple = Tuple(regularize_field_boundary_conditions(bcs, grid, tracer_names, name)
                                      for (name, bcs) in zip(boundary_conditions_names, boundary_conditions))
    boundary_conditions = NamedTuple{boundary_conditions_names}(boundary_conditions_tuple)
    return boundary_conditions
end

regularize_field_boundary_conditions(::Missing, grid, tracer_names, field_name=nothing) = missing


text:
using Test
using QuantumOptics

@testset "stochastic_definitions" begin

n=20
b=FockBasis(n)
psi0 = fockstate(b, 0)
a = destroy(b)
ad = dagger(a)
H0 = ad*a

fdeterm, fstoch = stochastic.homodyne_carmichael(H0, a, 0.5π)
Y = 1.0im*(ad - a)
@test fdeterm(0.0, psi0) == H0 + expect(Y, psi0)*a - 0.5im*ad*a
@test fstoch(0.0, psi0)[1].data ≈ a.data

psi1 = 0.5*psi0
fdeterm1, fstoch = stochastic.homodyne_carmichael(H0, a, 0.5π; normalize_expect=false)
@test fdeterm1(0.0, psi1) == H0 + 0.5*expect(Y, psi0)*a - 0.5im*ad*a

end # testset


text:
# We register our dependence on tika jar as a data dependency
function register_tika_app_jar_dependency()
    register(DataDep(
        "tika-app-1.22.jar",
        "Apache Tika App Jar file version 1.22",
        ["http://it.apache.contactlab.it/tika/tika-app-1.22.jar"],
        "d7219709abc547136fa5fca17632a85fe1cd36dc08cb4031957e3c9a836543e2"))
end

withtxtext(fn) = first(splitext(fn)) * ".txt"

function pdf2text_apache_tika(filename::String,
    txtfilename::String = withtxtext(filename))

    tikajardir = datadep"tika-app-1.22.jar"
    tikajarpath = joinpath(tikajardir, "tika-app-1.22.jar")

    println("Converting $filename to $txtfilename using Apache Tika v1.22")
    try
        run(pipeline(`/usr/bin/env java -jar $tikajarpath -t $filename`, txtfilename))
    catch err
        error("Error when running Apache Tika to convert to txt file")
        @show err
    end

    if !isfile(txtfilename)
        error("No text file created by Apache Tika!")
    end

    return txtfilename
end

pdf2text(filename::String, txtfilename::String = withtxtext(filename)) =
    pdf2text_apache_tika(filename, txtfilename)

function walkfilesmatching(fn::Function, dir::String, regexp::Regex; recursive = true)
    for (root, dirs, files) in walkdir(dir)
        if recursive
            for dir in dirs
                walkfilesmatching(fn, joinpath(root, dir), regexp; recursive = recursive)
            end
        end
        for file in files
            filepath = joinpath(root, file)
            if occursin(regexp, filepath)
                fn(filepath)
            end
        end
    end
end

isnewer(filepath1, filepath2) = stat(filepath1).mtime > stat(filepath2).mtime

# Convert all pdf files in a dir to txt files. Iff recursive we traverse also into subdirs and
# iff onlyifnew we only convert if there is not already a text file with the same name
# which is older than the pdf file. Returns the paths to all text files.
function convert_pdf_files_to_txt(dirname::String; recursive = true, onlyifnew = true)
    count = converted = 0
    alltextfiles = String[]
    walkfilesmatching(dirname, r".*\.pdf$"; recursive = recursive) do pdffilepath
        count += 1
        txtfilepath = withtxtext(pdffilepath)
        push!(alltextfiles, txtfilepath)
        if onlyifnew == false || !isfile(txtfilepath) || isnewer(pdffilepath, txtfilepath)
            converted += 1
            pdf2text_apache_tika(pdffilepath)
        end
    end
    println("Found $count pdf files and converted $converted of them.")
    return alltextfiles
end


text:
module Tables

	export
		read_families,
		write_table,
		read_table

	using Taro
	Taro.init()

	using DataArrays
	using DataFrames
	using ZipFile

	#TODO Make some max i.e. end for range coding, like small x or whatever.
	function read_and_clean_table(file_name, sheet, range)
		t = DataArray(Taro.readxl(file_name, sheet, range, header=false));
		t = t[Bool[!all(isna(t[r, :])) for r in 1:size(t)[1]], :]
		t = t[:, Bool[!all(isna(t[:, c])) for c in 1:size(t)[2]]]
		return t
	end

	function parse_csv_file(f)
		# Bank table.
		bt = DataFrame(
			bank = [],
			pin_type = [],
			pin_num = [],
			byte_group_0_pin_num = [],
			byte_group_1_pin_num = [],
			byte_group_2_pin_num = [],
			byte_group_3_pin_num = [],
		)

		# Read pin stuff.
		ls = readlines(f)
		for l in ls[4:end-2]
			r = split(l, ',')
			bank = r[4]
			byte_group = r[3]
			pin_type = r[7]

			if !(bank in bt[:bank])
				push!(bt, [bank; pin_type; 0; 0;0;0;0])
			end
			i = findfirst(bt[:bank], bank)
			@assert i != 0

			if pin_type != "NA"
				bt[i, :pin_type] = pin_type
			end
			bt[i, :pin_num] += 1

			if byte_group != "NA"
				byte_group = parse(Int, byte_group)
			end
			if byte_group == 0
				bt[i, :byte_group_0_pin_num] += 1
			elseif byte_group == 1
				bt[i, :byte_group_1_pin_num] += 1
			elseif byte_group == 2
				bt[i, :byte_group_2_pin_num] += 1
			elseif byte_group == 3
				bt[i, :byte_group_3_pin_num] += 1
			end
		end

		return f.name, bt
	end

	function read_families()
		xls_files = []
		for (r, d, f) in walkdir("families/")
			for ff in f
				if endswith(ff, ".xls")
					push!(xls_files, joinpath(r, ff))
				end
			end
		end

		families = Dict{}()

		for xls_fn in xls_files
			name = basename(xls_fn)[1:end-length(".xls")]

			#TODO Headers.
			hs = read_and_clean_table(xls_fn, "Summary", "A1:Z100")
			h = hs[1, :][:]
			s = hs[2:end, :]

			c = size(s)[2]
			summary = DataFrame(Dict(zip(1:c, [s[:, i] for i in 1:c])))
			function prepare(s)
				s = replace(s, " ", "_")
				s = replace(s, "(", "_")
				s = replace(s, ")", "")
				s = replace(s, "/", "")
				if isnumber(s[1])
					s = "_" * s
				end
				return s
			end
			n = map(prepare, Vector(h))
			sn = Symbol[parse(nn) for nn in n]
			names!(summary.colindex, sn)

			devices = s[:, 1][:]

			p = read_and_clean_table(xls_fn, "Pins", "B1:AZ2")
			packages = p[!isna(p[:])]
			# Grouped by compatibility.
			grouped_packages = []
			for c in 1:size(p)[2]
				gp = p[:, c][:]
				gp = gp[!isna(gp)] # Remove NAs.
				push!(grouped_packages, gp)
			end

			sg = read_and_clean_table(xls_fn, "Memory_Speed", "C1:Z2")
			speed_grades = sg[!isna(sg[:])]
			simple_speed_grades = speed_grades[
				Bool[length(s) == 2 for s in speed_grades]
			]

			ds = read_and_clean_table(xls_fn, "Memory_Speed", "B5:Z10")
			# pin_type -> speed_grade -> ddr_speed
			ddr_speeds = Dict{String, Dict{String, Int}}()
			for r in 1:size(ds)[1]
				pin_type = ds[r, 1]
				d = Dict{String, Int}()
				@assert size(sg)[2] == size(ds)[2]-1
				for c2 in 1:size(sg)[2]
					s = Int(ds[r, c2+1])
					for r2 in 1:size(sg)[1]
						speed_grade = sg[r2, c2]
						if !isna(speed_grade)
							d[speed_grade] = s
						end
					end
				end
				ddr_speeds[pin_type] = d
			end

			pt = read_and_clean_table(xls_fn, "Pins", "B5:AZ5")
			pt_per_p = Int(length(pt)/length(grouped_packages))
			pin_types = pt[1:pt_per_p]

			combs = read_and_clean_table(xls_fn, "Pins", "B7:AZ30")

			# First column is device-package combination.
			# Others columns are pin count for pin_types, respectively.
			dev_pack_pins = DataFrame()
			dev_pack_pins[:dev_pack] = []
			for pt in pin_types
				dev_pack_pins[Symbol(pt)] = []
			end
			for (c, gp) in enumerate(grouped_packages)
				for (r, d) in enumerate(devices)
					pc = combs[r, (c-1)*pt_per_p+1:c*pt_per_p]
					pc[isna(pc)] = 0
					if sum(pc) > 0
						for p in gp
							dp = (d, p)
							push!(dev_pack_pins, [dp; pc])
						end
					end
				end
			end

			# Key is device-package combination, value is table of banks.
			zip_fn = xls_fn[1:end-length(".xls")] * ".zip"
			fn_banks = Dict()
			if isfile(zip_fn)
				z = ZipFile.Reader(zip_fn)
					for f in z.files
						if endswith(f.name, ".csv")
							fn, bt = parse_csv_file(f)
							fn_banks[fn] = bt
						end
					end
				close(z)
			end
			dev_pack_banks = Dict()
			for d in devices
				for gp in grouped_packages
					found = true

					for p in gp
						dp = (d, p)
						if dp in dev_pack_pins[:dev_pack]
							# If there is combination,
							# then need to have csv file.
							found = false
							fn = lowercase(dp[1] * dp[2] * "pkg.csv")
							if haskey(fn_banks, fn)
								for p in gp
									dp = (d, p)
									dev_pack_banks[dp] = fn_banks[fn]
								end
								delete!(fn_banks, fn)
								found = true
								break # Next gp
							end
						end
					end

					if !found
						warn("Do not have banks info for device: ", d,
							" packages: ", gp)
					end
				end
			end
			if length(fn_banks) != 0
				println("Not added dev-pack to bank data: ", keys(fn_banks))
			end

			tables = Dict(
				"summary" => summary,
				"devices" => devices,
				"speed_grades" => speed_grades,
				"simple_speed_grades" => simple_speed_grades,
				"ddr_speeds" => ddr_speeds,
				"packages" => packages,
				"grouped_packages" => grouped_packages,
				"pin_types" => pin_types,
				"dev_pack_pins" => dev_pack_pins,
				"dev_pack_banks" => dev_pack_banks
			)

			families[name] = tables

		end

		families
	end

	# Reading DataFrame per sheets to XLS with header as first row.
	function write_table(file_name, sheet, data_frame::DataFrame)
		w = Workbook()
		s = createSheet(w, sheet)

		r = createRow(s, 0)
		for ci in 1:size(data_frame)[2]
			c = createCell(r, ci-1)
			setCellValue(c, string(names(data_frame)[ci]))
		end

		for ri in 1:size(data_frame)[1]
			r = createRow(s, ri)
			for ci in 1:size(data_frame)[2]
				c = createCell(r, ci-1)
				setCellValue(c, data_frame[ri, ci])
			end
		end

		mkpath(dirname(file_name))

		write(file_name, w)
	end

	# Reading DataFrame per sheets to XLS with header as first row.
	function read_table(file_name, sheet)
		hf = read_and_clean_table(file_name, sheet, "A1:Z3000")
		h = hf[1, :][:]
		f = hf[2:end, :]

		c = size(f)[2]
		t = DataFrame(Dict(zip(1:c, [f[:, i] for i in 1:c])))
		n = Symbol[parse(hh) for hh in h]
		names!(t.colindex, n)

		return t
	end

end # Tables


text:
"""
    mutable struct GOCGenerator <: AbstractNodeLabel

Mutable structure descendant of AbstractNodeLabel

# Fields
- `busname::String`: identify the bus associated to the generator
- `id::String`: identify the generator for the bus `busname`
- `power_min::Complex128`: Smin = Pmin + im Qmin
- `power_max::Complex128`: Smax = Pmax + im Qmax
- `participation_factor::Float64`: coefficient to define active power coupling constraints
- `dict_obj_coeffs::SortedDict{Int64,Float64}`: dictionary for objective coefficients: degree => coeff
"""
mutable struct GOCGenerator <: AbstractNodeLabel
  busid::Int64
  id::Any
  power_min::Complex #Smin = Pmin + i Qmin
  power_max::Complex #Smax = Pmax + i Qmax
  participation_factor::Float64
  dict_obj_coeffs::SortedDict{Int64,Float64} # dict : degree => coeff
end

"""
    create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: GOCGenerator

Create variables for generator `elemid` at `bus`  in `bus_vars` if necessary.\n
If `elem_formulation == :NewVar`, create Sgen variable for the generator\n
If `elem_formulation == :GOCcoupling`, create Sgen variable for the generator, binary variables for comparison of voltages and create delta variable\n
Return nothing
"""
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: GOCGenerator
    if elem_formulation == :NewVar
        bus_vars[elemid] = Variable(variable_name("Sgen", bus, elemid, scenario), Complex)
    elseif elem_formulation == :GOCcoupling
        bus_vars[elemid] = Variable(variable_name("Sgen", bus, elemid, scenario), Complex)

        # Binary variables for coupling constraints
        Vbc_inf_Vsc = get_binInf_varname(basecase_scenario_name(), scenario, bus)
        bus_vars[Vbc_inf_Vsc] = Variable(Vbc_inf_Vsc, Bool)

        Vsc_inf_Vbc = get_binInf_varname(scenario, basecase_scenario_name(), bus)
        bus_vars[Vsc_inf_Vbc] = Variable(Vsc_inf_Vbc, Bool)

        Veq = get_binEq_varname(scenario, basecase_scenario_name(), bus)
        bus_vars[Veq] = Variable(Veq, Bool)

        # Scaling factor for current scenario
        delta_varname = get_delta_varname(scenario)
        bus_vars[delta_varname] = Variable(delta_varname, Real)
    end

    return
end

"""
    [cstrname, polynom, lb, ub] = Snodal(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: GOCGenerator

Return the polynomial contribution in power of generator `elemid` at `bus`(name, value, lower bound, upper bound). Will be used to construct power balance constraints in polynomial problem.\n
If `elem_formulation == :NbMinVar`, return generator bounds `["UNIT", Polynomial() ,Smin, Smax]`\n
If `elem_formulation == :NewVar` or `:GOCCoupling`, return -Sgen `["UNIT", -Sgen, 0, 0]`\n
Return no contribution `["", Polynomial(), 0, 0]` otherwise
"""
function Snodal(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: GOCGenerator
    cstrname = "UNIT"
    if elem_formulation == :NbMinVar
      lb = element.power_min
      ub = element.power_max
      return [cstrname, Polynomial(), lb, ub]
  elseif elem_formulation == :NewVar || elem_formulation == :GOCcoupling
      return [cstrname, Polynomial(-bus_vars[elemid]), 0, 0]
    end
    return ["", Polynomial(), 0, 0]
end

"""
    constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: GOCGenerator

Return all the constraints defined by generator `elemid` at `bus`. Will be used to construct constraints in polynomial problem.\n
If `elem_formulation == :NewVar`, return generator bounds : "Genbounds" => Smin <= Sgen <= Smax\n
If `elem_formulation == :GOCCoupling`, return generator bounds defined by coupling constraints\n
Return empty dictionary otherwise
"""
function constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: GOCGenerator
    cstrs = SortedDict{String, Constraint}()
    if elem_formulation == :NewVar
      lb = element.power_min
      ub = element.power_max
      cstrname = get_GenBounds_cstrname()
      cstrs[cstrname] = lb << bus_vars[elemid] << ub
    end

    ## Active / Reactive formulation
    if elem_formulation == :GOCcoupling
        # Problem variables and parameters
        Sgen = bus_vars[elemid]

        Volt_sc = bus_vars[volt_name()]
        Volt_bc = OPFpbs[basecase_scenario_name()].mp.node_vars[bus][volt_name()]

        Vbc_inf_Vsc = bus_vars[get_binInf_varname(basecase_scenario_name(), scenario, bus)]
        Vsc_inf_Vbc = bus_vars[get_binInf_varname(scenario, basecase_scenario_name(), bus)]
        Veq = bus_vars[get_binEq_varname(scenario, basecase_scenario_name(), bus)]

        Pmin, Qmin = real(element.power_min), imag(element.power_min)
        Pmax, Qmax = real(element.power_max), imag(element.power_max)
        ϵ, M = get_GOC_Volt_ϵ(), get_GOC_BigM()

        ## Definition of binary variables, upper constraint
        ccname_bindef_upper = get_VoltBinDef_upper()
        cstrs[ccname_bindef_upper] = (abs2(Volt_bc) - abs2(Volt_sc) - (Vbc_inf_Vsc * (-ϵ) + (1-Vbc_inf_Vsc+Veq) * M - 2*M*Veq)) << 0
        # cstrs[ccname_bindef_upper] = (abs2(Volt_bc) - abs2(Volt_sc) - (Vbc_inf_Vsc * (-ϵ) + (1-Vbc_inf_Vsc) * M - M*(1-Vbc_inf_Vsc - Vsc_inf_Vbc))) << 0
        ## Definition of binary variables, lower constraint
        ccname_bindef_lower = get_VoltBinDef_lower()
        ##TODO: verify big M constraints
        cstrs[ccname_bindef_lower] = (abs2(Volt_bc) - abs2(Volt_sc) - (Vsc_inf_Vbc * ϵ + (1-Vsc_inf_Vbc+Veq) * (-M) + 2*M*Veq)) >> 0
        # cstrs[ccname_bindef_lower] = (abs2(Volt_bc) - abs2(Volt_sc) - (Vsc_inf_Vbc * ϵ + (1-Vsc_inf_Vbc) * (-M) + M*(1-Vbc_inf_Vsc - Vsc_inf_Vbc))) >> 0
        ####
        # If Vsc < Vbc, lhs = abs2(Volt_bc) - abs2(Volt_sc) > 0
        ## If Vbc_inf_Vsc = 1, positive lhs < -ϵ => upper constraint not satisfied => NOT POSSIBLE
        ## If Veq = 1, positive lhs <  0 => upper constraint not satisfied => NOT POSSIBLE
        ## If Vsc_inf_Vbc = 1, positive lhs < M, upper constraint satisfied AND lhs > ϵ, lower constraint satisfied

        # If Vbc < Vsc, lhs = abs2(Volt_bc) - abs2(Volt_sc) < 0
        ## If Vbc_inf_Vsc = 1, lhs < -ϵ, upper constraint satisfied AND lhs > -M, lower constraint satisfied
        ## If Veq = 1,  lhs < 0, upper constraint satisfied AND lhs > 0, lower constraint not satisfied => NOT POSSIBLE
        ## If Vsc_inf_Vbc = 1, lhs < M, upper constraint satisfied AND lhs > ϵ, lower constraint not satisfied => NOT POSSIBLE

        # If |Vbc-Vsc|< ϵ, -ϵ < abs2(Volt_bc) - abs2(Volt_sc) < ϵ
        ## If Vbc_inf_Vsc = 1, lhs < -ϵ, upper constraint not satisfied => NOT POSSIBLE
        ## If Veq = 1,  lhs < 0, upper constraint satisfied AND lhs > 0, lower constraint satisfied
        ## If Vsc_inf_Vbc = 1, lhs < M, upper constraint satisfied AND lhs > ϵ, lower constraint not satisfied => NOT POSSIBLE

        ####
        ## Definition of binary variables, complementarity
        ccname_compl = get_VoltBinDef_complement()
        ## Three separate possibilities => only one binary variable == 1
        cstrs[ccname_compl] = (Vbc_inf_Vsc + Vsc_inf_Vbc + Veq == 1)

        P0 = real(OPFpbs[basecase_scenario_name()].mp.node_vars[bus][elemid])
        Pgen, Qgen = real(Sgen), imag(Sgen)
        deltavar = bus_vars[get_delta_varname(scenario)]
        α = element.participation_factor

        ## Coupling constraint active :
        cstrname = get_CC_active_cstrname()
        cstrs[cstrname] = (Pgen - (P0 + α * deltavar)) == 0

        cstrname = get_GenBounds_cstrname()
        cstrs[cstrname] = Pmin << Pgen << Pmax

        ## Coupling constraint reactive : upper bound
        ccname_upper = get_CC_reactiveupper_cstrname()
        cstrs[ccname_upper] = ((Qgen - Vbc_inf_Vsc * (Qmin - Qmax) - Qmax)) << 0

        ## Coupling constraint reactive : lower bound
        ccname_lower = get_CC_reactivelower_cstrname()
        cstrs[ccname_lower] = ((Qgen - Vsc_inf_Vbc * (Qmax - Qmin) - Qmin)) >> 0
        #
        # cstrs["Bin_Vbc_inf_Vsc"] = (Vbc_inf_Vsc^2 - Vbc_inf_Vsc == 0)
        # cstrs["Bin_Vsc_inf_Vbc"] = (Vsc_inf_Vbc^2 - Vsc_inf_Vbc == 0)
        # cstrs["Bin_Veq"] = (Veq^2 - Veq == 0)
    end

    return cstrs
end

"""
    cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, Snode::Polynomial, lb, ub) where T <: GOCGenerator

Return the polynomial contribution in objective generator `elemid` at `bus`.
"""
function cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, Snode::Polynomial, lb, ub) where T <: GOCGenerator
  Sgen = gencost = Polynomial()
  if elem_formulation == :NbMinVar
    _, S, lb_gen, ub_gen = Snodal(element, bus, elemid, elem_formulation, bus_vars)
    abs((lb-lb_gen) - (ub-ub_gen)) ≤ 5e-5 || warn("Power missing at $bus... $(abs((lb-lb_gen) - (ub-ub_gen)))")
    Sgen = Snode - S - (ub-ub_gen)
    # Sgen = add(Snode, -S - (ub-ub_gen))
  elseif elem_formulation == :NewVar
    Sgen::Polynomial = bus_vars[elemid]
    # Sgen = Polynomial(bus_vars[elemid])
  end
  add!(Sgen, conj(Sgen))
  Sactive = (Sgen)*0.5
  Sexp = 1 + Polynomial()
  dict_coeffs = element.dict_obj_coeffs
  degrees_sorted = sort(collect(keys(dict_coeffs)))
  #quartic objective
  imax = 3
  #quadratic objective :TODO
  #imax = 2
  if imax != 3
    warning("Objective is troncated to $imax degree")
  end
  for i=1:imax
        degree = degrees_sorted[i] ## NOTE: NOT normal behaviour, required for linear cost...
        add!(gencost, dict_coeffs[degree]*Sexp)
        Sexp = Sexp*Sactive
  end
  return gencost
end