text:
#Naive convolution method using for loops with some parallel using
#@simd and @inbounds

function convolve(cLayer::Conv1D, Ai::AbstractArray{T1,3}) where {T1}

    Aip = padding(cLayer, Ai)

    f_H = cLayer.f
    s_H = cLayer.s
    lastDim = ndims(Ai)
    paddedS = paddedSize(cLayer, Ai)[1:end-2]
    s = cLayer.s
    f = cLayer.f
    outputS = outDims(cLayer, Ai)[1:end-2]
    ci, m = cLayer.inputS[end-1:end]
    co = cLayer.channels
    Z = zeros(eltype(Ai), outputS..., co, m)
    n_H, c, m = size(Z)
    W = cLayer.W
    B = cLayer.B
    # @simd
    for mi = 1:m
        # @simd for
        Threads.@threads for ci = 1:c
            # @simd
            for hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                @inbounds ai = view(Aip,h_start:h_end, :, mi)
                @inbounds Z[hi, ci, mi] = W[:, :, ci] ⋅ ai + B[1,1,ci]
                # ai = nothing
                # Base.GC.gc()
            end #for mi=1:m, hi=1:n_H
        end #for ci=1:c
    end #for mi=1:m
    # @inbounds Z .+= B[:, 1, :]
    actFun = cLayer.actFun
    Ao = eval(:($actFun($Z)))

    # Base.GC.gc()

    return Dict(:Z => Z, :A => Ao)
end #function convolve(cLayer::Conv1D

function convolve(cLayer::Conv2D, Ai::AbstractArray{T1,4}) where {T1}

    Aip = padding(cLayer, Ai)

    f_H, f_W = cLayer.f
    s_H, s_W = cLayer.s
    lastDim = ndims(Ai)
    cLayer.inputS = size(Ai)
    paddedS = paddedSize(cLayer, Ai)[1:end-2]
    s = cLayer.s
    f = cLayer.f
    outputS = outDims(cLayer, Ai)[1:end-2]
    ci, m = cLayer.inputS[end-1:end]
    co = cLayer.channels
    Z = zeros(eltype(Ai), outputS..., co, m)
    n_H, n_W, c, m = size(Z)
    W = cLayer.W
    B = cLayer.B
    # @simd
    for mi = 1:m
        # @simd
        Threads.@threads for ci = 1:c
            # @simd for
            # wi = 1:n_W
            for wi = 1:n_W, hi = 1:n_H
                # @simd
                # for hi = 1:n_H
                    h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                    h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                    w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                    w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                    @inbounds ai = view(Aip,h_start:h_end, w_start:w_end, :, mi)
                    # @inbounds ai = Aip[h_start:h_end, w_start:w_end, :, mi]
                    @inbounds Z[hi, wi, ci, mi] = W[:, :, :, ci] ⋅ ai + B[1,1,1,ci]
                # end #for hi=1:n_H
            end #for wi=1:n_W, hi=1:n_H
        end #for ci=1:c
    end #for mi=1:m
    # @inbounds Z .+= B[:, :, 1, :]
    actFun = cLayer.actFun
    Ao = eval(:($actFun($Z)))

    # Base.GC.gc()

    return Dict(:Z => Z, :A => Ao)
end #function convolve(cLayer::Conv2D

function convolve(cLayer::Conv3D, Ai::AbstractArray{T1,5}) where {T1}

    Aip = padding(cLayer, Ai)

    f_H, f_W, f_D = cLayer.f
    s_H, s_W, s_D = cLayer.s
    lastDim = ndims(Ai)
    paddedS = paddedSize(cLayer, Ai)[1:end-2]
    s = cLayer.s
    f = cLayer.f
    outputS = outDims(cLayer, Ai)[1:end-2]
    ci, m = cLayer.inputS[end-1:end]
    co = cLayer.channels
    Z = zeros(eltype(Ai), outputS..., co, m)
    n_H, n_W, n_D, c, m = size(Z)
    W = cLayer.W
    B = cLayer.B
    # @simd
    for mi = 1:m
        # @simd for
        Threads.@threads for ci = 1:c
            # @simd for

                # @simd for

                    for di = 1:n_D, wi = 1:n_W, hi = 1:n_H
                        h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                        h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                        w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                        w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                        d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
                        d_end = di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
                        @inbounds ai = view(Aip,
                            h_start:h_end,
                            w_start:w_end,
                            d_start:d_end,
                            :,
                            mi,
                        )
                        @inbounds Z[hi, wi, di, ci, mi] =
                            W[:, :, :, :, ci] ⋅ ai + B[1,1,1,1,ci]
                        # ai = nothing
                        # Base.GC.gc()
                    end #for hi=1:n_H
        #         end #for wi=1:n_W
        #     end #for di=1:n_D
        end #for ci=1:c
    end #for mi=1:m
    # @inbounds Z .+= B[:, :, :, 1, :]
    actFun = cLayer.actFun
    Ao = eval(:($actFun($Z)))

    # Base.GC.gc()

    return Dict(:Z => Z, :A => Ao)
end #function convolve(cLayer::Conv3D

export convolve

### dconvolve

#Naive convolution method using for loops with some parallel using
#@simd and @inbounds

function dconvolve!(
    cLayer::Conv1D,
    Ai::AbstractArray{T1,3},
    dAi::AbstractArray{T2,3},
    dZ::AbstractArray{T3,3},
) where {T1,T2,T3}

    Aip = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    dAip = zeros(promote_type(eltype(dZ),eltype(Ai)), size(Aip))

    f_H = cLayer.f
    s_H = cLayer.s
    lastDim = ndims(Ai)
    n_H, c, m = size(dZ)
    W = cLayer.W
    cLayer.dW = similar(W)
    cLayer.dW .= 0
    B = cLayer.B
    cLayer.dB = similar(B)
    cLayer.dB .= 0
    # @simd
    for mi = 1:m
        # @simd
        Threads.@threads for ci = 1:c
            for hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                @inbounds ai = view(Aip, h_start:h_end, :, mi)
                @inbounds dAip[h_start:h_end, :, mi] .+=
                    dZ[hi, ci, mi] .* W[:, :, ci]

                @inbounds cLayer.dW[:, :, ci] .+= (ai .* dZ[hi, ci, mi])
                @inbounds cLayer.dB[:, :, ci] .+= dZ[hi, ci, mi]
            end #for
        end #for ci=1:c
    end #for mi=1:m

    dAi .= dAip[1+padS[1]:end-padS[2], :, :]

    return dAi
end #function dconvolve(cLayer::Conv1D

function dconvolve!(
    cLayer::Conv2D,
    Ai::AbstractArray{T1,4},
    dAi::AbstractArray{T2,4},
    dZ::AbstractArray{T3,4},
) where {T1,T2,T3}

    Aip = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    dAip = zeros(promote_type(eltype(dZ),eltype(Ai)), size(Aip))

    f_H, f_W = cLayer.f
    s_H, s_W = cLayer.s
    lastDim = ndims(Ai)
    n_H, n_W, c, m = size(dZ)
    W = cLayer.W
    cLayer.dW = similar(W)
    cLayer.dW .= 0
    B = cLayer.B
    cLayer.dB = similar(B)
    cLayer.dB .= 0
    # @simd
    for mi = 1:m
        # @simd
        Threads.@threads for ci = 1:c
            for wi = 1:n_W, hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                @inbounds ai = view(Aip, h_start:h_end, w_start:w_end, :, mi)
                @inbounds dAip[h_start:h_end, w_start:w_end, :, mi] .+=
                    dZ[hi, wi, ci, mi] .* W[:, :, :, ci]

                @inbounds cLayer.dW[:, :, :, ci] .+= (ai .* dZ[hi, wi, ci, mi])
                @inbounds cLayer.dB[:, :, :, ci] .+= dZ[hi, wi, ci, mi]
            end #for
        end #for ci=1:c
    end #for mi=1:m

    dAi .= dAip[1+padS[1]:end-padS[2], 1+padS[3]:end-padS[4], :, :]

    return dAi
end #function dconvolve(cLayer::Conv2D

function dconvolve!(
    cLayer::Conv3D,
    Ai::AbstractArray{T1,5},
    dAi::AbstractArray{T2,5},
    dZ::AbstractArray{T3,5},
) where {T1,T2,T3}

    Aip = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    dAip = zeros(promote_type(eltype(dZ),eltype(Ai)), size(Aip))

    f_H, f_W, f_D = cLayer.f
    s_H, s_W, s_D = cLayer.s
    lastDim = ndims(Ai)
    n_H, n_W, n_D, c, m = size(dZ)
    W = cLayer.W
    cLayer.dW = similar(W)
    cLayer.dW .= 0
    B = cLayer.B
    cLayer.dB = similar(B)
    cLayer.dB .= 0
    # @simd
    for mi = 1:m
        # @simd
        Threads.@threads for ci = 1:c
            for di = 1:n_D, wi = 1:n_W, hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
                d_end = di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
                @inbounds ai =
                    view(Aip, h_start:h_end, w_start:w_end, d_start:d_end, :, mi)
                @inbounds dAip[
                    h_start:h_end,
                    w_start:w_end,
                    d_start:d_end,
                    :,
                    mi,
                ] .+= dZ[hi, wi, di, ci, mi] .* W[:, :, :, :, ci]

                @inbounds cLayer.dW[:, :, :, :, ci] .+=
                    (ai .* dZ[hi, wi, ci, mi])
                @inbounds cLayer.dB[:, :, :, :, ci] .+= dZ[hi, wi, di, ci, mi]
            end #for
        end #for ci=1:c
    end #for mi=1:m

    dAi .= dAip[1+padS[1]:end-padS[2],
        1+padS[3]:end-padS[4],
        1+padS[5]:end-padS[6],
        :,
        :,
    ]

    return nothing
end #function dconvolve(cLayer::Conv3D

export dconvolve!


text:
using Test
using MPI
using PETScX
using LinearAlgebra: norm

MPI.Initialized() || MPI.Init()

@testset "mpivec tests" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)

    # local first and last value
    n0 = sum(10 .+ (0:mpirank)) - 10
    n1 = sum(10 .+ (0:(mpirank + 1))) - 11

    # global last value
    ne = sum(10 .+ (0:mpisize)) - 11
    exact_length = ne + 1

    # exact norm
    exact_norm = norm(0:ne)

    for petsclib in PETScX.petsclibs
        PetscScalar = PETScX.scalartype(petsclib)
        PetscInt = PETScX.inttype(petsclib)

        for test_version in 1:3
            PETScX.Initialize(petsclib)
            if test_version == 1 # Test creation from array
                julia_x = PetscScalar.(n0:n1)
                petsc_x = PETScX.Vec(petsclib, comm, julia_x)
            elseif test_version == 2 # Test creation from local size
                petsc_x = PETScX.Vec(petsclib, comm, n1 - n0 + 1)

                julia_x = PetscScalar.(n0:n1)
                inds = PetscInt.(n0:n1)
                PETScX.setvalues!(petsc_x, inds, julia_x, PETScX.INSERT_VALUES)

            else # Test creation from global size
                petsc_x = PETScX.Vec(
                    petsclib,
                    comm,
                    PETScX.PETSC_DECIDE;
                    global_length = exact_length,
                )

                inds = PetscInt.(PETScX.ownershiprange(petsc_x))
                julia_x = PetscScalar.(inds)
                PETScX.setvalues!(petsc_x, inds, julia_x, PETScX.INSERT_VALUES)
            end

            PETScX.assemblybegin!(petsc_x)
            PETScX.assemblyend!(petsc_x)

            @test length(petsc_x) == exact_length

            vec_norm = norm(petsc_x)
            @test exact_norm ≈ vec_norm

            # We cannot get a local form for this type
            @test_throws PETScX.PETScX_NoLocalForm PETScX.getlocalform(petsc_x)

            PETScX.with_unsafe_localarray!(petsc_x) do array
                @test array == julia_x
            end

            PETScX.destroy(petsc_x)

            PETScX.Finalize(petsclib)
        end
    end
end

@testset "mpivec with ghost tests" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)

    # local first and last value
    n0 = sum(10 .+ (0:mpirank)) - 10
    n1 = sum(10 .+ (0:(mpirank + 1))) - 11
    local_length = n1 - n0 + 1

    # global last value
    ne = sum(10 .+ (0:mpisize)) - 11
    exact_length = ne + 1

    # exact norm
    exact_norm = norm(0:ne)

    for petsclib in PETScX.petsclibs
        PetscScalar = PETScX.scalartype(petsclib)
        PetscInt = PETScX.inttype(petsclib)

        # two ghost to left and right
        ghost = PetscInt.([])
        if mpirank != 0
            ghost = PetscInt.([ghost..., n0 - 2, n0 - 1])
        end
        if mpirank != mpisize - 1
            ghost = PetscInt.([ghost..., n1 + 1, n1 + 2])
        end
        for with_array in (true, false)
            PETScX.Initialize(petsclib)

            # Test creation with array
            if with_array
                julia_x = PetscScalar.(n0:n1)
                julia_x_with_ghost =
                    zeros(PetscScalar, local_length + length(ghost))
                julia_x_with_ghost[1:local_length] = julia_x
                petsc_x = PETScX.Vec(petsclib, comm, julia_x_with_ghost, ghost)
            else
                petsc_x = PETScX.Vec(petsclib, comm, local_length, ghost)

                inds = PetscInt.(PETScX.ownershiprange(petsc_x))
                julia_x = PetscScalar.(inds)
                PETScX.setvalues!(petsc_x, inds, julia_x, PETScX.INSERT_VALUES)
            end

            PETScX.assemblybegin!(petsc_x)
            PETScX.assemblyend!(petsc_x)

            @test length(petsc_x) == exact_length

            vec_norm = norm(petsc_x)
            @test exact_norm ≈ vec_norm

            l_x = PETScX.getlocalform(petsc_x)
            @test length(l_x) == local_length + length(ghost)

            if length(ghost) > 0
                vals = zeros(PetscScalar, length(ghost))
                inds = PetscInt.(local_length - 1 .+ (1:length(ghost)))
                PETScX.getvalues!(vals, l_x, inds)
                @test !(vals == ghost)

                PETScX.ghostupdatebegin!(petsc_x)
                PETScX.ghostupdateend!(petsc_x)

                vals = zeros(PetscScalar, length(ghost))
                inds = PetscInt.(local_length - 1 .+ (1:length(ghost)))
                PETScX.getvalues!(vals, l_x, inds)
                @test vals == ghost
            end

            PETScX.restorelocalform!(l_x)

            PETScX.withlocalform(petsc_x) do l_vec
                @test length(l_vec) == local_length + length(ghost)
            end

            PETScX.with_unsafe_localarray!(petsc_x) do array
                @test array == julia_x
            end

            PETScX.withlocalform(petsc_x) do l_vec
                PETScX.with_unsafe_localarray!(l_vec) do array
                    @test array[1:local_length] == julia_x
                    @test array[(local_length + 1):end] == ghost
                end
            end
            PETScX.with_unsafe_localarray!(petsc_x; write = true) do array
                array[1:local_length] .= 1
            end
            PETScX.ghostupdatebegin!(petsc_x)
            PETScX.ghostupdateend!(petsc_x)
            PETScX.with_unsafe_localarray!(petsc_x; read = true) do array
                @test all(array .== 1)
            end

            PETScX.destroy(petsc_x)

            PETScX.Finalize(petsclib)
        end
    end
end

ext: jl


text:
"""
$(TYPEDEF)

Epigraph of perspective function of halved squared Euclidean norm (AKA rotated
second-order) cone of dimension `dim`.

    $(FUNCTIONNAME){T}(dim::Int)
"""
mutable struct EpiPerSquare{T <: Real} <: Cone{T}
    dim::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    sqrt_hess_prod_updated::Bool
    inv_sqrt_hess_prod_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}

    dist::T
    rtdist::T
    denom::T
    sqrt_hess_vec::Vector{T}
    inv_sqrt_hess_vec::Vector{T}

    function EpiPerSquare{T}(dim::Int) where {T <: Real}
        @assert dim >= 3
        cone = new{T}()
        cone.dim = dim
        return cone
    end
end

use_dual_barrier(::EpiPerSquare) = false

reset_data(cone::EpiPerSquare) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.sqrt_hess_prod_updated =
    cone.inv_sqrt_hess_prod_updated = false)

use_sqrt_hess_oracles(::Int, cone::EpiPerSquare) = true

get_nu(cone::EpiPerSquare) = 2

function set_initial_point!(arr::AbstractVector, cone::EpiPerSquare)
    @views arr[1:2] .= 1
    @views arr[3:end] .= 0
    return arr
end

# TODO refac with dual feas check
function update_feas(cone::EpiPerSquare{T}) where T
    @assert !cone.feas_updated
    u = cone.point[1]
    v = cone.point[2]

    if u > eps(T) && v > eps(T)
        @views w = cone.point[3:end]
        cone.dist = u * v - sum(abs2, w) / 2
        cone.is_feas = (cone.dist > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::EpiPerSquare{T}) where T
    u = cone.dual_point[1]
    v = cone.dual_point[2]

    if u > eps(T) && v > eps(T)
        @views w = cone.dual_point[3:end]
        return (u * v - sum(abs2, w) / 2 > eps(T))
    end

    return false
end

function update_grad(cone::EpiPerSquare)
    @assert cone.is_feas

    @. cone.grad = cone.point / cone.dist
    g2 = cone.grad[2]
    cone.grad[2] = -cone.grad[1]
    cone.grad[1] = -g2

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::EpiPerSquare)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data

    mul!(H, cone.grad, cone.grad')
    inv_dist = inv(cone.dist)
    @inbounds for j in 3:cone.dim
        H[j, j] += inv_dist
    end
    H[1, 2] -= inv_dist

    cone.hess_updated = true
    return cone.hess
end

function update_inv_hess(cone::EpiPerSquare)
    @assert cone.is_feas
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    Hi = cone.inv_hess.data

    mul!(Hi, cone.point, cone.point')
    @inbounds for j in 3:cone.dim
        Hi[j, j] += cone.dist
    end
    Hi[1, 2] -= cone.dist

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiPerSquare,
    )
    u = cone.point[1]
    v = cone.point[2]
    w = @view cone.point[3:end]

    @inbounds for j in 1:size(prod, 2)
        uj = arr[1, j]
        vj = arr[2, j]
        @views wj = arr[3:end, j]
        ga = (dot(w, wj) - v * uj - u * vj) / cone.dist
        prod[1, j] = -ga * v - vj
        prod[2, j] = -ga * u - uj
        @. @views prod[3:end, j] = ga * w + wj
    end
    @. prod /= cone.dist

    return prod
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiPerSquare,
    )
    @assert cone.is_feas

    @inbounds @views for j in 1:size(prod, 2)
        pa = dot(cone.point, arr[:, j])
        @. prod[:, j] = pa * cone.point
    end
    @. @views prod[1, :] -= cone.dist * arr[2, :]
    @. @views prod[2, :] -= cone.dist * arr[1, :]
    @. @views prod[3:end, :] += cone.dist * arr[3:end, :]

    return prod
end

function update_sqrt_hess_prod(cone::EpiPerSquare{T}) where T
    @assert cone.is_feas
    @assert !cone.sqrt_hess_prod_updated
    if !isdefined(cone, :sqrt_hess_vec)
        cone.sqrt_hess_vec = zeros(T, cone.dim)
    end

    rtdist = cone.rtdist = sqrt(cone.dist)
    cone.denom = 2 * rtdist + cone.point[1] + cone.point[2]
    vec = cone.sqrt_hess_vec
    @. @views vec[3:end] = cone.point[3:end] / rtdist
    vec[1] = -cone.point[2] / rtdist - 1
    vec[2] = -cone.point[1] / rtdist - 1

    cone.sqrt_hess_prod_updated = true
    return
end

function update_inv_sqrt_hess_prod(cone::EpiPerSquare{T}) where T
    @assert cone.is_feas
    @assert !cone.inv_sqrt_hess_prod_updated
    if !isdefined(cone, :inv_sqrt_hess_vec)
        cone.inv_sqrt_hess_vec = zeros(T, cone.dim)
    end

    rtdist = cone.rtdist = sqrt(cone.dist)
    cone.denom = 2 * rtdist + cone.point[1] + cone.point[2]
    vec = cone.inv_sqrt_hess_vec
    copyto!(vec, cone.point)
    vec[1:2] .+= rtdist

    cone.inv_sqrt_hess_prod_updated = true
    return
end

function sqrt_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiPerSquare{T},
    ) where {T <: Real}
    if !cone.sqrt_hess_prod_updated
        update_sqrt_hess_prod(cone)
    end
    vec = cone.sqrt_hess_vec
    rtdist = cone.rtdist

    @inbounds @views for j in 1:size(arr, 2)
        dotj = dot(vec, arr[:, j]) / cone.denom
        @. prod[:, j] = dotj * vec
    end
    @. @views prod[1, :] -= arr[2, :] / rtdist
    @. @views prod[2, :] -= arr[1, :] / rtdist
    @. @views prod[3:end, :] += arr[3:end, :] / rtdist

    return prod
end

function inv_sqrt_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiPerSquare{T},
    ) where {T <: Real}
    if !cone.inv_sqrt_hess_prod_updated
        update_inv_sqrt_hess_prod(cone)
    end
    vec = cone.inv_sqrt_hess_vec
    rtdist = cone.rtdist

    @inbounds @views for j in 1:size(arr, 2)
        dotj = dot(vec, arr[:, j]) / cone.denom
        @. prod[:, j] = dotj * vec
    end
    @. @views prod[1, :] -= arr[2, :] * rtdist
    @. @views prod[2, :] -= arr[1, :] * rtdist
    @. @views prod[3:end, :] += arr[3:end, :] * rtdist

    return prod
end

function dder3(cone::EpiPerSquare, dir::AbstractVector)
    @assert cone.grad_updated
    dim = cone.dim
    dder3 = cone.dder3
    point = cone.point
    u = point[1]
    v = point[2]
    u_dir = dir[1]
    v_dir = dir[2]
    @views w = point[3:end]
    @views w_dir = dir[3:end]

    jdotpd = u * v_dir + v * u_dir - dot(w, w_dir)
    hess_prod!(dder3, dir, cone)
    dotdHd = -dot(dir, dder3)
    dotpHd = dot(point, dder3)
    dder3 .*= jdotpd
    @. @views dder3[3:end] += dotdHd * w + dotpHd * w_dir
    dder3[1] += -dotdHd * v - dotpHd * v_dir
    dder3[2] += -dotdHd * u - dotpHd * u_dir
    dder3 ./= 2 * cone.dist

    return dder3
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

const _default_PT = Float32 # Default type/eltype for variate/parameter values

abstract type AbstractVarBounds end

function Base.intersect(a::AbstractVarBounds, b::AbstractVarBounds)
    totalndof(a) != totalndof(b) && throw(ArgumentError("Can't intersect bounds with different number of DOF"))
    unsafe_intersect(a, b)
end

@inline float_iseven(n::T) where {T<:AbstractFloat} = (n - T(2) * floor((n + T(0.5)) * T(0.5))) < T(0.5)

struct NoVarBounds <: AbstractVarBounds
    ndims::Int
end

# TODO: Find better solution to determine value type if no bounds:
Base.eltype(bounds::NoVarBounds) = _default_PT

Base.in(x::Any, bounds::NoVarBounds) = true

Base.isinf(bounds::NoVarBounds) = true

ValueShapes.totalndof(b::NoVarBounds) = b.ndims

unsafe_intersect(a::NoVarBounds, b::NoVarBounds) = a
unsafe_intersect(a::AbstractVarBounds, b::NoVarBounds) = a
unsafe_intersect(a::NoVarBounds, b::AbstractVarBounds) = b

abstract type VarVolumeBounds{T<:Real, V<:SpatialVolume{T}} <: AbstractVarBounds end

Base.in(x::Any, bounds::VarVolumeBounds) = in(x, bounds.vol)

# Random.rand(rng::AbstractRNG, bounds::VarVolumeBounds) =
#     rand!(rng, bounds, Vector{float(eltype(bounds))}(totalndof(bounds)))

# Random.rand(rng::AbstractRNG, bounds::VarVolumeBounds, n::Integer) =
#     rand!(rng, bounds, Matrix{float(eltype(bounds))}(npartotalndofams(bounds), n))
#
# Random.rand!(rng::AbstractRNG, bounds::VarVolumeBounds, x::StridedVecOrMat{<:Real}) = rand!(rng, spatialvolume(bounds), x)

ValueShapes.totalndof(b::VarVolumeBounds) = ndims(b.vol)

"""
    spatialvolume(b::VarVolumeBounds)::SpatialVolume

*BAT-internal, not part of stable public API.*

Returns the spatial volume that defines the variate/parameter bounds.
"""
function spatialvolume end

struct HyperRectBounds{T<:Real} <: VarVolumeBounds{T, HyperRectVolume{T}}
    vol::HyperRectVolume{T}

    function HyperRectBounds{T}(vol::HyperRectVolume{T}) where {T<:Real}
        nd = ndims(vol)
        (nd > 0) && isempty(vol) && throw(ArgumentError("Cannot create bounds with emtpy volume of $nd dimensions"))
        new{T}(vol)
    end
end

HyperRectBounds(vol::HyperRectVolume{T}) where {T<:Real} = HyperRectBounds{T}(vol)
HyperRectBounds(lo::AbstractVector{T}, hi::AbstractVector{T}) where {T<:Real} = HyperRectBounds(HyperRectVolume(lo, hi))
HyperRectBounds(intervals::AbstractVector{<:ClosedInterval{<:Real}}) = HyperRectBounds(minimum.(intervals), maximum.(intervals))
HyperRectBounds(bounds::AbstractInterval) = HyperRectBounds([bounds.left], [bounds.right])

function HyperRectBounds(bounds::Vector{<:AbstractInterval})
    lo = [b.left for b in bounds]
    hi = [b.right for b in bounds]

    return HyperRectBounds(lo, hi)
end

Base.similar(bounds::HyperRectBounds) = HyperRectBounds(
        HyperRectVolume(
            fill!(similar(bounds.vol.lo), zero(eltype(bounds.vol.lo))),
            fill!(similar(bounds.vol.hi), one(eltype(bounds.vol.hi)))
        )
    )

Base.eltype(bounds::HyperRectBounds{T}) where {T} = T

function Base.intersect(a::HyperRectBounds, b::HyperRectBounds)
    c = similar(a)
    for i in eachindex(a.vol.lo, a.vol.hi, b.vol.lo, b.vol.hi)
        iv_a = a.vol.lo[i]..a.vol.hi[i]
        iv_b = b.vol.lo[i]..b.vol.hi[i]
        iv_c = iv_a ∩ iv_b
        c.vol.lo[i] = minimum(iv_c)
        c.vol.hi[i] = maximum(iv_c)
    end
    c
end

Base.isinf(bounds::HyperRectBounds) = isinf(bounds.vol)

function Base.vcat(xs::HyperRectBounds...)
    lo = vcat(map(x -> x.vol.lo, xs)...)
    hi = vcat(map(x -> x.vol.hi, xs)...)
    HyperRectBounds(lo, hi)
end

spatialvolume(bounds::HyperRectBounds) = bounds.vol


text:
struct CLTestStruct
    f1::NTuple{3, Float32}
    f2::Nothing
    f3::Float32
end

@testset "OpenCL.Kernel" begin

    test_source = "
    __kernel void sum(__global const float *a,
                      __global const float *b,
                      __global float *c,
                      const unsigned int count)
    {
      unsigned int gid = get_global_id(0);
      if (gid < count) {
          c[gid] = a[gid] + b[gid];
      }
    }
    "

    #TODO: tests for invalid kernel build error && logs...

    @testset "OpenCL.Kernel constructor" begin
        for device in cl.devices()
            if device[:platform][:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Kernel constructor for " *
                     "Portable Computing Language Platform")
                continue
            end
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source)
            @test_throws ArgumentError cl.Kernel(prg, "sum")
            cl.build!(prg)
            @test cl.Kernel(prg, "sum") != nothing
        end
    end

    @testset "OpenCL.Kernel info" begin
        for device in cl.devices()
            if device[:platform][:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Kernel info for Portable Computing Language Platform")
                continue
            end
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source)
            cl.build!(prg)
            k = cl.Kernel(prg, "sum")
            @test k[:name] == "sum"
            @test k[:num_args] == 4
            @test k[:reference_count] > 0
            @test k[:program] == prg
            @test typeof(k[:attributes]) == String
        end
    end

    @testset "OpenCL.Kernel mem/workgroup size" begin
        for device in cl.devices()
            if device[:platform][:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Kernel mem/workgroup size for Portable Computing Language Platform")
                continue
            end
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source)
            cl.build!(prg)
            k = cl.Kernel(prg, "sum")
            for (sf, clf) in [(:size, cl.CL_KERNEL_WORK_GROUP_SIZE),
                              (:compile_size, cl.CL_KERNEL_COMPILE_WORK_GROUP_SIZE),
                              (:local_mem_size, cl.CL_KERNEL_LOCAL_MEM_SIZE),
                              (:private_mem_size, cl.CL_KERNEL_PRIVATE_MEM_SIZE),
                              (:prefered_size_multiple, cl.CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE)]
                @test cl.work_group_info(k, sf, device) != nothing
                @test cl.work_group_info(k, clf, device) != nothing
                if sf != :compile_size
                    @test cl.work_group_info(k, sf, device) == cl.work_group_info(k, clf, device)
                end
            end
        end
    end

    @testset "OpenCL.Kernel set_arg!/set_args!" begin
         for device in cl.devices()

            if device[:platform][:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Kernel mem/workgroup size for Portable Computing Language Platform")
                continue
            end

            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)

            prg = cl.Program(ctx, source=test_source) |> cl.build!
            k = cl.Kernel(prg, "sum")

            count  = 1024
            nbytes = count * sizeof(Float32)

            h_ones = ones(Float32, count)

            A = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_ones)
            B = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=h_ones)
            C = cl.Buffer(Float32, ctx, :w, count)

            # sizeof mem object for buffer in bytes
            @test sizeof(A) == nbytes
            @test sizeof(B) == nbytes
            @test sizeof(C) == nbytes

            # we use julia's index by one convention
            @test cl.set_arg!(k, 1, A) != nothing
            @test cl.set_arg!(k, 2, B) != nothing
            @test cl.set_arg!(k, 3, C) != nothing
            @test cl.set_arg!(k, 4, UInt32(count)) != nothing

            cl.enqueue_kernel(queue, k, count) |> cl.wait
            r = cl.read(queue, C)

            @test all(x -> x == 2.0, r)
            cl.flush(queue)

            # test set_args with new kernel
            k2 = cl.Kernel(prg, "sum")
            cl.set_args!(k2, A, B, C, UInt32(count))

            h_twos = fill(2f0, count)
            cl.copy!(queue, A, h_twos)
            cl.copy!(queue, B, h_twos)

            #TODO: check for ocl version, fill is opencl v1.2
            #cl.enqueue_fill(queue, A, 2f0)
            #cl.enqueue_fill(queue, B, 2f0)

            cl.enqueue_kernel(queue, k, count)
            cl.finish(queue)

            r = cl.read(queue, C)

            @test all(x -> x == 4.0, r)
        end
    end

    @testset "OpenCL.Kernel enqueue_kernel" begin
        for device in cl.devices()
            if device[:platform][:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Kernel mem/workgroup size for Portable Computing Language Platform")
                continue
            end

            simple_kernel = "
                __kernel void test(__global float *i) {
                    *i += 1;
                };"

            ctx = cl.Context(device)

            h_buff = Float32[1,]
            d_buff = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=h_buff)

            p = cl.Program(ctx, source=simple_kernel) |> cl.build!
            k = cl.Kernel(p, "test")
            q = cl.CmdQueue(ctx)

            # dimensions must be the same size
            @test_throws ArgumentError q(k, (1,), (1,1), d_buff)
            @test_throws ArgumentError q(k, (1,1), (1,), d_buff)

            # dimensions are bounded
            max_work_dim = device[:max_work_item_dims]
            bad = tuple([1 for _ in 1:(max_work_dim + 1)])
            @test_throws MethodError q(k, bad, d_buff)

            # devices have finite work sizes
            @test_throws MethodError q(k, (typemax(Int),), d_buff)

            # blocking call to kernel finishes cmd queue
            q(k, 1, 1, d_buff)

            r = cl.read(q, d_buff)
            @test r[1] == 2

            # alternative kernel call syntax
            k[q, (1,), (1,)](d_buff)
            r = cl.read(q, d_buff)
            @test r[1] == 3

            # enqueue task is an alias for calling
            # a kernel with a global/local size of 1
            evt = cl.enqueue_task(q, k)
            r = cl.read(q, d_buff)
            @test r[1] == 4
        end
    end

    test_source = "
    struct __attribute__((packed)) Test2{
        long f1;
        int __attribute__((aligned (8))) f2;
    };
    __kernel void structest(__global float *out, struct Test2 b){
        out[0] = b.f1;
        out[1] = b.f2;
    }
    "
    for device in cl.devices()
        if device[:platform][:name] == "Portable Computing Language"
            @warn("Skipping OpenCL.Kernel constructor for " *
                 "Portable Computing Language Platform")
            continue
        end
        Sys.isapple() && continue
        ctx = cl.Context(device)
        prg = cl.Program(ctx, source = test_source)
        queue = cl.CmdQueue(ctx)
        cl.build!(prg)
        structkernel = cl.Kernel(prg, "structest")
        out = cl.Buffer(Float32, ctx, :w, 2)
        bstruct = (1, Int32(4))
        structkernel[queue, (1,)](out, bstruct)
        r = cl.read(queue, out)
        @test r  == [1f0, 4f0]
    end

    test_source = "
    //packed
    struct __attribute__((packed)) Test{
        float3 f1;
        int f2; // empty type gets replaced with Int32 (no empty types allowed in OpenCL)
        // you might need to define the alignement of fields to match julia's layout
        float f3; // for the types used here the alignement matches though!
    };
    __kernel void structest(__global float *out, struct Test a){
        out[0] = a.f1.x;
        out[1] = a.f1.y;
        out[2] = a.f1.z;
        out[3] = a.f3;
    }
    "

    for device in cl.devices()
        if device[:platform][:name] == "Portable Computing Language"
            @warn("Skipping OpenCL.Kernel constructor for " *
                 "Portable Computing Language Platform")
            continue
        end
        ctx = cl.Context(device)
        prg = cl.Program(ctx, source = test_source)
        queue = cl.CmdQueue(ctx)
        cl.build!(prg)
        structkernel = cl.Kernel(prg, "structest")
        out = cl.Buffer(Float32, ctx, :w, 4)
        astruct = CLTestStruct((1f0, 2f0, 3f0), nothing, 22f0)
        structkernel[queue, (1,)](out, astruct)
        r = cl.read(queue, out)
        @test r == [1f0, 2f0, 3f0, 22f0]
    end
end


text:
function main()
	const s = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"
	N = 13
	bestval = 0
	for i = 1:length(s)-N

		n = 1
		for j = 1:N
			n *= parseint(s[i+j-1])
		end
		if n > bestval
			bestval = n
		end
	end
	println(bestval)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
## ExoplanetsSysSim/src/corbits.jl
## (c) 2015 Eric B. Ford

if VERSION >= v"0.7.0-"
   using Libdl
   using Pkg
   using CORBITS
   #global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[".",joinpath(Pkg.devdir(),"ExoplanetsSysSim/"),joinpath(Pkg.devdir(),"CORBITS/"),"/usr/local/lib"])  # WARNING: Assumes can find libcorbits.so
   global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[".",joinpath(dirname(pathof(ExoplanetsSysSim))),joinpath(dirname(pathof(CORBITS))),"/usr/local/lib"])  # WARNING: Assumes can find libcorbits.so
elseif VERSION >= v"0.4.0-"
   global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[".",joinpath(Pkg.dir(),"ExoplanetsSysSim/"),joinpath(Pkg.dir(),"CORBITS/"),"/usr/local/lib"])  # WARNING: Assumes can find libcorbits.so
else
   global const LIB_CORBITS = find_library(["libcorbits.so"],[".",joinpath(Pkg.devdir(),"ExoplanetsSysSim/"),joinpath(Pkg.devdir(),"CORBITS/"),"/usr/local/lib"])  # WARNING: Assumes can find libcorbits.so
end

#if !isdefined(is_corbits_initialized)
#  @windows_only println("# Warning:  Can't call external C libraries from Windows.  CORBITS won't work")
#end
#is_corbits_initialized = true

# Call CORBITS's function prob_of_transits_approx_arrays(a, r_star, r, e, Omega, omega, inc, use)
# Returns a Cdouble (aka Float64)
# For documentation see https://github.com/jbrakensiek/CORBITS
function prob_of_transits_approx(a::Vector{Cdouble},r_star::Cdouble,r::Vector{Cdouble}, e::Vector{Cdouble},
                                 Omega::Vector{Cdouble}, omega::Vector{Cdouble}, inc::Vector{Cdouble}, use::Vector{Cint})
  @assert(length(a) == length(r) == length(e) == length(Omega) == length(omega) == length(inc) >= length(use) )
  @assert(length(use) >=1 )
  num_pl = length(use)
  return ccall( (:prob_of_transits_approx_arrays, LIB_CORBITS), Cdouble,
      (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint),
       a, r_star, r, e, Omega, omega, inc, use, num_pl)
end

function test_corbits()
  #@compat is_windows() ? error("# CORBITS won't work on windows") : nothing
  a =  Cdouble[0.05, 0.15]
  r_star = convert(Cdouble,0.005)
  r = Cdouble[0.0001,0.0001]
  ecc = Cdouble[0.02, 0.1]
  Omega = Cdouble[0.0, 0.0]
  omega = Cdouble[ 0.0, 0.5]
  inc = Cdouble[pi/2, pi/2]
  use_pl = Cint[1,1]
  prob_of_transits_approx(a, r_star, r, ecc, Omega, omega, inc, use_pl)
end


text:
a = 2
function bar() println(a); if false a = 3 end end
# rien à voir avec le fait d'être exécuté ou non
# la présence *statique* de l'affectation a=3 en fait une variable locale
bar()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function τpeak(T::Array{Float64,1}, X::Array{Float64,1})::Float64

    # find the index of the max value -
    idx = argmax(X)
    return T[idx]
end

function τlag(rule::Function, T::Array{Float64,1}, X::Array{Float64,1})::Float64

    # filter -
    idx = findfirst(rule, X)

    if (isnothing(idx) == true)
        return 1000.0
    end

    # return -
    return T[idx]
end

function auc(T::Array{Float64,1}, X::Array{Float64,1})::Float64
    return integrate(T, X)
end

function max_FIIa_rate(T::Array{Float64,1}, X::Array{Float64,1})::Float64

    # intialize -
    N = 600 # sim to 6 min
    diff_array = Array{Float64,1}(undef, N - 1)
    ΔT = T[2] - T[1]

    # compute the ΔX/ΔT -
    for i ∈ 1:(N-1)

        # compute diff -
        ΔX = (X[i+1] - X[i])*(1/ΔT)
        diff_array[i] = ΔX
    end

    # find the max -
    return maximum(diff_array)
end

function model_output_vector(T::Array{Float64,1}, X::Array{Float64,1})::Array{Float64,1}

    # initialize -
    output_vector = Array{Float64,1}()

    # Lagtime,Peak,T.Peak,Max,AUC
    lagtime = τlag(x->x>10.0, T, X)
    FIIa_peak = maximum(X)
    time_to_peak = τpeak(T, X)
    FIIa_rate = max_FIIa_rate(T, X)
    area_under_curve = auc(T, X)

    # package -
    push!(output_vector, lagtime)
    push!(output_vector, FIIa_peak)
    push!(output_vector, time_to_peak)
    push!(output_vector, FIIa_rate)
    push!(output_vector, area_under_curve)

    # return -
    return output_vector
end

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
export inequalities, basicsemialgebraicset

struct BasicSemialgebraicSet{T, PT<:APL{T}, AT <: AbstractAlgebraicSet} <: AbstractBasicSemialgebraicSet
    V::AT
    p::Vector{PT}
end
function BasicSemialgebraicSet{T, PT}() where {T, PT<:APL{T}}
    BasicSemialgebraicSet(AlgebraicSet{T, PT}(), PT[])
end
function BasicSemialgebraicSet(V::AlgebraicSet{T, PT, A, S}, p::Vector{PT}) where {T, PT<:APL{T}, A, S<:AbstractAlgebraicSolver}
    BasicSemialgebraicSet{T, PT, typeof(V)}(V, p)
end
function BasicSemialgebraicSet(V::AlgebraicSet{T, PT, A, ST}, p::Vector{PS}) where {T, PT<:APL{T}, S, PS<:APL{S}, A, ST<:AbstractAlgebraicSolver}
    U = promote_type(T, S)
    PU = promote_type(PT, PS)
    BasicSemialgebraicSet(convert(AlgebraicSet{U, PU, A, ST}, V), Vector{PU}(p))
end
#BasicSemialgebraicSet{T, PT<:APL{T}}(V::AlgebraicSet{T, PT}, p::Vector{PT}) = BasicSemialgebraicSet{T, PT}(V, p)
function basicsemialgebraicset(V, p)
    BasicSemialgebraicSet(V, p)
end

MP.changecoefficienttype(::Type{BasicSemialgebraicSet{S, PS, AT}}, T::Type) where {S, PS, AT} = BasicSemialgebraicSet{T, MP.changecoefficienttype(PS, T), MP.changecoefficienttype(AT, T)}

function Base.convert(::Type{BasicSemialgebraicSet{T, PT, AT}}, set::BasicSemialgebraicSet) where {T, PT, AT}
    return BasicSemialgebraicSet{T, PT, AT}(set.V, set.p)
end

MP.variables(S::BasicSemialgebraicSet{T, PT, FullSpace}) where {T, PT<:APL{T}} = MP.variables(S.p)
MP.variables(S::BasicSemialgebraicSet) = sort(union(MP.variables(S.V), MP.variables(S.p)), rev=true)
nequalities(S::BasicSemialgebraicSet) = nequalities(S.V)
equalities(S::BasicSemialgebraicSet) = equalities(S.V)
addequality!(S::BasicSemialgebraicSet, p) = addequality!(S.V, p)
ninequalities(S::BasicSemialgebraicSet) = length(S.p)
inequalities(S::BasicSemialgebraicSet) = S.p
addinequality!(S::BasicSemialgebraicSet, p) = push!(S.p, p)

Base.intersect(S::BasicSemialgebraicSet, T::BasicSemialgebraicSet) = BasicSemialgebraicSet(S.V ∩ T.V, [S.p; T.p])
Base.intersect(S::BasicSemialgebraicSet, T::AbstractAlgebraicSet) = BasicSemialgebraicSet(S.V ∩ T, copy(S.p))
Base.intersect(S::BasicSemialgebraicSet, ::FullSpace) = S
Base.intersect(T::AbstractAlgebraicSet, S::BasicSemialgebraicSet) = intersect(S, T)

function Base.show(io::IO, V::BasicSemialgebraicSet)
    print(io, "{ (", join(variables(V), ", "), ") | ",
          join(string.(equalities(V)) .* " = 0", ", "))
    if nequalities(V) > 0
        print(io, ", ")
    end
    print(io, join(string.(inequalities(V)) .* " ≥ 0", ", "), " }")
end
function Base.show(io::IO, mime::MIME"text/plain", V::BasicSemialgebraicSet)
    print(io, "Basic semialgebraic Set defined by ")
    _show_els(io, "equalit", nequalities(V), equalities(V), "=")
    _show_els(io, "inequalit", ninequalities(V), inequalities(V), "≥")
end


text:
# random functions that dispatch either to CURAND or GPUArrays' generic RNG

using Random

export rand_logn!, rand_poisson!

# native RNG

"""
    CUDA.RNG()

A random number generator using `rand()` in a device kernel.

See also: [CUDA.Philox2x32](@ref)
"""
mutable struct RNG <: AbstractRNG
    seed::UInt32
    counter::UInt32

    function RNG(seed::Integer)
        new(seed%UInt32, 0)
    end
end

RNG() = RNG(Random.rand(UInt32))

function Random.seed!(rng::RNG, seed::Integer)
    rng.seed = seed % UInt32
    rng.counter = 0
end

Random.seed!(rng::RNG) = seed!(rng, Random.rand(UInt32))

function Random.rand!(rng::RNG, A::AnyCuArray)
    function kernel(A::AbstractArray{T}, seed::UInt32, counter::UInt32) where {T}
        device_rng = Random.default_rng()

        # initialize the state
        @inbounds Random.seed!(device_rng, seed, counter)

        # grid-stride loop
        threadId = threadIdx().x
        window = (blockDim().x - 1) * gridDim().x
        offset = (blockIdx().x - 1) * blockDim().x
        while offset < length(A)
            i = threadId + offset
            if i <= length(A)
                @inbounds A[i] = Random.rand(device_rng, T)
            end

            offset += window
        end

        return
    end

    kernel = @cuda launch=false name="rand!" kernel(A, rng.seed, rng.counter)
    config = launch_configuration(kernel.fun; max_threads=64)
    threads = max(32, min(config.threads, length(A)))
    blocks = min(config.blocks, cld(length(A), threads))
    kernel(A, rng.seed, rng.counter; threads=threads, blocks=blocks)

    new_counter = Int64(rng.counter) + length(A)
    overflow, remainder = fldmod(new_counter, typemax(UInt32))
    rng.seed += overflow     # XXX: is this OK?
    rng.counter = remainder

    A
end

function Random.randn!(rng::RNG, A::AnyCuArray{<:T}) where {T<:Union{AbstractFloat,Complex{<:AbstractFloat}}}
    function kernel(A::AbstractArray{T}, seed::UInt32, counter::UInt32)
        device_rng = Random.default_rng()

        # initialize the state
        @inbounds Random.seed!(device_rng, seed, counter)

        # grid-stride loop
        threadId = threadIdx().x
        window = (blockDim().x - 1) * gridDim().x
        offset = (blockIdx().x - 1) * blockDim().x
        while offset < length(A)
            i = threadId + offset
            j = threadId + offset + window
            if i <= length(A)
                # Box–Muller transform
                U1 = Random.rand(device_rng, T)
                U2 = Random.rand(device_rng, T)
                Z0 = sqrt(T(-2.0)*log(U1))*cos(T(2pi)*U2)
                Z1 = sqrt(T(-2.0)*log(U1))*sin(T(2pi)*U2)
                @inbounds A[i] = Z0
                if j <= length(A)
                    @inbounds A[j] = Z1
                end
            end

            offset += 2*window
        end

        return
    end

    kernel = @cuda launch=false name="rand!" kernel(A, rng.seed, rng.counter)
    config = launch_configuration(kernel.fun; max_threads=64)
    threads = max(32, min(config.threads, length(A)÷2))
    blocks = min(config.blocks, cld(length(A)÷2, threads))
    kernel(A, rng.seed, rng.counter; threads=threads, blocks=blocks)

    new_counter = Int64(rng.counter) + length(A)
    overflow, remainder = fldmod(new_counter, typemax(UInt32))
    rng.seed += overflow     # XXX: is this OK?
    rng.counter = remainder

    A
end

function default_rng()
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{rng::RNG}
    states = get!(task_local_storage(), :RNG) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        # CUDA RNG objects are cheap, so we don't need to cache them
        (; rng=RNG())
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    return state.rng
end

# old native RNG

# we keep this for the GPUArrays.jl tests

const idle_gpuarray_rngs = HandleCache{CuContext,GPUArrays.RNG}()

function GPUArrays.default_rng(::Type{<:CuArray})
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{rng::GPUArrays.RNG}
    states = get!(task_local_storage(), :GPUArraysRNG) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        new_rng = pop!(idle_gpuarray_rngs, cuda.context) do
            N = attribute(cuda.device, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
            buf = CuArray{NTuple{4, UInt32}}(undef, N)
            GPUArrays.RNG(buf)
        end

        finalizer(current_task()) do task
            push!(idle_gpuarray_rngs, cuda.context, new_rng) do
                # no need to do anything, as the RNG is collected by its finalizer
            end
        end

        Random.seed!(new_rng)
        (; rng=new_rng)
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    return state.rng
end

# generic functionality

function Random.rand!(rng::Union{RNG,CURAND.RNG,GPUArrays.RNG}, A::AbstractArray{T}) where {T}
    B = CuArray{T}(undef, size(A))
    Random.rand!(rng, B)
    copyto!(A, B)
end

function Random.rand(rng::Union{RNG,CURAND.RNG,GPUArrays.RNG}, T::Type)
    assertscalar("scalar rand")
    A = CuArray{T}(undef, 1)
    Random.rand!(rng, A)
    A[]
end

# RNG-less interface

cuda_rng() = default_rng()
curand_rng() = CURAND.default_rng()

function seed!(seed=Base.rand(UInt64))
    Random.seed!(cuda_rng(), seed)
    Random.seed!(curand_rng(), seed)
end

# CURAND in-place
Random.rand!(A::CURAND.UniformArray) = Random.rand!(curand_rng(), A)
Random.randn!(A::CURAND.NormalArray; kwargs...) = Random.randn!(curand_rng(), A; kwargs...)
rand_logn!(A::CURAND.LognormalArray; kwargs...) = CURAND.rand_logn!(curand_rng(), A; kwargs...)
rand_poisson!(A::CURAND.PoissonArray; kwargs...) = CURAND.rand_poisson!(curand_rng(), A; kwargs...)

# CURAND out-of-place
rand(T::CURAND.UniformType, dims::Dims) = Random.rand(curand_rng(), T, dims)
randn(T::CURAND.NormalType, dims::Dims; kwargs...) = Random.randn(curand_rng(), T, dims; kwargs...)
rand_logn(T::CURAND.LognormalType, dims::Dims; kwargs...) = CURAND.rand_logn(curand_rng(), T, dims; kwargs...)
rand_poisson(T::CURAND.PoissonType, dims::Dims; kwargs...) = CURAND.rand_poisson(curand_rng(), T, dims; kwargs...)

# support all dimension specifications
rand(T::CURAND.UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), T, Dims((dim1, dims...)))
randn(T::CURAND.NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_logn(T::CURAND.LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_poisson(T::CURAND.PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), T, Dims((dim1, dims...)); kwargs...)

# native in-place
Random.rand!(A::AnyCuArray) = Random.rand!(cuda_rng(), A)
Random.randn!(A::AnyCuArray) = Random.randn!(cuda_rng(), A)
rand_logn!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating lognormally-distributed random numbers of type $(eltype(A))")
rand_poisson!(A::AnyCuArray; kwargs...) =
    error("CUDA.jl does not support generating Poisson-distributed random numbers of type $(eltype(A))")

# native out-of-place
rand(T::Type, dims::Dims) = Random.rand!(CuArray{T}(undef, dims...))
randn(T::Type, dims::Dims; kwargs...) = Random.randn!(CuArray{T}(undef, dims...); kwargs...)
rand_logn(T::Type, dims::Dims; kwargs...) = rand_logn!(CuArray{T}(undef, dims...); kwargs...)
rand_poisson(T::Type, dims::Dims; kwargs...) = rand_poisson!(CuArray{T}(undef, dims...); kwargs...)

# support all dimension specifications
rand(T::Type, dim1::Integer, dims::Integer...) =
    Random.rand!(CuArray{T}(undef, dim1, dims...))
randn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_logn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_poisson(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson!(CuArray{T}(undef, dim1, dims...); kwargs...)

# untyped out-of-place
rand(dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), Dims((dim1, dims...)))
randn(dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_logn(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_poisson(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), Dims((dim1, dims...)); kwargs...)


text:
import Random
using AbstractAlgebra

function test_params(rng, nvariables, exps, nterms,
                        npolys, grounds, coeffssize, orderings)

    for n in nvariables
        for e in exps
            for nt in nterms
                for np in npolys
                    for gr in grounds
                        for ord in orderings
                            for csz in coeffssize
                                set = Groebner.generate_set(
                                    n, e, nt, np, csz, rng, gr, ord
                                )
                                isempty(set) && continue
                                gb = Groebner.groebner(set)
                                @test Groebner._isgroebner_reference(gb, initial_gens=set)

                                if !Groebner._isgroebner_reference(gb, initial_gens=set)
                                    @error "Beda!" n e nt np gr ord
                                    println("Rng ", rng)
                                    println("Set ", set)
                                    println("Gb ", gb)
                                    error("exit tests")
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

@testset "f4 random stress tests" begin

    rng = Random.MersenneTwister(5)

    # :degrevlex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 4:5]
    nterms     = [1:1, 1:2, 4:5]
    npolys     = [1:1, 1:3, 4:5]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:degrevlex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for f4"
    test_params(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :lex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:lex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for f4"
    test_params(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :deglex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:deglex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for f4"
    test_params(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :degrevlex rational case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [QQ]
    coeffssize = [3, 1000, 2^31 - 1, 9223372036854775837]
    orderings  = [:degrevlex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) rational tests for f4"
    test_params(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :deglex rational case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [QQ]
    coeffssize = [3, 1000, 2^31 - 1, 9223372036854775837]
    orderings  = [:deglex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) rational tests for f4"
    test_params(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ProtoBuf.google.protobuf

const HTTPExtension_Verb = (;[
    Symbol("NONE") => Int32(0),
    Symbol("GET") => Int32(1),
    Symbol("HEAD") => Int32(2),
    Symbol("POST") => Int32(3),
    Symbol("PUT") => Int32(4),
    Symbol("DELETE") => Int32(5),
    Symbol("CONNECT") => Int32(6),
    Symbol("OPTIONS") => Int32(7),
    Symbol("TRACE") => Int32(8),
]...)

mutable struct HTTPExtension <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HTTPExtension(; kwargs...)
        obj = new(meta(HTTPExtension), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct HTTPExtension
const __meta_HTTPExtension = Ref{ProtoMeta}()
function meta(::Type{HTTPExtension})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HTTPExtension)
            __meta_HTTPExtension[] = target = ProtoMeta(HTTPExtension)
            allflds = Pair{Symbol,Union{Type,String}}[:verb => Int32, :querystring => AbstractString]
            meta(target, HTTPExtension, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_HTTPExtension[]
    end
end
function Base.getproperty(obj::HTTPExtension, name::Symbol)
    if name === :verb
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :querystring
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct InvokeRequest <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function InvokeRequest(; kwargs...)
        obj = new(meta(InvokeRequest), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct InvokeRequest
const __meta_InvokeRequest = Ref{ProtoMeta}()
function meta(::Type{InvokeRequest})
    ProtoBuf.metalock() do
        if !isassigned(__meta_InvokeRequest)
            __meta_InvokeRequest[] = target = ProtoMeta(InvokeRequest)
            allflds = Pair{Symbol,Union{Type,String}}[:method => AbstractString, :data => ProtoBuf.google.protobuf._Any, :content_type => AbstractString, :http_extension => HTTPExtension]
            meta(target, InvokeRequest, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_InvokeRequest[]
    end
end
function Base.getproperty(obj::InvokeRequest, name::Symbol)
    if name === :method
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :data
        return (obj.__protobuf_jl_internal_values[name])::ProtoBuf.google.protobuf._Any
    elseif name === :content_type
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :http_extension
        return (obj.__protobuf_jl_internal_values[name])::HTTPExtension
    else
        getfield(obj, name)
    end
end

mutable struct InvokeResponse <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function InvokeResponse(; kwargs...)
        obj = new(meta(InvokeResponse), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct InvokeResponse
const __meta_InvokeResponse = Ref{ProtoMeta}()
function meta(::Type{InvokeResponse})
    ProtoBuf.metalock() do
        if !isassigned(__meta_InvokeResponse)
            __meta_InvokeResponse[] = target = ProtoMeta(InvokeResponse)
            allflds = Pair{Symbol,Union{Type,String}}[:data => ProtoBuf.google.protobuf._Any, :content_type => AbstractString]
            meta(target, InvokeResponse, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_InvokeResponse[]
    end
end
function Base.getproperty(obj::InvokeResponse, name::Symbol)
    if name === :data
        return (obj.__protobuf_jl_internal_values[name])::ProtoBuf.google.protobuf._Any
    elseif name === :content_type
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct Etag <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Etag(; kwargs...)
        obj = new(meta(Etag), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Etag
const __meta_Etag = Ref{ProtoMeta}()
function meta(::Type{Etag})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Etag)
            __meta_Etag[] = target = ProtoMeta(Etag)
            allflds = Pair{Symbol,Union{Type,String}}[:value => AbstractString]
            meta(target, Etag, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Etag[]
    end
end
function Base.getproperty(obj::Etag, name::Symbol)
    if name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

const StateOptions_StateConcurrency = (;[
    Symbol("CONCURRENCY_UNSPECIFIED") => Int32(0),
    Symbol("CONCURRENCY_FIRST_WRITE") => Int32(1),
    Symbol("CONCURRENCY_LAST_WRITE") => Int32(2),
]...)

const StateOptions_StateConsistency = (;[
    Symbol("CONSISTENCY_UNSPECIFIED") => Int32(0),
    Symbol("CONSISTENCY_EVENTUAL") => Int32(1),
    Symbol("CONSISTENCY_STRONG") => Int32(2),
]...)

mutable struct StateOptions <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StateOptions(; kwargs...)
        obj = new(meta(StateOptions), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StateOptions
const __meta_StateOptions = Ref{ProtoMeta}()
function meta(::Type{StateOptions})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StateOptions)
            __meta_StateOptions[] = target = ProtoMeta(StateOptions)
            allflds = Pair{Symbol,Union{Type,String}}[:concurrency => Int32, :consistency => Int32]
            meta(target, StateOptions, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StateOptions[]
    end
end
function Base.getproperty(obj::StateOptions, name::Symbol)
    if name === :concurrency
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :consistency
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct StateItem_MetadataEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StateItem_MetadataEntry(; kwargs...)
        obj = new(meta(StateItem_MetadataEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StateItem_MetadataEntry (mapentry)
const __meta_StateItem_MetadataEntry = Ref{ProtoMeta}()
function meta(::Type{StateItem_MetadataEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StateItem_MetadataEntry)
            __meta_StateItem_MetadataEntry[] = target = ProtoMeta(StateItem_MetadataEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => AbstractString]
            meta(target, StateItem_MetadataEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StateItem_MetadataEntry[]
    end
end
function Base.getproperty(obj::StateItem_MetadataEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct StateItem <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StateItem(; kwargs...)
        obj = new(meta(StateItem), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StateItem
const __meta_StateItem = Ref{ProtoMeta}()
function meta(::Type{StateItem})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StateItem)
            __meta_StateItem[] = target = ProtoMeta(StateItem)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Vector{UInt8}, :etag => Etag, :metadata => Base.Dict{AbstractString,AbstractString}, :options => StateOptions]
            meta(target, StateItem, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StateItem[]
    end
end
function Base.getproperty(obj::StateItem, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :etag
        return (obj.__protobuf_jl_internal_values[name])::Etag
    elseif name === :metadata
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,AbstractString}
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::StateOptions
    else
        getfield(obj, name)
    end
end

export HTTPExtension_Verb, HTTPExtension, InvokeRequest, InvokeResponse, StateItem_MetadataEntry, StateItem, Etag, StateOptions_StateConcurrency, StateOptions_StateConsistency, StateOptions
# mapentries: "StateItem_MetadataEntry" => ("AbstractString", "AbstractString")


text:
@testset "▶ types/drawing2              " begin
    # Boxplot
    stats = randn(5, 3)
    bp = G.Boxplot(stats, 3)
    @test bp.stats == stats
    @test bp.nobj == 3
    @test length(bp.boxstyles) == 3
    @test bp.boxstyles[1] isa G.BoxplotStyle
    @test bp.horiz == false

    # Heatmap
    data = randn(5, 5)
    h = (heatmap(data)).drawing
    @test h.zmin == minimum(data)
    @test h.zmax == maximum(data)
    @test h.cmap == colormap("RdBu", 10)
    @test h.cmiss == c"white"
    @test h.transpose == false
end

@testset "▶ /drawing2                   " begin
    # -------
    # BOXPLOT
    # -------
    f = Figure()

    # error on empty
    @test_throws ArgumentError boxplot()

    nobj = 3
    X = randn(20, nobj)
    q00, q25, q50, q75, q100 = quantile(X[:,1], [.0,.25,.5,.75,1.0])
    iqr = q75-q25
    mu = Statistics.mean(X[:,1])
    wlow = q25 - 1.5*iqr
    whigh = q75 + 1.5*iqr
    dhb = boxplot(X)
    @test dhb.drawing isa G.Boxplot
    @test dhb.drawing.nobj == nobj
    @test dhb.drawing.stats[1, :] == [wlow,q25,q50,q75,whigh,mu]
    @test gca().xaxis.min == 0
    @test gca().xaxis.max == nobj+1
    @test gca().xaxis.ticks.places == collect(1:nobj)

    global axmin, axmax
    axmin = Inf
    axmax = -Inf
    for c in 1:nobj
        global axmin, axmax
        q00, q25, q50, q75, q100 = quantile(X[:, c], [.0, .25, .5, .75, 1.0])
        iqr = q75 - q25
        mu = Statistics.mean(X[:, c])
        wlow  = q25 - 1.5*iqr
        whigh = q75 + 1.5*iqr

        axmin_ = min(q00, wlow)
        axmin_ -= 0.5abs(axmin_)
        axmin_ < axmin && (axmin = axmin_)

        axmax_ = max(q100, whigh)
        axmax_ += 0.5abs(axmax_)
        axmax_ > axmax && (axmax = axmax_)
    end
    @test gca().yaxis.min == axmin
    @test gca().yaxis.max == axmax

    dhb = boxplot(X, whiskers=Inf, horiz=true)
    q00, q25, q50, q75, q100 = quantile(X[:,1], [.0,.25,.5,.75,1.0])
    mu = Statistics.mean(X[:,1])
    @test dhb.drawing.stats[1, :] == [q00,q25,q50,q75,q100,mu]
end

@testset "▶ set_prop/drawing2           " begin
end

@testset "▶ apply_gle/drawing2          " begin
    f = Figure()
    begin # BOXPLOT
        clf()
        X = [0.214554    0.204954  -0.106862;
             -1.96864     0.441998   0.862069;
              0.355096    1.34095    0.94196;
             -1.04456    -0.159273   0.542868;
              0.698542    0.68037    0.130701;
              0.136469    0.878217   1.19285;
              0.811287   -0.526169   0.0629558;
              0.498855    0.833732   0.668236;
             -0.0956764   1.04196    0.548804;
             -1.09126     1.27123   -1.58358;
              5.0        -5.0        6.0  ] # forces outlier

        boxplot(X; box_lw=0.1, box_cols="red", med_ls="--", med_cols="seagreen",
                   mean_show=true, mean_markers="diamond", out_show=true,
                   out_markers="+", out_msize=0.2, out_mcols="red")
        s = G.assemble_figure(gcf(), debug=true)
        for i ∈ 1:3
            q25,q50,q75 = quantile(X[:,i], [0.25,0.50,0.75])
            iqr = q75-q25
            μ = sum(X[:,i])/length(X[:,i])
            iqr15 = 1.5iqr
            isin(s, "draw bp_vert $i $(q25-iqr15) $q25 $q50 $q75 $(q75+iqr15) $μ")
            isin(s, "0.6 0.3 1 0.1 \"rgba(1.0,0.0,0.0,1.0)\"") # box
            isin(s, "9 0.0 \"rgba(0.18,0.545,0.341,1.0)\" 1 diamond") # med,mean
        end
        # outliers
        for i ∈ 4:6
            isin(s, "d$i marker plus msize 0.2 color rgba(1.0,0.0,0.0,1.0)")
        end

        X[1,1] = NaN
        @test_throws ArgumentError boxplot(X)

    end # > BOXPLOT

    begin # HEATMAP
        clf()
        X = randn(5, 5)
        heatmap(X)
        s = G.assemble_figure(gcf(), debug=true)
        isin(s, "d1=c0,c1 d2=c0,c2 d3=c0,c3 d4=c0,c4 d5=c0,c5")
        isin(s, "draw hm_")
        for i ∈ 1:5
            isin(s, " $i d$i 0.2 0.2")
        end
        @test_throws ArgumentError heatmap_ticks("xxx")

        cla()
        @test_throws ArgumentError heatmap_ticks("x") # empty

    end # > HEATMAP
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## File Operations (Libuv-based) ##

module Filesystem

const S_IRUSR = 0o400
const S_IWUSR = 0o200
const S_IXUSR = 0o100
const S_IRWXU = 0o700
const S_IRGRP = 0o040
const S_IWGRP = 0o020
const S_IXGRP = 0o010
const S_IRWXG = 0o070
const S_IROTH = 0o004
const S_IWOTH = 0o002
const S_IXOTH = 0o001
const S_IRWXO = 0o007

export File,
       StatStruct,
       # open,
       futime,
       write,
       JL_O_WRONLY,
       JL_O_RDONLY,
       JL_O_RDWR,
       JL_O_APPEND,
       JL_O_CREAT,
       JL_O_EXCL,
       JL_O_TRUNC,
       JL_O_TEMPORARY,
       JL_O_SHORT_LIVED,
       JL_O_SEQUENTIAL,
       JL_O_RANDOM,
       JL_O_NOCTTY,
       S_IRUSR, S_IWUSR, S_IXUSR, S_IRWXU,
       S_IRGRP, S_IWGRP, S_IXGRP, S_IRWXG,
       S_IROTH, S_IWOTH, S_IXOTH, S_IRWXO

import .Base:
    IOError, _UVError, _sizeof_uv_fs, check_open, close, eof, eventloop, fd, isopen,
    bytesavailable, position, read, read!, readavailable, seek, seekend, show,
    skip, stat, unsafe_read, unsafe_write, write, transcode, uv_error,
    rawhandle, OS_HANDLE, INVALID_OS_HANDLE, windowserror, filesize

import .Base.RefValue

if Sys.iswindows()
    import .Base: cwstring
end

# Average buffer size including null terminator for several filesystem operations.
# On Windows we use the MAX_PATH = 260 value on Win32.
const AVG_PATH = Sys.iswindows() ? 260 : 512

include("path.jl")
include("stat.jl")
include("file.jl")
include(string(length(Core.ARGS) >= 2 ? Core.ARGS[2] : "", "file_constants.jl"))  # include($BUILDROOT/base/file_constants.jl)

## Operations with File (fd) objects ##

abstract type AbstractFile <: IO end

mutable struct File <: AbstractFile
    open::Bool
    handle::OS_HANDLE
    File(fd::OS_HANDLE) = new(true, fd)
end
if OS_HANDLE !== RawFD
    File(fd::RawFD) = File(Libc._get_osfhandle(fd)) # TODO: calling close would now destroy the wrong handle
end

rawhandle(file::File) = file.handle

# Filesystem.open, not Base.open
function open(path::AbstractString, flags::Integer, mode::Integer=0)
    req = Libc.malloc(_sizeof_uv_fs)
    local handle
    try
        ret = ccall(:uv_fs_open, Int32,
                    (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Int32, Int32, Ptr{Cvoid}),
                    C_NULL, req, path, flags, mode, C_NULL)
        handle = ccall(:uv_fs_get_result, Cssize_t, (Ptr{Cvoid},), req)
        ccall(:uv_fs_req_cleanup, Cvoid, (Ptr{Cvoid},), req)
        ret < 0 && uv_error("open($(repr(path)), $flags, $mode)", ret)
    finally # conversion to Cstring could cause an exception
        Libc.free(req)
    end
    return File(OS_HANDLE(@static Sys.iswindows() ? Ptr{Cvoid}(handle) : Cint(handle)))
end

isopen(f::File) = f.open

function check_open(f::File)
    if !isopen(f)
        throw(ArgumentError("file is closed"))
    end
end

function close(f::File)
    if isopen(f)
        f.open = false
        err = ccall(:jl_fs_close, Int32, (OS_HANDLE,), f.handle)
        f.handle = INVALID_OS_HANDLE
        uv_error("close", err)
    end
    nothing
end

# sendfile is the most efficient way to copy from a file descriptor
function sendfile(dst::File, src::File, src_offset::Int64, bytes::Int)
    check_open(dst)
    check_open(src)
    while true
        result = ccall(:jl_fs_sendfile, Int32, (OS_HANDLE, OS_HANDLE, Int64, Csize_t),
                       src.handle, dst.handle, src_offset, bytes)
        uv_error("sendfile", result)
        nsent = result
        bytes -= nsent
        src_offset += nsent
        bytes <= 0 && break
    end
    nothing
end

function unsafe_write(f::File, buf::Ptr{UInt8}, len::UInt, offset::Int64=Int64(-1))
    check_open(f)
    err = ccall(:jl_fs_write, Int32, (OS_HANDLE, Ptr{UInt8}, Csize_t, Int64),
                f.handle, buf, len, offset)
    uv_error("write", err)
    return len
end

write(f::File, c::UInt8) = write(f, Ref{UInt8}(c))

function truncate(f::File, n::Integer)
    check_open(f)
    req = Libc.malloc(_sizeof_uv_fs)
    err = ccall(:uv_fs_ftruncate, Int32,
                (Ptr{Cvoid}, Ptr{Cvoid}, OS_HANDLE, Int64, Ptr{Cvoid}),
                C_NULL, req, f.handle, n, C_NULL)
    Libc.free(req)
    uv_error("ftruncate", err)
    return f
end

function futime(f::File, atime::Float64, mtime::Float64)
    check_open(f)
    req = Libc.malloc(_sizeof_uv_fs)
    err = ccall(:uv_fs_futime, Int32,
                (Ptr{Cvoid}, Ptr{Cvoid}, OS_HANDLE, Float64, Float64, Ptr{Cvoid}),
                C_NULL, req, f.handle, atime, mtime, C_NULL)
    Libc.free(req)
    uv_error("futime", err)
    return f
end

function read(f::File, ::Type{UInt8})
    check_open(f)
    ret = ccall(:jl_fs_read_byte, Int32, (OS_HANDLE,), f.handle)
    uv_error("read", ret)
    return ret % UInt8
end

function read(f::File, ::Type{Char})
    b0 = read(f, UInt8)
    l = 8 * (4 - leading_ones(b0))
    c = UInt32(b0) << 24
    if l < 24
        s = 16
        while s ≥ l && !eof(f)
            p = position(f)
            b = read(f, UInt8)
            if b & 0xc0 != 0x80
                seek(f, p)
                break
            end
            c |= UInt32(b) << s
            s -= 8
        end
    end
    return reinterpret(Char, c)
end

read(f::File, ::Type{T}) where {T<:AbstractChar} = T(read(f, Char)) # fallback

function unsafe_read(f::File, p::Ptr{UInt8}, nel::UInt)
    check_open(f)
    ret = ccall(:jl_fs_read, Int32, (OS_HANDLE, Ptr{Cvoid}, Csize_t),
                f.handle, p, nel)
    uv_error("read", ret)
    ret == nel || throw(EOFError())
    nothing
end

bytesavailable(f::File) = max(0, filesize(f) - position(f)) # position can be > filesize

eof(f::File) = bytesavailable(f) == 0

function readbytes!(f::File, b::Array{UInt8}, nb=length(b))
    nr = min(nb, bytesavailable(f))
    if length(b) < nr
        resize!(b, nr)
    end
    ret = ccall(:jl_fs_read, Int32, (OS_HANDLE, Ptr{Cvoid}, Csize_t),
                f.handle, b, nr)
    uv_error("read", ret)
    return ret
end
read(io::File) = read!(io, Base.StringVector(bytesavailable(io)))
readavailable(io::File) = read(io)
read(io::File, nb::Integer) = read!(io, Base.StringVector(min(nb, bytesavailable(io))))

const SEEK_SET = Int32(0)
const SEEK_CUR = Int32(1)
const SEEK_END = Int32(2)

function seek(f::File, n::Integer)
    ret = ccall(:jl_lseek, Int64, (OS_HANDLE, Int64, Int32), f.handle, n, SEEK_SET)
    ret == -1 && (@static Sys.iswindows() ? windowserror : systemerror)("seek")
    return f
end

function seekend(f::File)
    ret = ccall(:jl_lseek, Int64, (OS_HANDLE, Int64, Int32), f.handle, 0, SEEK_END)
    ret == -1 && (@static Sys.iswindows() ? windowserror : systemerror)("seekend")
    return f
end

function skip(f::File, n::Integer)
    ret = ccall(:jl_lseek, Int64, (OS_HANDLE, Int64, Int32), f.handle, n, SEEK_CUR)
    ret == -1 && (@static Sys.iswindows() ? windowserror : systemerror)("skip")
    return f
end

function position(f::File)
    check_open(f)
    ret = ccall(:jl_lseek, Int64, (OS_HANDLE, Int64, Int32), f.handle, 0, SEEK_CUR)
    ret == -1 && (@static Sys.iswindows() ? windowserror : systemerror)("lseek")
    return ret
end

fd(f::File) = f.handle
stat(f::File) = stat(f.handle)

end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using DifferentialEquations
using BenchmarkTools
using MRIgeneralizedBloch
using MRIgeneralizedBloch: apply_hamiltonian_gbloch!, hamiltonian_linear
using Test

## set variables
α = (.01:.01:1) * π
TRF = 500e-6 # s
ω1 = α ./ TRF # rad/s
m0s = 0.1
m0f = 1-m0s
R1f = 0.3 # 1/s
R1s = 2.0 # 1/s
R2f = 1 / 50e-3 # 1/s
T2s = 10e-6 # s
Rx = 70 # 1/s

## pre-calcualtions
G = interpolate_greens_function(greens_superlorentzian, 0, maximum(TRF) / T2s)
mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5)

print("Time to pre-compute saturation: ")
R2slT = @time precompute_R2sl()
R2sl = R2slT[1]

## init
u0_5D = [0,0,m0f,  m0s,1]
u0_6D = [0,0,m0f,0,m0s,1]

## simulate and plot magnetization at the end of RF-pulses with different flip angles
M_full = zeros(length(ω1), 4)
M_appx = similar(M_full)
# Threads.@threads
for i in eachindex(ω1)
    M_full[i,:] = solve(DDEProblem(apply_hamiltonian_gbloch!, u0_5D, mfun, (0.0, TRF), (ω1[i], 1, 0, m0s, R1f, R2f, Rx, R1s, T2s, G)), MethodOfSteps(DP8()))[end][1:4]
    u = exp(hamiltonian_linear(ω1[i], 1, 0, TRF, m0s, R1f, R2f, Rx, R1s, R2sl(TRF, α[i], 1, T2s))) * u0_6D
    M_appx[i,:] = u[[1:3;5]]
end

@test M_appx ≈ M_full rtol = 1e-4

## benchmark the different solvers (excecute one line at a time to provoke individual results to be printed)
print("Time to solve the full gene. Bloch IDE for 100us π-pulse:")
@btime solve(DDEProblem(MRIgeneralizedBloch.apply_hamiltonian_gbloch!, u0_5D, mfun, (0.0, TRF), (ω1[end], 1, 0, m0s, R1f, R2f, Rx, R1s, T2s, G)), MethodOfSteps(DP8()))

print("Time to solve the linear approximation for 100us π-pulse:")
@btime exp(hamiltonian_linear(ω1[end-1], 1, 0, TRF, m0s, R1f, R2f, Rx, R1s, R2sl(TRF, α[end-1], 1, T2s))) * u0_6D

## ##########################################################################################################
# Test gradients
#############################################################################################################
R2sl, dR2sldT2s, dR2sldB1, dR2sldω1, dR2sldTRF, dR2sldT2sdω1, dR2sldB1dω1, dR2sldT2sTRF, dR2sldB1dTRF = R2slT

B1 = 1
α = π/2
TRF = 400e-6
ω1 = α / TRF

rtol = 1e-4

Δ = 1e-9
@test dR2sldT2s(TRF, ω1*TRF, B1, T2s) ≈ (R2sl(TRF, ω1*TRF, B1, T2s + Δ) - R2sl(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol
@test dR2sldTRF(TRF, ω1*TRF, B1, T2s) ≈ (R2sl(TRF+Δ, ω1*(TRF+Δ), B1, T2s) - R2sl(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol
@test dR2sldT2sTRF(TRF, ω1*TRF, B1, T2s) ≈ (dR2sldT2s(TRF+Δ, ω1*(TRF+Δ), B1, T2s) - dR2sldT2s(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol
@test dR2sldB1dTRF(TRF, ω1*TRF, B1, T2s) ≈ (dR2sldB1(TRF+Δ, ω1*(TRF+Δ), B1, T2s) - dR2sldB1(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol

Δ = 1e-6
@test dR2sldB1(TRF, ω1*TRF, B1, T2s) ≈ (R2sl(TRF, ω1*TRF, B1 + Δ, T2s) - R2sl(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol
@test dR2sldω1(TRF, ω1*TRF, B1, T2s) ≈ (R2sl(TRF, (ω1+Δ)*TRF, B1, T2s) - R2sl(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol

Δ = 1e-2
@test dR2sldT2sdω1(TRF, ω1*TRF, B1, T2s) ≈ (dR2sldT2s(TRF, (ω1+Δ)*TRF, B1, T2s) - dR2sldT2s(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol
@test dR2sldB1dω1(TRF, ω1*TRF, B1, T2s) ≈ (dR2sldB1(TRF, (ω1+Δ)*TRF, B1, T2s) - dR2sldB1(TRF, ω1*TRF, B1, T2s)) / Δ rtol = rtol


text:
using LinearRegression
using Test, DataFrames
using StatsModels

include("test_sweep_operator.jl")
include("test_utilities.jl")
include("test_LinearRegression.jl")
include("test_cooksd.jl")
include("test_lessthanfullrank.jl")
include("test_noint.jl")
include("test_heteroscedasticity.jl")
include("test_kfold.jl")


text:
# # POVM simulation
# This notebook shows how we can check how much depolarizing noise a qubit positive operator-valued measure (POVM) can take before it becomes simulable by projective measurements. The general method is described in [arXiv:1609.06139](https://arxiv.org/abs/1609.06139). The question of simulability by projective measurements boils down to an SDP problem. Eq. (8) from the paper defines the noisy POVM that we obtain subjecting a POVM $\mathbf{M}$ to a depolarizing channel $\Phi_t$:
#
# $$
# \left[\Phi_t\left(\mathbf{M}\right)\right]_i := t M_i + (1-t)\frac{\mathrm{tr}(M_i)}{d} \mathbb{1}.
# $$
#
# If this visibility $t\in[0,1]$ is one, the POVM $\mathbf{M}$ is simulable.
#
# We will use Convex.jl to solve the SDP problem.

using Convex, SCS, LinearAlgebra
if VERSION < v"1.2.0-DEV.0"
    (I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))
     LinearAlgebra.diagm(v::AbstractVector) = diagm(0 => v)
end

# For the qubit case, a four outcome qubit POVM $\mathbf{M} \in\mathcal{P}(2,4)$ is simulable if and only if
#
# $M_{1}=N_{12}^{+}+N_{13}^{+}+N_{14}^{+},$
#
# $M_{2}=N_{12}^{-}+N_{23}^{+}+N_{24}^{+},$
#
# $M_{3}=N_{13}^{-}+N_{23}^{-}+N_{34}^{+},$
#
# $M_{4}=N_{14}^{-}+N_{24}^{-}+N_{34}^{-},$
#
# where Hermitian operators $N_{ij}^{\pm}$ satisfy $N_{ij}^{\pm}\geq0$ and $N_{ij}^{+}+N_{ij}^{-}=p_{ij}\mathbb{1}$, where $i<j$ , $i,j=1,2,3,4$ and $p_{ij}\geq0$ as well as $\sum_{i<j}p_{ij}=1$, that is, the $p_{ij}$ values form a probability vector. This forms an SDP feasibility problem, which we can rephrase as an optimization problem by adding depolarizing noise to the left-hand side of the above equations and maximizing the visibility $t$:
#
# $\max_{t\in[0,1]} t$
#
# such that
#
# $t\,M_{1}+(1-t)\,\mathrm{tr}(M_{1})\frac{\mathbb{1}}{2}=N_{12}^{+}+N_{13}^{+}+N_{14}^{+},$
#
# $t\,M_{2}+(1-t)\,\mathrm{tr}(M_{2})\frac{\mathbb{1}}{2}=N_{12}^{-}+N_{23}^{+}+N_{24}^{+},$
#
# $t\,M_{3}+(1-t)\,\mathrm{tr}(M_{3})\frac{\mathbb{1}}{2}=N_{13}^{-}+N_{23}^{-}+N_{34}^{+},$
#
# $t\,M_{4}+(1-t)\,\mathrm{tr}(M_{4})\frac{\mathbb{1}}{2}=N_{14}^{-}+N_{24}^{-}+N_{34}^{-}$.
#
# We organize these constraints in a function that takes a four-output qubit POVM as its argument:

function get_visibility(K)
    noise = real([tr(K[i])*I(2)/2 for i=1:size(K, 1)])
    P = [[ComplexVariable(2, 2) for i=1:2] for j=1:6]
    q = Variable(6, Positive())
    t = Variable(1, Positive())
    constraints = [P[i][j] in :SDP for i=1:6 for j=1:2]
    constraints += sum(q)==1
    constraints += t<=1
    constraints += [P[i][1]+P[i][2] == q[i]*I(2) for i=1:6]
    constraints += t*K[1] + (1-t)*noise[1] == P[1][1] + P[2][1] + P[3][1]
    constraints += t*K[2] + (1-t)*noise[2] == P[1][2] + P[4][1] + P[5][1]
    constraints += t*K[3] + (1-t)*noise[3] == P[2][2] + P[4][2] + P[6][1]
    constraints += t*K[4] + (1-t)*noise[4] == P[3][2] + P[5][2] + P[6][2]
    p = maximize(t, constraints)
    solve!(p, SCS.Optimizer(verbose=0))
    return p.optval
end

# We check this function using the tetrahedron measurement (see Appendix B in [arXiv:quant-ph/0702021](https://arxiv.org/abs/quant-ph/0702021)). This measurement is non-simulable, so we expect a value below one.

function dp(v)
    I(2) + v[1]*[0 1; 1 0] + v[2]*[0 -im; im 0] + v[3]*[1 0; 0 -1]
end
b = [ 1  1  1;
     -1 -1  1;
     -1  1 -1;
      1 -1 -1]/sqrt(3)
M = [dp(b[i, :]) for i=1:size(b,1)]/4;
get_visibility(M)

# This value matches the one [we obtained](https://github.com/peterwittek/ipython-notebooks/blob/master/Simulating_POVMs.ipynb) using [PICOS](http://picos.zib.de/).


text:
abstract type DataSet end

struct TrainDataSet <: DataSet
    rootdir::AbstractString
end

struct ValDataSet <: DataSet
    rootdir::AbstractString
end


text:
using Documenter, CI_Test

makedocs(;
    modules=[CI_Test],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/EdwardStables/CI_Test.jl/blob/{commit}{path}#L{line}",
    sitename="CI_Test.jl",
    authors="EdwardStables <edward.stables1198@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/EdwardStables/CI_Test.jl",
)