text:
# Automatically-generated headers for cudnn
#
# DO NOT EDIT THIS FILE DIRECTLY.

# Automatically generated using Clang.jl

const CUDNN_MAJOR = 7
const CUDNN_MINOR = 6
const CUDNN_PATCHLEVEL = 5
const CUDNN_VERSION = CUDNN_MAJOR * 1000 + CUDNN_MINOR * 100 + CUDNN_PATCHLEVEL
const CUDNN_DIM_MAX = 8
const CUDNN_LRN_MIN_N = 1
const CUDNN_LRN_MAX_N = 16
const CUDNN_LRN_MIN_K = 1.0e-5
const CUDNN_LRN_MIN_BETA = 0.01
const CUDNN_BN_MIN_EPSILON = 0.0
const CUDNN_SEQDATA_DIM_COUNT = 4
const CUDNN_ATTN_QUERYMAP_ALL_TO_ONE = 0
const CUDNN_ATTN_QUERYMAP_ONE_TO_ONE = UInt32(1) << 0
const CUDNN_ATTN_DISABLE_PROJ_BIASES = 0
const CUDNN_ATTN_ENABLE_PROJ_BIASES = UInt32(1) << 1
const CUDNN_ATTN_WKIND_COUNT = 8

@cenum cudnnSeverity_t::UInt32 begin
    CUDNN_SEV_FATAL = 0
    CUDNN_SEV_ERROR = 1
    CUDNN_SEV_WARNING = 2
    CUDNN_SEV_INFO = 3
end

const cudnnContext = Cvoid
const cudnnHandle_t = Ptr{cudnnContext}

@cenum cudnnStatus_t::UInt32 begin
    CUDNN_STATUS_SUCCESS = 0
    CUDNN_STATUS_NOT_INITIALIZED = 1
    CUDNN_STATUS_ALLOC_FAILED = 2
    CUDNN_STATUS_BAD_PARAM = 3
    CUDNN_STATUS_INTERNAL_ERROR = 4
    CUDNN_STATUS_INVALID_VALUE = 5
    CUDNN_STATUS_ARCH_MISMATCH = 6
    CUDNN_STATUS_MAPPING_ERROR = 7
    CUDNN_STATUS_EXECUTION_FAILED = 8
    CUDNN_STATUS_NOT_SUPPORTED = 9
    CUDNN_STATUS_LICENSE_ERROR = 10
    CUDNN_STATUS_RUNTIME_PREREQUISITE_MISSING = 11
    CUDNN_STATUS_RUNTIME_IN_PROGRESS = 12
    CUDNN_STATUS_RUNTIME_FP_OVERFLOW = 13
end

const cudnnRuntimeTag_t = Cvoid

@cenum cudnnErrQueryMode_t::UInt32 begin
    CUDNN_ERRQUERY_RAWCODE = 0
    CUDNN_ERRQUERY_NONBLOCKING = 1
    CUDNN_ERRQUERY_BLOCKING = 2
end

const cudnnTensorStruct = Cvoid
const cudnnTensorDescriptor_t = Ptr{cudnnTensorStruct}
const cudnnConvolutionStruct = Cvoid
const cudnnConvolutionDescriptor_t = Ptr{cudnnConvolutionStruct}
const cudnnPoolingStruct = Cvoid
const cudnnPoolingDescriptor_t = Ptr{cudnnPoolingStruct}
const cudnnFilterStruct = Cvoid
const cudnnFilterDescriptor_t = Ptr{cudnnFilterStruct}
const cudnnLRNStruct = Cvoid
const cudnnLRNDescriptor_t = Ptr{cudnnLRNStruct}
const cudnnActivationStruct = Cvoid
const cudnnActivationDescriptor_t = Ptr{cudnnActivationStruct}
const cudnnSpatialTransformerStruct = Cvoid
const cudnnSpatialTransformerDescriptor_t = Ptr{cudnnSpatialTransformerStruct}
const cudnnOpTensorStruct = Cvoid
const cudnnOpTensorDescriptor_t = Ptr{cudnnOpTensorStruct}
const cudnnReduceTensorStruct = Cvoid
const cudnnReduceTensorDescriptor_t = Ptr{cudnnReduceTensorStruct}
const cudnnCTCLossStruct = Cvoid
const cudnnCTCLossDescriptor_t = Ptr{cudnnCTCLossStruct}
const cudnnTensorTransformStruct = Cvoid
const cudnnTensorTransformDescriptor_t = Ptr{cudnnTensorTransformStruct}

@cenum cudnnDataType_t::UInt32 begin
    CUDNN_DATA_FLOAT = 0
    CUDNN_DATA_DOUBLE = 1
    CUDNN_DATA_HALF = 2
    CUDNN_DATA_INT8 = 3
    CUDNN_DATA_INT32 = 4
    CUDNN_DATA_INT8x4 = 5
    CUDNN_DATA_UINT8 = 6
    CUDNN_DATA_UINT8x4 = 7
    CUDNN_DATA_INT8x32 = 8
end

@cenum cudnnMathType_t::UInt32 begin
    CUDNN_DEFAULT_MATH = 0
    CUDNN_TENSOR_OP_MATH = 1
    CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION = 2
end

@cenum cudnnNanPropagation_t::UInt32 begin
    CUDNN_NOT_PROPAGATE_NAN = 0
    CUDNN_PROPAGATE_NAN = 1
end

@cenum cudnnDeterminism_t::UInt32 begin
    CUDNN_NON_DETERMINISTIC = 0
    CUDNN_DETERMINISTIC = 1
end

@cenum cudnnReorderType_t::UInt32 begin
    CUDNN_DEFAULT_REORDER = 0
    CUDNN_NO_REORDER = 1
end

@cenum cudnnTensorFormat_t::UInt32 begin
    CUDNN_TENSOR_NCHW = 0
    CUDNN_TENSOR_NHWC = 1
    CUDNN_TENSOR_NCHW_VECT_C = 2
end

@cenum cudnnFoldingDirection_t::UInt32 begin
    CUDNN_TRANSFORM_FOLD = 0
    CUDNN_TRANSFORM_UNFOLD = 1
end

@cenum cudnnOpTensorOp_t::UInt32 begin
    CUDNN_OP_TENSOR_ADD = 0
    CUDNN_OP_TENSOR_MUL = 1
    CUDNN_OP_TENSOR_MIN = 2
    CUDNN_OP_TENSOR_MAX = 3
    CUDNN_OP_TENSOR_SQRT = 4
    CUDNN_OP_TENSOR_NOT = 5
end

@cenum cudnnReduceTensorOp_t::UInt32 begin
    CUDNN_REDUCE_TENSOR_ADD = 0
    CUDNN_REDUCE_TENSOR_MUL = 1
    CUDNN_REDUCE_TENSOR_MIN = 2
    CUDNN_REDUCE_TENSOR_MAX = 3
    CUDNN_REDUCE_TENSOR_AMAX = 4
    CUDNN_REDUCE_TENSOR_AVG = 5
    CUDNN_REDUCE_TENSOR_NORM1 = 6
    CUDNN_REDUCE_TENSOR_NORM2 = 7
    CUDNN_REDUCE_TENSOR_MUL_NO_ZEROS = 8
end

@cenum cudnnReduceTensorIndices_t::UInt32 begin
    CUDNN_REDUCE_TENSOR_NO_INDICES = 0
    CUDNN_REDUCE_TENSOR_FLATTENED_INDICES = 1
end

@cenum cudnnIndicesType_t::UInt32 begin
    CUDNN_32BIT_INDICES = 0
    CUDNN_64BIT_INDICES = 1
    CUDNN_16BIT_INDICES = 2
    CUDNN_8BIT_INDICES = 3
end

@cenum cudnnConvolutionMode_t::UInt32 begin
    CUDNN_CONVOLUTION = 0
    CUDNN_CROSS_CORRELATION = 1
end

@cenum cudnnConvolutionFwdPreference_t::UInt32 begin
    CUDNN_CONVOLUTION_FWD_NO_WORKSPACE = 0
    CUDNN_CONVOLUTION_FWD_PREFER_FASTEST = 1
    CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT = 2
end

@cenum cudnnConvolutionFwdAlgo_t::UInt32 begin
    CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM = 0
    CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM = 1
    CUDNN_CONVOLUTION_FWD_ALGO_GEMM = 2
    CUDNN_CONVOLUTION_FWD_ALGO_DIRECT = 3
    CUDNN_CONVOLUTION_FWD_ALGO_FFT = 4
    CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING = 5
    CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD = 6
    CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED = 7
    CUDNN_CONVOLUTION_FWD_ALGO_COUNT = 8
end

struct cudnnConvolutionFwdAlgoPerf_t
    algo::cudnnConvolutionFwdAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Csize_t
    determinism::cudnnDeterminism_t
    mathType::cudnnMathType_t
    reserved::NTuple{3, Cint}
end

@cenum cudnnConvolutionBwdFilterPreference_t::UInt32 begin
    CUDNN_CONVOLUTION_BWD_FILTER_NO_WORKSPACE = 0
    CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST = 1
    CUDNN_CONVOLUTION_BWD_FILTER_SPECIFY_WORKSPACE_LIMIT = 2
end

@cenum cudnnConvolutionBwdFilterAlgo_t::UInt32 begin
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0 = 0
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1 = 1
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT = 2
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3 = 3
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD = 4
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED = 5
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING = 6
    CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT = 7
end

struct cudnnConvolutionBwdFilterAlgoPerf_t
    algo::cudnnConvolutionBwdFilterAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Csize_t
    determinism::cudnnDeterminism_t
    mathType::cudnnMathType_t
    reserved::NTuple{3, Cint}
end

@cenum cudnnConvolutionBwdDataPreference_t::UInt32 begin
    CUDNN_CONVOLUTION_BWD_DATA_NO_WORKSPACE = 0
    CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST = 1
    CUDNN_CONVOLUTION_BWD_DATA_SPECIFY_WORKSPACE_LIMIT = 2
end

@cenum cudnnConvolutionBwdDataAlgo_t::UInt32 begin
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_0 = 0
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_1 = 1
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT = 2
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING = 3
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD = 4
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED = 5
    CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT = 6
end

struct cudnnConvolutionBwdDataAlgoPerf_t
    algo::cudnnConvolutionBwdDataAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Csize_t
    determinism::cudnnDeterminism_t
    mathType::cudnnMathType_t
    reserved::NTuple{3, Cint}
end

@cenum cudnnSoftmaxAlgorithm_t::UInt32 begin
    CUDNN_SOFTMAX_FAST = 0
    CUDNN_SOFTMAX_ACCURATE = 1
    CUDNN_SOFTMAX_LOG = 2
end

@cenum cudnnSoftmaxMode_t::UInt32 begin
    CUDNN_SOFTMAX_MODE_INSTANCE = 0
    CUDNN_SOFTMAX_MODE_CHANNEL = 1
end

@cenum cudnnPoolingMode_t::UInt32 begin
    CUDNN_POOLING_MAX = 0
    CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING = 1
    CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING = 2
    CUDNN_POOLING_MAX_DETERMINISTIC = 3
end

@cenum cudnnActivationMode_t::UInt32 begin
    CUDNN_ACTIVATION_SIGMOID = 0
    CUDNN_ACTIVATION_RELU = 1
    CUDNN_ACTIVATION_TANH = 2
    CUDNN_ACTIVATION_CLIPPED_RELU = 3
    CUDNN_ACTIVATION_ELU = 4
    CUDNN_ACTIVATION_IDENTITY = 5
end

@cenum cudnnLRNMode_t::UInt32 begin
    CUDNN_LRN_CROSS_CHANNEL_DIM1 = 0
end

@cenum cudnnDivNormMode_t::UInt32 begin
    CUDNN_DIVNORM_PRECOMPUTED_MEANS = 0
end

@cenum cudnnBatchNormMode_t::UInt32 begin
    CUDNN_BATCHNORM_PER_ACTIVATION = 0
    CUDNN_BATCHNORM_SPATIAL = 1
    CUDNN_BATCHNORM_SPATIAL_PERSISTENT = 2
end

@cenum cudnnBatchNormOps_t::UInt32 begin
    CUDNN_BATCHNORM_OPS_BN = 0
    CUDNN_BATCHNORM_OPS_BN_ACTIVATION = 1
    CUDNN_BATCHNORM_OPS_BN_ADD_ACTIVATION = 2
end

@cenum cudnnSamplerType_t::UInt32 begin
    CUDNN_SAMPLER_BILINEAR = 0
end

const cudnnDropoutStruct = Cvoid
const cudnnDropoutDescriptor_t = Ptr{cudnnDropoutStruct}

@cenum cudnnRNNAlgo_t::UInt32 begin
    CUDNN_RNN_ALGO_STANDARD = 0
    CUDNN_RNN_ALGO_PERSIST_STATIC = 1
    CUDNN_RNN_ALGO_PERSIST_DYNAMIC = 2
    CUDNN_RNN_ALGO_COUNT = 3
end

@cenum cudnnRNNMode_t::UInt32 begin
    CUDNN_RNN_RELU = 0
    CUDNN_RNN_TANH = 1
    CUDNN_LSTM = 2
    CUDNN_GRU = 3
end

@cenum cudnnRNNBiasMode_t::UInt32 begin
    CUDNN_RNN_NO_BIAS = 0
    CUDNN_RNN_SINGLE_INP_BIAS = 1
    CUDNN_RNN_DOUBLE_BIAS = 2
    CUDNN_RNN_SINGLE_REC_BIAS = 3
end

@cenum cudnnDirectionMode_t::UInt32 begin
    CUDNN_UNIDIRECTIONAL = 0
    CUDNN_BIDIRECTIONAL = 1
end

@cenum cudnnRNNInputMode_t::UInt32 begin
    CUDNN_LINEAR_INPUT = 0
    CUDNN_SKIP_INPUT = 1
end

@cenum cudnnRNNClipMode_t::UInt32 begin
    CUDNN_RNN_CLIP_NONE = 0
    CUDNN_RNN_CLIP_MINMAX = 1
end

@cenum cudnnRNNDataLayout_t::UInt32 begin
    CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_UNPACKED = 0
    CUDNN_RNN_DATA_LAYOUT_SEQ_MAJOR_PACKED = 1
    CUDNN_RNN_DATA_LAYOUT_BATCH_MAJOR_UNPACKED = 2
end

@cenum cudnnRNNPaddingMode_t::UInt32 begin
    CUDNN_RNN_PADDED_IO_DISABLED = 0
    CUDNN_RNN_PADDED_IO_ENABLED = 1
end

const cudnnRNNStruct = Cvoid
const cudnnRNNDescriptor_t = Ptr{cudnnRNNStruct}
const cudnnPersistentRNNPlan = Cvoid
const cudnnPersistentRNNPlan_t = Ptr{cudnnPersistentRNNPlan}
const cudnnRNNDataStruct = Cvoid
const cudnnRNNDataDescriptor_t = Ptr{cudnnRNNDataStruct}
const cudnnAlgorithmStruct = Cvoid
const cudnnAlgorithmDescriptor_t = Ptr{cudnnAlgorithmStruct}
const cudnnAlgorithmPerformanceStruct = Cvoid
const cudnnAlgorithmPerformance_t = Ptr{cudnnAlgorithmPerformanceStruct}

@cenum cudnnSeqDataAxis_t::UInt32 begin
    CUDNN_SEQDATA_TIME_DIM = 0
    CUDNN_SEQDATA_BATCH_DIM = 1
    CUDNN_SEQDATA_BEAM_DIM = 2
    CUDNN_SEQDATA_VECT_DIM = 3
end

const cudnnSeqDataStruct = Cvoid
const cudnnSeqDataDescriptor_t = Ptr{cudnnSeqDataStruct}
const cudnnAttnQueryMap_t = UInt32
const cudnnAttnStruct = Cvoid
const cudnnAttnDescriptor_t = Ptr{cudnnAttnStruct}

@cenum cudnnMultiHeadAttnWeightKind_t::UInt32 begin
    CUDNN_MH_ATTN_Q_WEIGHTS = 0
    CUDNN_MH_ATTN_K_WEIGHTS = 1
    CUDNN_MH_ATTN_V_WEIGHTS = 2
    CUDNN_MH_ATTN_O_WEIGHTS = 3
    CUDNN_MH_ATTN_Q_BIASES = 4
    CUDNN_MH_ATTN_K_BIASES = 5
    CUDNN_MH_ATTN_V_BIASES = 6
    CUDNN_MH_ATTN_O_BIASES = 7
end

@cenum cudnnWgradMode_t::UInt32 begin
    CUDNN_WGRAD_MODE_ADD = 0
    CUDNN_WGRAD_MODE_SET = 1
end

@cenum cudnnCTCLossAlgo_t::UInt32 begin
    CUDNN_CTC_LOSS_ALGO_DETERMINISTIC = 0
    CUDNN_CTC_LOSS_ALGO_NON_DETERMINISTIC = 1
end

@cenum cudnnLossNormalizationMode_t::UInt32 begin
    CUDNN_LOSS_NORMALIZATION_NONE = 0
    CUDNN_LOSS_NORMALIZATION_SOFTMAX = 1
end

# FIXME: can't use such a union as the type in a ccall expression
#Algorithm = Union{cudnnConvolutionFwdAlgo_t, cudnnConvolutionBwdFilterAlgo_t, cudnnConvolutionBwdDataAlgo_t, cudnnRNNAlgo_t, cudnnCTCLossAlgo_t}
#struct cudnnAlgorithm_t
#    algo::Algorithm
#end
cudnnAlgorithm_t = Cint

struct cudnnDebug_t
    cudnn_version::UInt32
    cudnnStatus::cudnnStatus_t
    time_sec::UInt32
    time_usec::UInt32
    time_delta::UInt32
    handle::cudnnHandle_t
    stream::CUstream
    pid::Culonglong
    tid::Culonglong
    cudaDeviceId::Cint
    reserved::NTuple{15, Cint}
end

const cudnnCallback_t = Ptr{Cvoid}
const cudnnFusedOpsConstParamStruct = Cvoid
const cudnnFusedOpsConstParamPack_t = Ptr{cudnnFusedOpsConstParamStruct}
const cudnnFusedOpsVariantParamStruct = Cvoid
const cudnnFusedOpsVariantParamPack_t = Ptr{cudnnFusedOpsVariantParamStruct}
const cudnnFusedOpsPlanStruct = Cvoid
const cudnnFusedOpsPlan_t = Ptr{cudnnFusedOpsPlanStruct}

@cenum cudnnFusedOps_t::UInt32 begin
    CUDNN_FUSED_SCALE_BIAS_ACTIVATION_CONV_BNSTATS = 0
    CUDNN_FUSED_SCALE_BIAS_ACTIVATION_WGRAD = 1
    CUDNN_FUSED_BN_FINALIZE_STATISTICS_TRAINING = 2
    CUDNN_FUSED_BN_FINALIZE_STATISTICS_INFERENCE = 3
    CUDNN_FUSED_CONV_SCALE_BIAS_ADD_ACTIVATION = 4
    CUDNN_FUSED_SCALE_BIAS_ADD_ACTIVATION_GEN_BITMASK = 5
    CUDNN_FUSED_DACTIVATION_FORK_DBATCHNORM = 6
end

@cenum cudnnFusedOpsConstParamLabel_t::UInt32 begin
    CUDNN_PARAM_XDESC = 0
    CUDNN_PARAM_XDATA_PLACEHOLDER = 1
    CUDNN_PARAM_BN_MODE = 2
    CUDNN_PARAM_BN_EQSCALEBIAS_DESC = 3
    CUDNN_PARAM_BN_EQSCALE_PLACEHOLDER = 4
    CUDNN_PARAM_BN_EQBIAS_PLACEHOLDER = 5
    CUDNN_PARAM_ACTIVATION_DESC = 6
    CUDNN_PARAM_CONV_DESC = 7
    CUDNN_PARAM_WDESC = 8
    CUDNN_PARAM_WDATA_PLACEHOLDER = 9
    CUDNN_PARAM_DWDESC = 10
    CUDNN_PARAM_DWDATA_PLACEHOLDER = 11
    CUDNN_PARAM_YDESC = 12
    CUDNN_PARAM_YDATA_PLACEHOLDER = 13
    CUDNN_PARAM_DYDESC = 14
    CUDNN_PARAM_DYDATA_PLACEHOLDER = 15
    CUDNN_PARAM_YSTATS_DESC = 16
    CUDNN_PARAM_YSUM_PLACEHOLDER = 17
    CUDNN_PARAM_YSQSUM_PLACEHOLDER = 18
    CUDNN_PARAM_BN_SCALEBIAS_MEANVAR_DESC = 19
    CUDNN_PARAM_BN_SCALE_PLACEHOLDER = 20
    CUDNN_PARAM_BN_BIAS_PLACEHOLDER = 21
    CUDNN_PARAM_BN_SAVED_MEAN_PLACEHOLDER = 22
    CUDNN_PARAM_BN_SAVED_INVSTD_PLACEHOLDER = 23
    CUDNN_PARAM_BN_RUNNING_MEAN_PLACEHOLDER = 24
    CUDNN_PARAM_BN_RUNNING_VAR_PLACEHOLDER = 25
    CUDNN_PARAM_ZDESC = 26
    CUDNN_PARAM_ZDATA_PLACEHOLDER = 27
    CUDNN_PARAM_BN_Z_EQSCALEBIAS_DESC = 28
    CUDNN_PARAM_BN_Z_EQSCALE_PLACEHOLDER = 29
    CUDNN_PARAM_BN_Z_EQBIAS_PLACEHOLDER = 30
    CUDNN_PARAM_ACTIVATION_BITMASK_DESC = 31
    CUDNN_PARAM_ACTIVATION_BITMASK_PLACEHOLDER = 32
    CUDNN_PARAM_DXDESC = 33
    CUDNN_PARAM_DXDATA_PLACEHOLDER = 34
    CUDNN_PARAM_DZDESC = 35
    CUDNN_PARAM_DZDATA_PLACEHOLDER = 36
    CUDNN_PARAM_BN_DSCALE_PLACEHOLDER = 37
    CUDNN_PARAM_BN_DBIAS_PLACEHOLDER = 38
end

@cenum cudnnFusedOpsPointerPlaceHolder_t::UInt32 begin
    CUDNN_PTR_NULL = 0
    CUDNN_PTR_ELEM_ALIGNED = 1
    CUDNN_PTR_16B_ALIGNED = 2
end

@cenum cudnnFusedOpsVariantParamLabel_t::UInt32 begin
    CUDNN_PTR_XDATA = 0
    CUDNN_PTR_BN_EQSCALE = 1
    CUDNN_PTR_BN_EQBIAS = 2
    CUDNN_PTR_WDATA = 3
    CUDNN_PTR_DWDATA = 4
    CUDNN_PTR_YDATA = 5
    CUDNN_PTR_DYDATA = 6
    CUDNN_PTR_YSUM = 7
    CUDNN_PTR_YSQSUM = 8
    CUDNN_PTR_WORKSPACE = 9
    CUDNN_PTR_BN_SCALE = 10
    CUDNN_PTR_BN_BIAS = 11
    CUDNN_PTR_BN_SAVED_MEAN = 12
    CUDNN_PTR_BN_SAVED_INVSTD = 13
    CUDNN_PTR_BN_RUNNING_MEAN = 14
    CUDNN_PTR_BN_RUNNING_VAR = 15
    CUDNN_PTR_ZDATA = 16
    CUDNN_PTR_BN_Z_EQSCALE = 17
    CUDNN_PTR_BN_Z_EQBIAS = 18
    CUDNN_PTR_ACTIVATION_BITMASK = 19
    CUDNN_PTR_DXDATA = 20
    CUDNN_PTR_DZDATA = 21
    CUDNN_PTR_BN_DSCALE = 22
    CUDNN_PTR_BN_DBIAS = 23
    CUDNN_SCALAR_SIZE_T_WORKSPACE_SIZE_IN_BYTES = 100
    CUDNN_SCALAR_INT64_T_BN_ACCUMULATION_COUNT = 101
    CUDNN_SCALAR_DOUBLE_BN_EXP_AVG_FACTOR = 102
    CUDNN_SCALAR_DOUBLE_BN_EPSILON = 103
end


text:
using SparseArrays
using DelimitedFiles: readdlm
using ProgressBars
using PyPlot

function fast_kcore(nodes::Vector{Int64}, edges::Matrix{Int64})
    n = length(nodes)
    adjacency_list = Dict{Int64, Vector{Int64}}()
    for i in 1:n
        adjacency_list[i] = Vector{Int64}()
    end
    for edge in eachrow(edges)
         neighbors = get!(Vector{Int64}, adjacency_list, edge[1])
         push!(neighbors, edge[2])
         neighbors = get!(Vector{Int64}, adjacency_list, edge[2])
         push!(neighbors, edge[1])
    end

    degrees = zeros(Int64, n)

    for edge in eachrow(edges)
        degrees[edge[1]] += 1
        degrees[edge[2]] += 1
    end

    max_degree = maximum(degrees)
    eta = Dict{Int64, Int64}()
    c = 0
    for i in ProgressBar(1:n)
        v = argmin(degrees)
        c = max(c, degrees[v])
        eta[v] = c
        degrees[v] = max_degree + 1
        for neighbor in adjacency_list[v]
            degrees[neighbor] -= 1
        end
    end

    return eta
end

filename = "com-amazon.ungraph.txt"
edges = readdlm(filename, '\t', Int, comments=true)
nodes = Vector(1:maximum(edges))
@time core_numbers = fast_kcore(nodes, edges)
println("Core value for ", filename, " = ", maximum(values(core_numbers)))

filename = "com-lj.ungraph.txt"
edges = readdlm(filename, '\t', Int, comments=true)
edges = edges .+ 1
nodes = Vector(1:maximum(edges))
@time core_numbers = fast_kcore(nodes, edges)
println("Core value for ", filename, " = ", maximum(values(core_numbers)))

filename = "scholar/net.txt"
edges = readdlm(filename, '\t', Int, comments=true)
edges = edges .+ 1
nodes = Vector(1:maximum(edges))
@time core_numbers = fast_kcore(nodes, edges)

n = length(nodes)
degrees = zeros(Int64, n)
for edge in eachrow(edges)
    degrees[edge[1]] += 1
    degrees[edge[2]] += 1
end

coreness = zeros(Int64, n)
for key in keys(core_numbers)
    coreness[key] = core_numbers[key]
end

density = zeros(Int64, maximum(degrees), maximum(coreness))
for i in 1:n
    density_matrix[degrees[i], coreness[i]] += 1
end

a = Vector()
b = Vector()
c = Vector()

for element in findall(!iszero, density)
    push!(a, element[1])
    push!(b, element[2])
    push!(c, density[element])
end

figure()
title("Google scholar dataset")
scatter(x=a, y=b, c=c, s = c/maximum(c) * 2000)
xscale("log")
xlabel("Degree")
ylabel("Coreness")
colorbar()
plot(1:14, 1:14, c="black")

ext: jl

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
#############################################################################
#This file contains the implementation of the the functions that create the
#multiresolution ordering and sparsity pattern from an input consisting of
#a set of points
#############################################################################

#############################################################################
#Implementation of the Heap
#############################################################################

#HeapNode containing the global ID and the distance squared to the points added
mutable struct HeapNode
  Id::Int64
  dist2::Float64
end

#Heap containing its size (for convenience), an array of HeapNodes,
#and and array mapping the Id's of the nodes to their position in the heap.
struct Heap
  N::Int64
  nodes::Array{HeapNode,1}
  lookup::Array{Int64,1}
end

#Funktion swapping the nodes a and b, as identified by their respective heap ID,
#in the heap. Presently implemented by swapping the two values of the nodes.
function _swap!( h::Heap, aheapInd::Int64, bheapInd::Int64 )
  #Assinging the corresponding new values to the lookup table
  h.lookup[ h.nodes[aheapInd].Id ] = bheapInd
  h.lookup[ h.nodes[bheapInd].Id ] = aheapInd
  #Introducing temporary variables for the swapping.
  tempDist2::Float64 = h.nodes[ aheapInd ].dist2
  tempId::Int64 = h.nodes[ aheapInd ].Id
  #Swapping the entries of the corresponding node structs
  h.nodes[ aheapInd ].dist2 = h.nodes[ bheapInd ].dist2
  h.nodes[ aheapInd ].Id = h.nodes[ bheapInd ].Id
  h.nodes[ bheapInd ].dist2 = tempDist2
  h.nodes[ bheapInd ].Id = tempId
end

#Compares the element with heap Index hInd to its two children and if it is smaller
#than one of them, swaps it with the bigger one.
#The output is the the new hindex of the node, if a swap occured, and -1, otherwise.
function _moveDown!( h::Heap, hInd::Int64 )
  #Saving the squared distance of the node that is being swapped down.
  dist2::Float64 = h.nodes[ hInd ].dist2
  #If the both children exist:
  if 2 * hInd + 1 <= h.N
    #If the left child is larger:
    if h.nodes[ 2 * hInd ].dist2 >= h.nodes[ 2 * hInd + 1 ].dist2
      #Check if the Child is larger than the parent:
      if h.nodes[2*hInd ].dist2 >= dist2
        _swap!( h, hInd, 2 * hInd )
        return 2 * hInd
      else
        #No swap occuring:
        return -1
      end
    #If the left child is larger:
    else
      #Check if the Child is larger than the parent:
      if h.nodes[ 2*hInd + 1 ].dist2 >= dist2
        _swap!( h, hInd, 2 * hInd + 1 )
        return 2 * hInd + 1
      else
        #No swap occuring:
        return -1
      end
    end
  #If only one child exists:
  elseif 2 * hInd <= h.N
    #Check if the Child is larger than the parent:
    if h.nodes[ 2 * hInd ].dist2 > dist2
      _swap!( h, hInd, 2 * hInd )
      return 2 * hInd
    end
  end
  #No swap occuring:
  return -1
end

function topId( h::Heap )
  return h.nodes[1].Id
end

function topDist2( h::Heap )
  return h.nodes[1].dist2
end

#Updates (decreases) an element of the heap and restores the heap property.
function update!( h::Heap, Id::Int64, newDist2::Float64 )
  tempInd::Int64 = h.lookup[ Id ]
  h.nodes[ tempInd ].dist2 = newDist2
  while ( tempInd <= h.N ) && ( tempInd >= 0 )
    tempInd = _moveDown!( h, tempInd )
  end
end

#############################################################################
#Introducing the "daycare", which keeps track of the descendants of every node
#############################################################################

mutable struct Daycare
  NParents::Int64
  NChildren::Int64
  NBuffer::Int64

  #This array gives for any parent node the corresponding ID
  IdParents::Array{Int64, 1}
  firstChildren::Array{Int64, 1}
  children::Array{Int64, 1}
  #This array has $N$ entries, that given node Id, the corresponding entry in the list
  #of parents
  lookup::Array{Int64, 1}
end

function newparent( dc::Daycare, IdParent::Int64 )
  dc.NParents = dc.NParents + 1
  dc.IdParents[ dc.NParents ] = IdParent
  dc.firstChildren[ dc.NParents ] = dc.NChildren + 1
  dc.firstChildren[ dc.NParents + 1 ] = dc.NChildren + 1
  dc.lookup[ IdParent ] = dc.NParents
end

function newson( dc::Daycare, IdSon::Int64 )
  if dc.NChildren >= dc.NBuffer
    dc.NBuffer = 2 * dc.NBuffer
    resize!( dc.children, dc.NBuffer )
  end
  dc.NChildren = dc.NChildren + 1
  dc.firstChildren[ dc.NParents + 1 ] = dc.firstChildren[ dc.NParents + 1 ] + 1
  dc.children[ dc.NChildren ] = IdSon
end

#############################################################################
#Implementation of the ordering algorithm
#############################################################################

#############################################################################
#Version without levels
#############################################################################
function _determineChildren!( h::Heap,
                             dc::Daycare,
                             parents::Array{Int64, 1},
                             x::Array{Float64, 2},
                             ID::Int64 )
#Function to determine the children of a node in the ordering
#Inputs:
# h:
#   Heap to keep track of the elements and their distance to the points added
#   the ordering
# dc:
#   The "daycare" keeping track of the nodes and their children
# parents:
#   An array with N elements that in its i-th position contains the Id of the point
#   with Id i.
# x:
#   A d times N array containing the positions of the data points
# ID:
#   The ID of the point, the children of which are being determined

  #caching the Id of the father node
  father::Int64 = parents[ ID ]
  #caching the point
  y::Array{Float64, 1} = x[ :, ID ]
  #adding the a new node to the dc
  newparent( dc, ID )
  #storing the ambient dimension
  d::Int64 = size( x, 1 )

  #caching the distance of the newest elemnt
  pivotDist2::Float64 = topDist2( h )

  #We need to add the children to the new parent
  for index = dc.firstChildren[ dc.lookup[ father ] ]  : ( dc.firstChildren[ dc.lookup[ father ] + 1 ] - 1 )
    j::Int64 = dc.children[ index ]
    dist2::Float64 = dist2eval( x, j , y, d )

    if (dist2 <= pivotDist2)&&( h.nodes[ h.lookup[ j ] ].dist2 >= 0. )
      if dist2 < h.nodes[ h.lookup[ j ] ].dist2
        update!( h, j, dist2 )
      end
      newson( dc, j )
      if sqrt.( h.nodes[ h.lookup[ j ] ].dist2 ) .+ sqrt.( dist2 ) <= sqrt.( pivotDist2 )
        parents[ j ] = ID
      end
    end
  end

  #updating the newest element to distance -1, to indicate that it has been
  #included in the hierarchy already (important in case multiple measurements
  #have the same position
  update!( h, ID, -1. )

end

#function sortPoints( x::Array{Float64, 2}, boundDist::Array{Float64, 1} )
##Function to sort the points in a multicscale way
##Input:
## x:
##   A d times N array containing the measurement points
#  Nbuffer0::Int64 = 10
#  N::Int64 = size( x, 2 )
#  parents::Array{Int64, 1} = Array{Int64, 1}(N)
#  #Create heap:
#  h::Heap = Heap( N, Array{HeapNode, 1}( N ), 1 : 1 : N )
#  for i = 1 : N
#    h.nodes[ i ] = HeapNode( i, 1000000000000 )
#  end
#  #Create daycare
#  dc::Daycare = Daycare( 0, 0, Nbuffer0, Array{Int64, 1}(N), Array{Int64, 1}(N+1),  Array{Int64, 1}(Nbuffer0), Array{Int64, 1}(N))
#  #TODO remove ( only debugging )
#  dc.IdParents .= 0
#  dc.firstChildren .= 0
#  dc.children .= 0
#  dc.lookup .= 0
#  parents .= 0
#  #distances will contain as an i-th entry the distance upon elimination of the i-th
#  #degree of freedom
#  distances::Array{Float64} = zeros( N )
#
#
#  #Setting the first element:
#  #rootID::Int64 = randperm( N )[ 1 ]
#  rootID::Int64 = 1;
#  rootPoint::Array{Float64, 1} = x[ :, rootID ]
#  parents .= rootID
#  newparent( dc, rootID )
#  for i = 1 : N
#    newson( dc, i )
#    distances[1] = max( distances[1], sqrt.( dist2eval( x, i, rootPoint, size( x, 1 ) ) ) )
#    update!( h, i, dist2eval( x, i, rootPoint, size( x, 1 ) ) )
#  end
#
#  #updating the elements
#  for iter = 2 : N
#    pivotId::Int64 = topId( h )
#    distances[ iter ] = sqrt( topDist2( h ) )
#    _determineChildren!( h, dc, parents, x, pivotId )
#  end
#
#  #Output "P, revP, levels"
#  #Adding a vector with the distances, at which the
#
#  return dc.IdParents, dc.lookup, distances
#end

#Evaluates the distance of the point y in d-dimensional space to the jth point
#of the point cloud x
function dist2eval( x::Array{Float64, 2}, j::Int64, y::Array{Float64, 1}, d::Int64 )
  retVal::Float64 = 0.0;
  for i::Int64 = 1 : d
    retVal += ( x[ i, j ] - y[ i ] )^2
  end
  return retVal
end

#############################################################################
#Implementation of the sparsity pattern
#############################################################################

function _determineChildren!( revP::Array{Int64, 1}, distances::Array{Float64, 1}, dc::Daycare, parents::Array{Int64, 1}, x::Array{Float64, 2}, ID::Int64, rho::Float64 )
  #caches the father node
  father::Int64 = parents[ ID ]
  #caches the pivot point
  y::Array{Float64, 1} = x[ :, ID ]
  #Adds the new parent (row) to be added to the sparsity pattern
  newparent( dc, ID )
  #caches the ambient dimension
  d::Int64 = size( x, 1 )
  #caches the distance coresponding to the pivot
  pivotDist::Float64 = distances[ revP[ ID ] ]
  for index = dc.firstChildren[ dc.lookup[ father ] ]  : ( dc.firstChildren[ dc.lookup[ father ] + 1 ] - 1 )
    j::Int64 = dc.children[ index ]
    dist2::Float64 = dist2eval( x, j , y, d )

    if ( revP[ ID ] <= revP[ j ] ) && ( dist2 <= rho^2 * ( pivotDist.^2 ) )
      newson( dc, j )
      if sqrt.( dist2 ) + rho * distances[ revP[ j ] ] <= rho * pivotDist
        parents[ j ] = ID
      end
    end
  end
end

function sortPoints( x::Array{Float64, 2} )
#Function to sort the points in a multicscale way
#Input:
# x:
#   A d times N array containing the measurement points
  Nbuffer0::Int64 = 10
  N::Int64 = size( x, 2 )
  parents::Array{Int64, 1} = Array{Int64, 1}(N)
  #Create heap:
  h::Heap = Heap( N, Array{HeapNode, 1}( N ), 1 : 1 : N )
  for i = 1 : N
    h.nodes[ i ] = HeapNode( i, 1000000000000 )
  end
  #Create daycare
  dc::Daycare = Daycare( 0, 0, Nbuffer0, Array{Int64, 1}(N), Array{Int64, 1}(N+1),  Array{Int64, 1}(Nbuffer0), Array{Int64, 1}(N))
  #TODO remove ( only debugging )
  dc.IdParents .= 0
  dc.firstChildren .= 0
  dc.children .= 0
  dc.lookup .= 0
  parents .= 0
  #distances will contain as an i-th entry the distance upon elimnation of the i-th
  #degree of freedom
  distances::Array{Float64} = zeros( N )

  #Setting the first element:
  #rootID::Int64 = randperm( N )[ 1 ]
  rootID::Int64 = 1;
  rootPoint::Array{Float64, 1} = x[ :, rootID ]
  parents .= rootID
  newparent( dc, rootID )
  for i = 1 : N
    newson( dc, i )
    distances[1] = max( distances[1], sqrt.( dist2eval( x, i, rootPoint, size( x, 1 ) ) ) )
    update!( h, i, dist2eval( x, i, rootPoint, size( x, 1 ) ) )
  end
  update!( h, rootID, -1. )

  #updating the elements
  for iter = 2 : N
    pivotId::Int64 = topId( h )
    distances[ iter ] = sqrt( topDist2( h ) )
    _determineChildren!( h, dc, parents, x, pivotId )
  end

  #Output "P, revP, levels"
  #Adding a vector with the distances, at which the
  return dc.IdParents, dc.lookup, distances
end

#############################################################################
#Implementation of the sparsity pattern
#############################################################################

function sparsityPattern( x::Array{Float64, 2}, P::Array{Int64, 1}, revP::Array{Int64, 1}, distances::Array{Float64, 1}, rho::Float64 )
  d::Int64 = size( x, 1 )
  N::Int64 = size( x, 2 )
  Nbuffer0::Int64 = 10

  parents::Array{Int64, 1} = Array{Int64, 1}(N)
  dc::Daycare = Daycare( 0, 0, Nbuffer0, Array{Int64, 1}(N), Array{Int64, 1}(N+1),  Array{Int64, 1}(Nbuffer0), Array{Int64, 1}(N))

  dc.IdParents .= 0
  dc.firstChildren .= 0
  dc.children .= 0
  dc.lookup .= 0
  parents .= 0

  parents .= P[ 1 ]
  newparent( dc, P[ 1 ] )
  for i = 1 : N
    newson( dc, i )
  end

  for k = 2 : N
    _determineChildren!( revP, distances, dc, parents, x, P[ k ], rho)
  end

  #Sort the jind according to their ordering:
  for i = 1 : N
    dc.children[ dc.firstChildren[ i ] : ( dc.firstChildren[ i + 1 ] - 1 ) ] =
      P[ sort!( revP[ dc.children[ dc.firstChildren[ i ] : ( dc.firstChildren[ i + 1 ] - 1 ) ] ] ) ]
  end

  return dc.firstChildren, revP[ dc.children[ 1 : dc.NChildren ] ]
end

max_stars_repo_name: f-t-s/Spelliptic

max_issues_repo_name: f-t-s/Spelliptic

max_forks_repo_name: f-t-s/Spelliptic


text:
using Test
using HTTP

function generate_test_body()
    Vector{UInt8}("----------------------------918073721150061572809433\r\nContent-Disposition: form-data; name=\"namevalue\"; filename=\"multipart.txt\"\r\nContent-Type: text/plain\r\n\r\nnot much to say\n\r\n----------------------------918073721150061572809433\r\nContent-Disposition: form-data; name=\"key1\"\r\n\r\n1\r\n----------------------------918073721150061572809433\r\nContent-Disposition: form-data; name=\"key2\"\r\n\r\nkey the second\r\n----------------------------918073721150061572809433\r\nContent-Disposition: form-data; name=\"namevalue2\"; filename=\"multipart-leading-newline.txt\"\r\nContent-Type: text/plain\r\n\r\n\nfile with leading newline\n\r\n----------------------------918073721150061572809433--\r\n")
end

function generate_test_request()
    headers = [
        "User-Agent" => "PostmanRuntime/7.15.2",
        "Accept" => "*/*",
        "Cache-Control" => "no-cache",
        "Postman-Token" => "288c2481-1837-4ba9-add3-f23d380fa440",
        "Host" => "localhost:8888",
        "Accept-Encoding" => "gzip, deflate",
        "Accept-Encoding" => "gzip, deflate",
        "Content-Type" => "multipart/form-data; boundary=--------------------------918073721150061572809433",
        "Content-Length" => "657",
        "Connection" => "keep-alive",
    ]

    HTTP.Request("POST", "/", headers, generate_test_body())
end

function generate_non_multi_test_request()
    headers = [
        "User-Agent" => "PostmanRuntime/7.15.2",
        "Accept" => "*/*",
        "Cache-Control" => "no-cache",
        "Postman-Token" => "288c2481-1837-4ba9-add3-f23d380fa440",
        "Host" => "localhost:8888",
        "Accept-Encoding" => "gzip, deflate",
        "Accept-Encoding" => "gzip, deflate",
        "Content-Length" => "657",
        "Connection" => "keep-alive",
    ]

    HTTP.Request("POST", "/", headers, Vector{UInt8}())
end

@testset "parse multipart form-data" begin
    @testset "find_multipart_boundary" begin
        request = generate_test_request()

        # NOTE: this is the start of a "boundary delimiter line" and has two leading
        # '-' characters prepended to the boundary delimiter from Content-Type header
        delimiter = Vector{UInt8}("----------------------------918073721150061572809433")
        body = generate_test_body()
        # length of the delimiter, CRLF, and -1 for the end index to be the LF character
        endIndexOffset = length(delimiter) + 2 - 1

        (isTerminatingDelimiter, startIndex, endIndex) = HTTP.find_multipart_boundary(body, delimiter)
        @test !isTerminatingDelimiter
        @test 1 == startIndex
        @test (startIndex + endIndexOffset) == endIndex

        # the remaining "boundary delimiter lines" will have a CRLF preceding them
        endIndexOffset += 2

        (isTerminatingDelimiter, startIndex, endIndex) = HTTP.find_multipart_boundary(body, delimiter, start = startIndex + 1)
        @test !isTerminatingDelimiter
        @test 175 == startIndex
        @test (startIndex + endIndexOffset) == endIndex

        (isTerminatingDelimiter, startIndex, endIndex) = HTTP.find_multipart_boundary(body, delimiter, start = startIndex + 3)
        @test !isTerminatingDelimiter
        @test 279 == startIndex
        @test (startIndex + endIndexOffset) == endIndex

        (isTerminatingDelimiter, startIndex, endIndex) = HTTP.find_multipart_boundary(body, delimiter, start = startIndex + 3)
        @test !isTerminatingDelimiter
        @test 396 == startIndex
        @test (startIndex + endIndexOffset) == endIndex

        (isTerminatingDelimiter, startIndex, endIndex) = HTTP.find_multipart_boundary(body, delimiter, start = startIndex + 3)
        @test isTerminatingDelimiter
        @test 600 == startIndex
        # +2 because of the two additional '--' characters
        @test (startIndex + endIndexOffset + 2) == endIndex
    end

    @testset "parse_multipart_form" begin
        @test HTTP.parse_multipart_form(generate_non_multi_test_request()) === nothing

        multiparts = HTTP.parse_multipart_form(generate_test_request())
        @test 4 == length(multiparts)

        @test "multipart.txt" === multiparts[1].filename
        @test "namevalue" === multiparts[1].name
        @test "text/plain" === multiparts[1].contenttype
        @test "not much to say\n" === String(read(multiparts[1].data))

        @test multiparts[2].filename === nothing
        @test "key1" === multiparts[2].name
        @test "text/plain" === multiparts[2].contenttype
        @test "1" === String(read(multiparts[2].data))

        @test multiparts[3].filename === nothing
        @test "key2" === multiparts[3].name
        @test "text/plain" === multiparts[3].contenttype
        @test "key the second" === String(read(multiparts[3].data))

        @test "multipart-leading-newline.txt" === multiparts[4].filename
        @test "namevalue2" === multiparts[4].name
        @test "text/plain" === multiparts[4].contenttype
        @test "\nfile with leading newline\n" === String(read(multiparts[4].data))
    end
end

@testset "content_disposition_extract($(v[1])" for v in (
        ("; filename=abc.txt ; name = xyz", "xyz", "abc.txt"),
        ("; name=abc ; filename = xyz", "abc", "xyz"),
        ("""; mno;filename="abc";name=xyz""", "xyz", "abc"),
        (""";filename="abc";mno;name=xyz""", "xyz", "abc"),
        (""";filename=   "abc"   ;mno;name=xyz""", "xyz", "abc"),
        ("; filename=abc.txt ; name = xyz ;", "xyz", "abc.txt"),
        ("; filename=abc.txt ; name = xyz;", "xyz", "abc.txt"),
        ("; filename=abc.txt ; name = xyz ; mno", "xyz", "abc.txt"),
        ("; filename=abc.txt ; name = xyz ; mno ;", "xyz", "abc.txt"),
        (";name=\"ab\\\"cdef\"","ab\\\"cdef", nothing),
        (";filename=abc\\;xyz", nothing, "abc\\;xyz"),
        (";filename=\\\"abc;name=xyz", "xyz", "\\\"abc"),
        (";name=xyz;filename=;mno", "xyz", nothing),
        (";name=\"xy;z\";filename=;mno", "xy;z", nothing),
        (";name=\"x=z\";filename=bbb", "x=z", "bbb")
        )
    name = nothing
    filename = nothing

    for (pair, key, value) in HTTP.content_disposition_extract(v[1])
        if pair && key == "name"
            name = value
        elseif pair && key == "filename"
            filename = value
        end
    end

    @test name == v[2]
    @test filename == v[3]
end


text:
"""
Datatype for holding usearch and mafft file paths
"""
type Paths
    usearch::String
    mafft::String
end

"""
Default file paths.
"""
const PATHS = Paths("usearch", "mafft")


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

# test default packages

cd(joinpath(JULIA_HOME,"..","share","julia","site","v$(VERSION.major).$(VERSION.minor)")) do
    pkgs = readdir()
    if startswith(string(Sys.ARCH), "arm")
        # Remove profile from default tests on ARM since it currently segfaults
        # Allow explicitly adding it for testing
        warn("Skipping Profile tests")
        filter!(x -> (x != "Profile"), pkgs)
    end
    Pkg.Entry.test(convert(Vector{AbstractString}, pkgs))
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
# generate LLVM wrappers

using LLVM_full_jll

using Clang

function wrap()
    includedir = LLVM_full_jll.llvm_config() do config
        readchomp(`$config --includedir`)
    end
    cppflags = LLVM_full_jll.llvm_config() do config
        split(readchomp(`$config --cppflags`))
    end

    # Set-up arguments to clang
    clang_includes  = map(x->x[3:end], filter( x->startswith(x,"-I"), cppflags))
    clang_extraargs =                  filter(x->!startswith(x,"-I"), cppflags)

    # Recursively discover LLVM C API headers (files ending in .h)
    header_dirs = String[joinpath(includedir, "llvm-c")]
    header_files = String[]
    while !isempty(header_dirs)
        parent = pop!(header_dirs)
        children = readdir(parent)
        for child in children
            path = joinpath(parent, child)
            if isdir(path)
                push!(header_dirs, path)
            elseif isfile(path) && endswith(path, ".h")
                push!(header_files, path)
            end
        end
    end

    context = init(;
                    headers = header_files,
                    output_file = "libLLVM.jl",
                    common_file = "libLLVM_common.jl",
                    clang_includes = convert(Vector{String}, clang_includes),
                    clang_args = convert(Vector{String}, clang_extraargs),
                    header_library = x->"libllvm",
                    header_wrapped = (top,cursor)->occursin("include/llvm", cursor)
                  )

    run(context)
end

function main()
    cd(joinpath(dirname(@__DIR__), "lib")) do
        wrap()
    end
end

isinteractive() || main()

# Manual clean-up:
# - remove build-host details (HAVE_INTTYPES_H, LLVM_DEFAULT_TARGET_TRIPLE etc) in libLLVM_common.jl
# - remove LLVMInitializeAll and LLVMInitializeNative wrappers (these are macros)
# - remove "# Skipping ..." comments by Clang.jl
# - replace `ccall\(\((:.+), libllvm\), (.*)\)` with `@runtime_ccall($1, $2)`
# - replace `const (LLVMOpaque.*) = Cvoid` with `struct $1 end`
# - use `gawk -i inplace '/^[[:blank:]]*$/ { print; next; }; {cur = seen[$0]; if(!seen[$0]++ || (/^end$/ && !prev) || /^.*Clang.*$/) print $0; prev=cur}' libLLVM_h.jl` to remove duplicates
# - use `cat -s` to remove duplicate empty lines


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

using Base.Test

@test @__LINE__ == 5

include("test_sourcepath.jl")
thefname = "the fname!//\\&\1*"
include_string_test_func = include_string("include_string_test() = @__FILE__", thefname)
@test include_string_test_func() == Base.source_path()
@test include_string("Base.source_path()", thefname) == Base.source_path()
@test basename(@__FILE__) == "loading.jl"
@test isabspath(@__FILE__)

@test isdir(@__DIR__)
@test @__DIR__() == dirname(@__FILE__)

# Issue #5789 and PR #13542:
mktempdir() do dir
    cd(dir) do
        let true_filename = "cAsEtEsT.jl", lowered_filename="casetest.jl"
            touch(true_filename)
            @test Base.isfile_casesensitive(true_filename)
            @test !Base.isfile_casesensitive(lowered_filename)

            # check that case-sensitivity only applies to basename of a path:
            if isfile(lowered_filename) # case-insensitive filesystem
                mkdir("cAsEtEsT")
                touch(joinpath("cAsEtEsT", true_filename))
                @test Base.isfile_casesensitive(joinpath("casetest", true_filename))
                @test !Base.isfile_casesensitive(joinpath("casetest", lowered_filename))
            end
        end

        # Test Unicode normalization; pertinent for OS X
        let nfc_name = "\U00F4.jl"
            touch(nfc_name)
            @test Base.isfile_casesensitive(nfc_name)
        end
    end
end

SAVED_LOAD_PATH = copy(LOAD_PATH)
empty!(LOAD_PATH)
dir = abspath(@__DIR__)
push!(LOAD_PATH, dir)

@test Base.find_in_path("test_sourcepath") == joinpath(dir, "test_sourcepath.jl")
@test Base.find_in_path(GenericString("test_sourcepath")) == joinpath(dir, "test_sourcepath.jl")
LOAD_PATH[end] = GenericString(LOAD_PATH[end])
@test Base.find_in_path("test_sourcepath") == joinpath(dir, "test_sourcepath.jl")

struct CustomLoader
    path::String
end
push!(LOAD_PATH, CustomLoader("abc"))
let name = randstring(20)
    @test_throws ArgumentError Base.find_in_path(name, nothing)
    Base.load_hook(prefix::CustomLoader, name::String, found) = joinpath(prefix.path, name)
    @test Base.find_in_path(name, nothing) == joinpath("abc", name)
end
@test Base.find_in_path("test_sourcepath", nothing) == joinpath("abc", "test_sourcepath")
Base.load_hook(prefix::CustomLoader, name::String, found::String) = found
@test Base.find_in_path("test_sourcepath", nothing) == joinpath(dir, "test_sourcepath.jl")

empty!(LOAD_PATH)
append!(LOAD_PATH, SAVED_LOAD_PATH)
@test LOAD_PATH == SAVED_LOAD_PATH

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
# enable precompilation
__precompile__(true)

module Robotics

import Base: show, *, size, +

#=------------------------------------------------------------------------------
                           Includes
------------------------------------------------------------------------------=#
include("transformations.jl")
include("robots_types.jl")
include("utility.jl")
include("graphics.jl")
include("types.jl")
#=------------------------------------------------------------------------------
                              Exports
------------------------------------------------------------------------------=#
export   #types
         Point2D,
         #methods
         ξ,
         rot2,
         trot2,
         rotx,
         roty,
         rotz,
         show,
         euler2rot,
         rot2euler,
         angle_vector2rot,
         rot2trans,
         trotx,
         skew,
         skewa,
         vex,
         # graphics
         plot_frame,
         # utility
         generate_random_pose
end
#
# Robotics.jl ends here

max_stars_repo_name: elsuizo/Robot

max_issues_repo_name: elsuizo/Robot

max_forks_repo_name: elsuizo/Robot


text:
# deep copying between cpu/gpu
#
# I duplicated the code from base/deepcopy.jl to cpucopy and gpucopy:
#
# deepcopy(x)
#   deepcopy_internal(x,s)
#     _deepcopy_t(x,T,s)
#   deepcopy_internal(x::Array,s)
#     _deepcopy_array_t(x,T,s)
#
# cpucopy and gpucopy differ from deepcopy only in how Arrays,
# CudaArrays and other types based on them are handled.  New types can
# override deepcopy_internal, cpucopy_internal and gpucopy_internal.

using CUDArt
using Base: arrayset
VERSION < v"0.4-" && (nfields(a)=length(names(a)))

cpucopy(x) = cpucopy_internal(x, ObjectIdDict())

cpucopy_internal(x::Union{Symbol,LambdaStaticData,TopNode,GlobalRef,
                 DataType,Union,Task},
                 stackdict::ObjectIdDict) = x
cpucopy_internal(x::Tuple, stackdict::ObjectIdDict) =
    ntuple(i->cpucopy_internal(x[i], stackdict), length(x))
cpucopy_internal(x::Module, stackdict::ObjectIdDict) = error("cpucopy of Modules not supported")

function cpucopy_internal(x::Function, stackdict::ObjectIdDict)
    if isa(x.env, Union{MethodTable, Symbol}) || x.env === ()
        return x
    end
    invoke(cpucopy_internal, Tuple{Any, ObjectIdDict}, x, stackdict)
end

function cpucopy_internal(x, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _cpucopy_t(x, typeof(x), stackdict)
end

function _cpucopy_t(x, T::DataType, stackdict::ObjectIdDict)
    nf = nfields(T)
    (isbits(T) || nf == 0) && return x
    if T.mutable
        y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
        stackdict[x] = y
        for i in 1:nf
            if isdefined(x,i)
                y.(i) = cpucopy_internal(x.(i), stackdict)
            end
        end
    else
        fields = Any[cpucopy_internal(x.(i), stackdict) for i in 1:nf]
        y = ccall(:jl_new_structv, Any, (Any, Ptr{Void}, UInt32),
                  T, pointer(fields), length(fields))
    end
    return y::T
end

function cpucopy_internal(x::Array, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _cpucopy_array_t(x, eltype(x), stackdict)
end

function _cpucopy_array_t(x, T, stackdict::ObjectIdDict)
    if isbits(T)
        return (stackdict[x]=copy(x))
    end
    dest = similar(x)
    stackdict[x] = dest
    for i=1:length(x)
        if isdefined(x,i)
            arrayset(dest, cpucopy_internal(x[i], stackdict), i)
        end
    end
    return dest
end

gpucopy(x) = gpucopy_internal(x, ObjectIdDict())

gpucopy_internal(x::Union{Symbol,LambdaStaticData,TopNode, # GlobalRef,
                            DataType,Union,Task},
                   stackdict::ObjectIdDict) = x
gpucopy_internal(x::Tuple, stackdict::ObjectIdDict) =
    ntuple(i->gpucopy_internal(x[i], stackdict), length(x))
gpucopy_internal(x::Module, stackdict::ObjectIdDict) = error("gpucopy of Modules not supported")

function gpucopy_internal(x::Function, stackdict::ObjectIdDict)
    if isa(x.env, Union{MethodTable, Symbol}) || x.env === ()
        return x
    end
    invoke(gpucopy_internal, Tuple{Any, ObjectIdDict}, x, stackdict)
end

function gpucopy_internal(x, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _gpucopy_t(x, typeof(x), stackdict)
end

function _gpucopy_t(x, T::DataType, stackdict::ObjectIdDict)
    nf = nfields(T)
    (isbits(T) || nf == 0) && return x
    if T.mutable
        y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
        stackdict[x] = y
        for i in 1:nf
            if isdefined(x,i)
                y.(i) = gpucopy_internal(x.(i), stackdict)
            end
        end
    else
        fields = Any[gpucopy_internal(x.(i), stackdict) for i in 1:nf]
        y = ccall(:jl_new_structv, Any, (Any, Ptr{Void}, UInt32),
                  T, pointer(fields), length(fields))
    end
    return y::T
end

function gpucopy_internal(x::Array, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _gpucopy_array_t(x, eltype(x), stackdict)
end

function _gpucopy_array_t(x, T, stackdict::ObjectIdDict)
    if isbits(T)
        return (stackdict[x]=copy(x))
    end
    dest = similar(x)
    stackdict[x] = dest
    for i=1:length(x)
        if isdefined(x,i)
            arrayset(dest, gpucopy_internal(x[i], stackdict), i)
        end
    end
    return dest
end


text:
module Evaluator
export evaluate

using Combinatorics
using Random
using ..PokerCard
using ..PokerDeck
using ..Lookup

unsuited_lookup = Lookup.unsuited_lookup
flush_lookup = Lookup.flush_lookup

# check we've covered all of the unique hands
@assert length(unsuited_lookup) + length(flush_lookup) == 7462

function five_card(cards::Array)
    # first we check for a flush
    cond = cards[1] & cards[2] & cards[3] & cards[4] & cards[5] & 61440 # 0xF000
    if sum(cond) > 0
        handOR = (cards[1] | cards[2] | cards[3] | cards[4] | cards[5]) >> 16
        prime = PokerCard.prime_product_from_rankbits(handOR)
        return Lookup.flush_lookup[prime]
    else
        prime = PokerCard.prime_product_from_hand(cards)
        return Lookup.unsuited_lookup[prime]
    end
end

function six_and_seven_card(cards::Array)
    """
    Performs five_card_eval() on all (6 choose 5) = 6 subsets
    of 5 cards in the set of 6 to determine the best ranking,
    and returns this ranking. Same for 7 cards.
    """
    minimum_score = Lookup.MAX_HIGH_CARD
    all5cardcombobs = collect(Combinatorics.combinations(cards, 5))
    for combo in all5cardcombobs
        score = five_card(combo)
        if score < minimum_score
            minimum_score = score
        end
    end
    return minimum_score
end

function evaluate(cards::Array)
    l = length(cards)
    @assert (4 < l < 8)
    if l > 5
        return six_and_seven_card(cards)
    else
        return five_card(cards)
    end
end

function get_rank_class(hr::Int)
    """
    Returns the class of hand given the hand hand_rank
    returned from evaluate.
    """
    if (hr >= 0) & (hr <= Lookup.MAX_STRAIGHT_FLUSH)
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_STRAIGHT_FLUSH]
    elseif hr <= Lookup.MAX_FOUR_OF_A_KIND
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_FOUR_OF_A_KIND]
    elseif hr <= Lookup.MAX_FULL_HOUSE
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_FULL_HOUSE]
    elseif hr <= Lookup.MAX_FLUSH
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_FLUSH]
    elseif hr <= Lookup.MAX_STRAIGHT
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_STRAIGHT]
    elseif hr <= Lookup.MAX_THREE_OF_A_KIND
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_THREE_OF_A_KIND]
    elseif hr <= Lookup.MAX_TWO_PAIR
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_TWO_PAIR]
    elseif hr <= Lookup.MAX_PAIR
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_PAIR]
    elseif hr <= Lookup.MAX_HIGH_CARD
        return Lookup.MAX_TO_RANK_CLASS[Lookup.MAX_HIGH_CARD]
    else
        throw(DomainError())
    end

end

function class_to_string(class_int::Int)
    """
    Converts the integer class hand score into a human-readable string.
    """
    return Lookup.RANK_CLASS_TO_STRING[class_int]
end

function get_five_card_rank_percentage(hand_rank::Int)
    """
    Scales the hand rank score to the [0.0, 1.0] range.
    """
    return float(hand_rank) / float(Lookup.MAX_HIGH_CARD)
end

end #Evaluator

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
using VisualDL
train_logger = VisualDLLogger("tmp", 1, "train")
test_logger = as_mode(train_logger, "test")

for i in 1:100
    with_logger(train_logger) do
       @log_histogram h0=(i, randn(100))
    end
end


text:
function nz_stream_foodweb()
  data_path = joinpath(@__DIR__, "../..", "data", "nz_stream")
  files = readdir(data_path)
  documents = filter(f -> endswith(f, ".csv"), files)
  Ns = UnipartiteNetwork{Bool,String}[]
  for doc in documents
    content = readdlm(joinpath(@__DIR__, "../..", "data", "nz_stream", doc), ',')
    species_names = convert(Array{String}, content[:,1][2:end])
    interaction_matrix = map(Int64, content[2:end,2:end]).>0
    push!(Ns, UnipartiteNetwork(interaction_matrix', species_names))
  end
  return Ns
end

function web_of_life(name)
  file_pathname = name * ".csv"
  data_path = joinpath(@__DIR__, "../..", "data", "weboflife")
  files = readdir(data_path)
  @assert file_pathname in files
  content = readdlm(joinpath(@__DIR__, "../..", "data", "weboflife", file_pathname), ',')
  bottom_species_names = convert(Array{String}, vec(content[:,1][2:end]))
  top_species_names = convert(Array{String}, vec(content[1,:][2:end]))
  int_mat = content[2:end,2:end]
  repl = findall(int_mat .== "#")
  int_mat[repl] .= 0.0
  interaction_matrix = map(Int64, round.(permutedims(int_mat)))
  ntype = BipartiteQuantitativeNetwork
  if maximum(interaction_matrix) == 1
    interaction_matrix = interaction_matrix .> 0
    ntype = BipartiteNetwork
  end
  this_net = ntype(interaction_matrix, top_species_names, bottom_species_names)
  forbidden_names = ["Numbers of flowers", "Num. of hosts sampled", "#", "Frequency of occurrences", "Number of flowers", "Number of droppings analysed"]
  sp_1 = filter(x -> !(x ∈ forbidden_names), species(this_net; dims=1))
  sp_2 = filter(x -> !(x ∈ forbidden_names), species(this_net; dims=2))
  return this_net[sp_1, sp_2]
end

function web_of_life()
  wol_ref = joinpath(@__DIR__, "../..", "data", "weboflife", "references.csv")
  wol_infos = readdlm(wol_ref, ',')
  names = [Symbol(replace(n, " " => "_")) for n in wol_infos[1,:]]
  infos = [NamedTuple{tuple(names...)}(tuple(wol_infos[i,:]...)) for i in 2:size(wol_infos,1)]
  return infos
end

function pajek(file)
    pajek_lines = readlines(file)
    filter!(line -> !startswith(line, "% "), pajek_lines)
    FLAG_vertices = false
    FLAG_arcs = false
    FLAG_read = false
    nodes = Dict{Any,Any}()
    arcs = Any[]
    for line in pajek_lines
        if startswith(line, "*network ")
            FLAG_read = true
        end
        if FLAG_read
            # What are we reading?
            if startswith(line, "*vertices ")
                FLAG_vertices = true
                FLAG_arcs = false
                continue
            end
            if startswith(line, "*arcs ")
                FLAG_vertices = false
                FLAG_arcs = true
                continue
            end
            # Reading the line
            if FLAG_vertices
                nodematch = match(r"\s+(\d+)\s+\"(.+)\"", line)
                if typeof(nodematch) <: Nothing
                    FLAG_vertices = false
                    FLAG_arcs = true
                    continue
                else
                    nodeid, nodename = nodematch.captures
                    nodes[nodeid] = nodename
                end
            end
            if FLAG_arcs
                arcmatch = match(r"\s+(\d+)\s+(\d+)\s+(.+)\s?", line)
                if typeof(arcmatch) <: Nothing
                    continue
                else
                    nodefrom, nodeto, nodestrength = arcmatch.captures
                    push!(arcs, (nodefrom, nodeto, nodestrength))
                end
            end
        end
    end
    S = length(nodes)
    A = zeros(Float64, (S, S))
    U = UnipartiteQuantitativeNetwork(A, String.(collect(values(nodes))))
    for arc in arcs
        U[string(nodes[string(arc[1])]), string(nodes[string(arc[2])])] = parse(Float64, arc[3])
    end
    return U
 end

 function pajek()
    data_path = joinpath(@__DIR__, "../..", "data", "pajek")
    files = readdir(data_path)
    filter!(f -> endswith(f, ".paj"), files)
    networks = Dict{Symbol,UnipartiteQuantitativeNetwork}()
    for f in files
       fname = Symbol(first(split(f, ".paj")))
       fpath = joinpath(data_path, f)
       networks[fname] = pajek(fpath)
    end
    return networks
 end


text:
@testset "SubSetSimulation" begin
    proposal = Normal()
    subset = SubSetSimulation(2000, 0.2, 10, proposal)

    @test isa(subset, SubSetSimulation)
    @test subset.n == 2000
    @test subset.target == 0.2
    @test subset.levels == 10
    @test subset.proposal == proposal

    @test_throws ErrorException("proposal must be a symmetric distribution") SubSetSimulation(
        2000, 0.2, 10, Exponential()
    )
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using Distributions

export ucb, poi, ei

"""
    ucb(x::Real, data; β::Real = 1)

The Upper Confidence Bound acquisition function.

The `β` is a hyperparameter.
Set `β` high to make the Bayesian Optimization more explorative
and set it lower to make it more exploitative.

Use this function as parameter for the `gaussian_process` function.

See also [`gaussian_process`](@ref), [`poi`](@ref), [`ei`](@ref).
"""
function ucb(x::Real, data; β::Real = 1)
    return func_est_μ(x, data) + β * func_est_σ(x, data)
end

"""
    poi(x::Real, data; τ::Real = maximum(data.points.y))

The Probability of Improvement acquisition function.

The `τ` is a hyperparameter.
This function describes the probability of improvement upon the function value `τ`.

Use this function as parameter for the `gaussian_process` function.

See also [`gaussian_process`](@ref), [`ucb`](@ref), [`ei`](@ref).
"""
function poi(x::Real, data; τ::Real = maximum(data.points.y))
    σ = func_est_σ(x, data)
    σ == 0 && return 0
    μ = func_est_μ(x, data)

    stdnorm = Normal()
    return cdf(stdnorm, (μ - τ) / σ)
end

"""
    ei(x::Real, data; τ::Real = maximum(data.points.y))

The Expected Improvement acquisition function.

The `τ` is a hyperparameter.
This function describes the expected improvement upon the function value `τ`.

Use this function as parameter for the `gaussian_process` function.

See also [`gaussian_process`](@ref), [`ucb`](@ref), [`poi`](@ref).
"""
function ei(x::Real, data; τ::Real = maximum(data.points.y))
    σ = func_est_σ(x, data)
    σ == 0 && return 0
    μ = func_est_μ(x, data)

    stdnorm = Normal()
    return (μ - τ) * cdf(stdnorm, (μ - τ) / σ) + σ * pdf(stdnorm, (μ - τ) / σ)
end


text:
function p54(data::Dict)

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
    @assert typeof(struc_el) == GenericSolid
  else
    println("No struc_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst
  nels = struc_el.nels
  nn = struc_el.nn

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  y_coords = zeros(nn)
  z_coords = zeros(nn)

  g_coord = zeros(ndim,nn)
  @assert :g_coord in keys(data)
  g_coord = data[:g_coord]'

  g_num = zeros(Int, fin_el.nod, nels)
  @assert :g_num in keys(data)
  g_num = reshape(data[:g_num]', fin_el.nod, nels)

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)
  println("There are $(neq) equations.\n")

  # Find global array sizes
  for iel in 1:nels
    num = g_num[:, iel]
    num_to_g!(num, nf, g)
    g_g[:, iel] = g
  end

  sample!(fin_el, points, weights)

  loads = OffsetArray(zeros(neq + 1), 0:neq)
  gravlo = OffsetArray(zeros(neq + 1), 0:neq)

  gsm = spzeros(neq, neq)
  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      beemat!(bee, deriv)
      km += (bee')*dee*bee*detm*weights[i]
      eld[nodof:nodof:ndof] += fun[:]*detm*weights[i]
    end
    fsparm!(gsm, g, km)
    gravlo[g] -= eld*prop[etype[iel], 3]
  end

  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end
  loads += gravlo

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      no[i] = nf[data[:fixed_freedoms][i][2], data[:fixed_freedoms][i][1]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i], no[i]] += penalty
      loads[no[i]] = gsm[no[i], no[i]] * value[i]
    end
  end

  cfgsm = cholesky(gsm)
  loads[1:neq] = cfgsm \ loads[1:neq]

  displacements = zeros(size(nf))
  for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if nf[i, j] > 0
        displacements[i,j] = loads[nf[i, j]]
      end
    end
  end
  displacements = displacements'

  if ndim == 3
    dis_df = DataFrame(
      x_disp = displacements[:, 1],
      y_disp = displacements[:, 2],
      z_disp = displacements[:, 3]
    )
  else
    dis_df = DataFrame(
      x_disp = displacements[:, 1],
      y_disp = displacements[:, 2]
    )
  end

  loads[0] = 0.0

  gc1 = Vector{Float64}()
  gc2 = Vector{Float64}()
  s1 = Vector{Float64}()
  s2 = Vector{Float64}()
  s3 = Vector{Float64}()

  if ndim ==  3
    gc3 = Vector{Float64}()
    t1 = Vector{Float64}()
    t2 = Vector{Float64}()
    t3 = Vector{Float64}()
  end

  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'
    g = g_g[:, iel]
    eld = loads[g]
    for i in 1:struc_el.nip
      shape_der!(der, points, i)
      shape_fun!(fun, points, i)
      gc = fun'*coord
      jac = inv(der*coord)
      deriv = jac*der
      beemat!(bee, deriv)
      sigma = dee*(bee*eld)
      gc1 = append!(gc1, gc[1])
      gc2 = append!(gc2, gc[2])
      s1 = append!(s1, sigma[1])
      s2 = append!(s2, sigma[2])
      s3 = append!(s3, sigma[3])
      if ndim == 3
        gc3 = append!(gc3, gc[3])
        t1 = append!(t1, sigma[4])
        t2 = append!(t2, sigma[5])
        t3 = append!(t3, sigma[6])
      end
    end
  end

  if ndim == 3
    sigma_df = DataFrame(
      x_coord = gc1,
      y_coord = gc2,
      z_coord = gc3,
      sig_x = s1,
      sig_y = s2,
      sig_z = s3,
      tau_xy = t1,
      tau_yz = t2,
      tau_zx = t3
    )
  else
    println("\nElement  x-coord   y-coord      sig_x        sig_y        tau_xy")
    sigma_df = DataFrame(
      x_coord = gc1,
      y_coord = gc2,
      sig_x = s1,
      sig_y = s2,
      tau_xy = s3
    )
  end

  fem = PtFEM.jFEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof,
    neq, penalty, etype, g, g_g, g_num, nf, no,
    node, num, sense, actions, bee, coord, gamma, dee,
    der, deriv, displacements, eld, fun, gc, g_coord, jac,
    km, mm, kg, cfgsm, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial)

    (fem, dis_df, sigma_df)
end

max_stars_repo_path: src/5 Elastic Solids/p54.jl

max_issues_repo_path: src/5 Elastic Solids/p54.jl

max_forks_repo_path: src/5 Elastic Solids/p54.jl


text:
using LinearAlgebra: axpy!
using SparseArrays

# exponential_utils.jl
# Contains functions related to the evaluation of scalar/matrix phi functions
# that are used by the exponential integrators.
#
# TODO: write a version of `exp!` that is non-allocating.

###################################################
# Dense algorithms
using LinearAlgebra: exp!

"""
    phi(z,k[;cache]) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the scalar phi functions for all orders up to k.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_k(z+1) = \\frac{\\varphi_k(z) - 1}{z}
```

Instead of using the recurrence relation, which is numerically unstable, a
formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software
package for computing matrix exponentials. ACM Transactions on Mathematical
Software (TOMS), 24(1), 130-156. Theorem 1).
"""
function phi(z::T, k::Integer; cache=nothing) where {T <: Number}
  # Construct the matrix
  if cache == nothing
    cache = fill(zero(T), k+1, k+1)
  else
    fill!(cache, zero(T))
  end
  cache[1,1] = z
  for i = 1:k
    cache[i,i+1] = one(T)
  end
  P = exp!(cache)
  return P[1,:]
end

"""
    phiv_dense(A,v,k[;cache]) -> [phi_0(A)v phi_1(A)v ... phi_k(A)v]

Compute the matrix-phi-vector products for small, dense `A`. `k`` >= 1.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_k(z+1) = \\frac{\\varphi_k(z) - 1}{z}
```

Instead of using the recurrence relation, which is numerically unstable, a
formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software
package for computing matrix exponentials. ACM Transactions on Mathematical
Software (TOMS), 24(1), 130-156. Theorem 1).
"""
function phiv_dense(A, v, k; cache=nothing)
  w = Matrix{eltype(A)}(undef, length(v), k+1)
  phiv_dense!(w, A, v, k; cache=cache)
end
"""
    phiv_dense!(w,A,v,k[;cache]) -> w

Non-allocating version of `phiv_dense`.
"""
function phiv_dense!(w::AbstractMatrix{T}, A::AbstractMatrix{T},
  v::AbstractVector{T}, k::Integer; cache=nothing) where {T <: Number}
  @assert size(w, 1) == size(A, 1) == size(A, 2) == length(v) "Dimension mismatch"
  @assert size(w, 2) == k+1 "Dimension mismatch"
  m = length(v)
  # Construct the extended matrix
  if cache == nothing
    cache = fill(zero(T), m+k, m+k)
  else
    @assert size(cache) == (m+k, m+k) "Dimension mismatch"
    fill!(cache, zero(T))
  end
  cache[1:m, 1:m] = A
  cache[1:m, m+1] = v
  for i = m+1:m+k-1
    cache[i, i+1] = one(T)
  end
  P = exp!(cache)
  # Extract results
  @views mul!(w[:, 1], P[1:m, 1:m], v)
  @inbounds for i = 1:k
    @inbounds for j = 1:m
      w[j, i+1] = P[j, m+i]
    end
  end
  return w
end

"""
    phi(A,k[;cache]) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the matrix phi functions for all orders up to k. `k` >= 1.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_k(z+1) = \\frac{\\varphi_k(z) - 1}{z}
```

Calls `phiv_dense` on each of the basis vectors to obtain the answer.
"""
function phi(A::AbstractMatrix{T}, k; caches=nothing) where {T <: Number}
  m = size(A, 1)
  out = [Matrix{T}(undef, m, m) for i = 1:k+1]
  phi!(out, A, k; caches=caches)
end
"""
    phi!(out,A,k[;caches]) -> out

Non-allocating version of `phi` for matrix inputs.
"""
function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches=nothing) where {T <: Number}
  m = size(A, 1)
  @assert length(out) == k + 1 && all(P -> size(P) == (m,m), out) "Dimension mismatch"
  if caches == nothing
    e = Vector{T}(undef, m)
    W = Matrix{T}(undef, m, k+1)
    C = Matrix{T}(undef, m+k, m+k)
  else
    e, W, C = caches
    @assert size(e) == (m,) && size(W) == (m, k+1) && size(C) == (m+k, m+k) "Dimension mismatch"
  end
  @inbounds for i = 1:m
    fill!(e, zero(T)); e[i] = one(T) # e is the ith basis vector
    phiv_dense!(W, A, e, k; cache=C) # W = [phi_0(A)*e phi_1(A)*e ... phi_k(A)*e]
    @inbounds for j = 1:k+1
      @inbounds for s = 1:m
        out[j][s, i] = W[s, j]
      end
    end
  end
  return out
end

##############################################
# Krylov algorithms
"""
    KrylovSubspace{T}(n,[maxiter=30]) -> Ks

Constructs an uninitialized Krylov subspace, which can be filled by `arnoldi!`.

The dimension of the subspace, `Ks.m`, can be dynamically altered but should
be smaller than `maxiter`, the maximum allowed arnoldi iterations.

    getV(Ks) -> V
    getH(Ks) -> H

Access methods for the (extended) orthonormal basis `V` and the (extended)
Gram-Schmidt coefficients `H`. Both methods return a view into the storage
arrays and has the correct dimensions as indicated by `Ks.m`.

    resize!(Ks, maxiter) -> Ks

Resize `Ks` to a different `maxiter`, destroying its contents.

This is an expensive operation and should be used scarsely.
"""
mutable struct KrylovSubspace{B, T}
  m::Int        # subspace dimension
  maxiter::Int  # maximum allowed subspace size
  beta::B       # norm(b,2)
  V::Matrix{T}  # orthonormal bases
  H::Matrix{T}  # Gram-Schmidt coefficients
  KrylovSubspace{T}(n::Integer, maxiter::Integer=30) where {T} = new{real(T), T}(
    maxiter, maxiter, zero(real(T)), Matrix{T}(undef, n, maxiter + 1),
    fill(zero(T), maxiter + 1, maxiter))
end
getH(Ks::KrylovSubspace) = @view(Ks.H[1:Ks.m + 1, 1:Ks.m])
getV(Ks::KrylovSubspace) = @view(Ks.V[:, 1:Ks.m + 1])
function Base.resize!(Ks::KrylovSubspace{B,T}, maxiter::Integer) where {B,T}
  V = Matrix{T}(undef, size(Ks.V, 1), maxiter + 1)
  H = fill(zero(T), maxiter + 1, maxiter)
  Ks.V = V; Ks.H = H
  Ks.m = Ks.maxiter = maxiter
  return Ks
end
function Base.show(io::IO, Ks::KrylovSubspace)
  println(io, "$(Ks.m)-dimensional Krylov subspace with fields")
  println(io, "beta: $(Ks.beta)")
  print(io, "V: ")
  println(IOContext(io, limit=true), getV(Ks))
  print(io, "H: ")
  println(IOContext(io, limit=true), getH(Ks))
end

"""
    arnoldi(A,b[;m,tol,opnorm,cache]) -> Ks

Performs `m` anoldi iterations to obtain the Krylov subspace K_m(A,b).

The n x (m + 1) basis vectors `getV(Ks)` and the (m + 1) x m upper Heisenberg
matrix `getH(Ks)` are related by the recurrence formula

```
v_1=b,\\quad Av_j = \\sum_{i=1}^{j+1}h_{ij}v_i\\quad(j = 1,2,\\ldots,m)
```

`iop` determines the length of the incomplete orthogonalization procedure [^1].
The default value of 0 indicates full Arnoldi. For symmetric/Hermitian `A`,
`iop` will be ignored and the Lanczos algorithm will be used instead.

Refer to `KrylovSubspace` for more information regarding the output.

Happy-breakdown occurs whenver `norm(v_j) < tol * opnorm(A, Inf)`, in this case
the dimension of `Ks` is smaller than `m`.

[^1]: Koskela, A. (2015). Approximating the matrix exponential of an
advection-diffusion operator using the incomplete orthogonalization method. In
Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353).
Springer, Cham.
"""
function arnoldi(A, b; m=min(30, size(A, 1)), tol=1e-7, opnorm=LinearAlgebra.opnorm,
  iop=0, cache=nothing)
  Ks = KrylovSubspace{eltype(b)}(length(b), m)
  arnoldi!(Ks, A, b; m=m, tol=tol, opnorm=opnorm, cache=cache, iop=iop)
end
"""
    arnoldi!(Ks,A,b[;tol,m,opnorm,cache]) -> Ks

Non-allocating version of `arnoldi`.
"""
function arnoldi!(Ks::KrylovSubspace{B, T}, A, b::AbstractVector{T}; tol::Real=1e-7,
  m::Int=min(Ks.maxiter, size(A, 1)), opnorm=LinearAlgebra.opnorm, iop::Int=0, cache=nothing) where {B, T <: Number}
  if ishermitian(A)
    return lanczos!(Ks, A, b; tol=tol, m=m, opnorm=opnorm, cache=cache)
  end
  if m > Ks.maxiter
    resize!(Ks, m)
  else
    Ks.m = m # might change if happy-breakdown occurs
  end
  V, H = getV(Ks), getH(Ks)
  vtol = tol * opnorm(A, Inf)
  if iop == 0
    iop = m
  end
  # Safe checks
  n = size(V, 1)
  @assert length(b) == size(A,1) == size(A,2) == n "Dimension mismatch"
  if cache == nothing
    cache = similar(b)
  else
    @assert size(cache) == (n,) "Dimension mismatch"
  end
  # Arnoldi iterations (with IOP)
  fill!(H, zero(T))
  Ks.beta = norm(b)
  @. V[:, 1] = b / Ks.beta
  @inbounds for j = 1:m
    mul!(cache, A, @view(V[:, j]))
    @inbounds for i = max(1, j - iop + 1):j
      alpha = dot(@view(V[:, i]), cache)
      H[i, j] = alpha
      axpy!(-alpha, @view(V[:, i]), cache)
    end
    beta = norm(cache)
    H[j+1, j] = beta
    @inbounds for i = 1:n
      V[i, j+1] = cache[i] / beta
    end
    if beta < vtol # happy-breakdown
      Ks.m = j
      break
    end
  end
  return Ks
end
"""
    lanczos!(Ks,A,b[;tol,m,opnorm,cache]) -> Ks

A variation of `arnoldi!` that uses the Lanczos algorithm for Hermitian matrices.
"""
function lanczos!(Ks::KrylovSubspace{B, T}, A, b::AbstractVector{T}; tol=1e-7,
  m=min(Ks.maxiter, size(A, 1)), opnorm=LinearAlgebra.opnorm, cache=nothing) where {B, T <: Number}
  if m > Ks.maxiter
    resize!(Ks, m)
  else
    Ks.m = m # might change if happy-breakdown occurs
  end
  V, H = getV(Ks), getH(Ks)
  vtol = tol * opnorm(A, Inf)
  # Safe checks
  n = size(V, 1)
  @assert length(b) == size(A,1) == size(A,2) == n "Dimension mismatch"
  if cache == nothing
    cache = similar(b)
  else
    @assert size(cache) == (n,) "Dimension mismatch"
  end
  # Lanczos iterations
  fill!(H, zero(T))
  Ks.beta = norm(b)
  @. V[:, 1] = b / Ks.beta
  @inbounds for j = 1:m
    vj = @view(V[:, j])
    mul!(cache, A, vj)
    alpha = dot(vj, cache)
    H[j, j] = alpha
    axpy!(-alpha, vj, cache)
    if j > 1
      axpy!(-H[j-1, j], @view(V[:, j-1]), cache)
    end
    beta = norm(cache)
    H[j+1, j] = beta
    if j < m
      H[j, j+1] = beta
    end
    @inbounds for i = 1:n
      V[i, j+1] = cache[i] / beta
    end
    if beta < vtol # happy-breakdown
      Ks.m = j
      break
    end
  end
  return Ks
end

# Cache type for expv
mutable struct ExpvCache{T}
  mem::Vector{T}
  ExpvCache{T}(maxiter::Int) where {T} = new{T}(Vector{T}(undef, maxiter^2))
end
function Base.resize!(C::ExpvCache{T}, maxiter::Int) where {T}
  C.mem = Vector{T}(undef, maxiter^2 * 2)
  return C
end
function get_cache(C::ExpvCache, m::Int)
  m^2 > length(C.mem) && resize!(C, m) # resize the cache if needed
  reshape(@view(C.mem[1:m^2]), m, m)
end
"""
    expv(t,A,b; kwargs) -> exp(tA)b

Compute the matrix-exponential-vector product using Krylov.

A Krylov subspace is constructed using `arnoldi` and `exp!` is called
on the Heisenberg matrix. Consult `arnoldi` for the values of the keyword
arguments.

    expv(t,Ks; cache) -> exp(tA)b

Compute the expv product using a pre-constructed Krylov subspace.
"""
function expv(t, A, b; m=min(30, size(A, 1)), tol=1e-7, opnorm=LinearAlgebra.opnorm, cache=nothing, iop=0)
  Ks = arnoldi(A, b; m=m, tol=tol, opnorm=opnorm, iop=iop)
  w = similar(b)
  expv!(w, t, Ks; cache=cache)
end
function expv(t, Ks::KrylovSubspace{B, T}; cache=nothing) where {B, T}
  n = size(getV(Ks), 1)
  w = Vector{T}(undef, n)
  expv!(w, t, Ks; cache=cache)
end
"""
    expv!(w,t,Ks[;cache]) -> w

Non-allocating version of `expv` that uses precomputed Krylov subspace `Ks`.
"""
function expv!(w::AbstractVector{T}, t::Number, Ks::KrylovSubspace{B, T};
  cache=nothing) where {B, T <: Number}
  m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
  @assert length(w) == size(V, 1) "Dimension mismatch"
  if cache == nothing
    cache = Matrix{T}(undef, m, m)
  elseif isa(cache, ExpvCache)
    cache = get_cache(cache, m)
  else
    throw(ArgumentError("Cache must be an ExpvCache"))
  end
  lmul!(t, copyto!(cache, @view(H[1:m, :])))
  if ishermitian(cache)
    # Optimize the case for symtridiagonal H
    F = eigen!(SymTridiagonal(cache))
    expHe = F.vectors * (exp.(F.values) .* @view(F.vectors[1, :]))
  else
    expH = exp!(cache)
    expHe = @view(expH[:, 1])
  end
  lmul!(beta, mul!(w, @view(V[:, 1:m]), expHe)) # exp(A) ≈ norm(b) * V * exp(H)e
end

# Cache type for phiv
mutable struct PhivCache{T}
  mem::Vector{T}
  function PhivCache{T}(maxiter::Int, p::Int) where {T}
    numelems = maxiter + maxiter^2 + (maxiter + p)^2 + maxiter*(p + 1)
    new{T}(Vector{T}(undef, numelems))
  end
end
function Base.resize!(C::PhivCache{T}, maxiter::Int, p::Int) where {T}
  numelems = maxiter + maxiter^2 + (maxiter + p)^2 + maxiter*(p + 1)
  C.mem = Vector{T}(undef, numelems * 2)
  return C
end
function get_caches(C::PhivCache, m::Int, p::Int)
  numelems = m + m^2 + (m + p)^2 + m*(p + 1)
  numelems^2 > length(C.mem) && resize!(C, m, p) # resize the cache if needed
  e = @view(C.mem[1:m]); offset = m
  Hcopy = reshape(@view(C.mem[offset + 1:offset + m^2]), m, m); offset += m^2
  C1 = reshape(@view(C.mem[offset + 1:offset + (m+p)^2]), m+p, m+p); offset += (m+p)^2
  C2 = reshape(@view(C.mem[offset + 1:offset + m*(p+1)]), m, p+1)
  return e, Hcopy, C1, C2
end
"""
    phiv(t,A,b,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using Krylov. `k` >= 1.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_k(z+1) = \\frac{\\varphi_k(z) - 1}{z}
```

A Krylov subspace is constructed using `arnoldi` and `phiv_dense` is called
on the Heisenberg matrix. If `correct=true`, then phi_0 through phi_k-1 are
updated using the last Arnoldi vector v_m+1 [^1]. If `errest=true` then an
additional error estimate for the second-to-last phi is also returned. For
the additional keyword arguments, consult `arnoldi`.

  phiv(t,Ks,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.

[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating
the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.
Formula (10).
"""
function phiv(t, A, b, k; m=min(30, size(A, 1)), tol=1e-7, opnorm=LinearAlgebra.opnorm, iop=0,
  cache=nothing, correct=false, errest=false)
  Ks = arnoldi(A, b; m=m, tol=tol, opnorm=opnorm, iop=iop)
  w = Matrix{eltype(b)}(undef, length(b), k+1)
  phiv!(w, t, Ks, k; cache=cache, correct=correct, errest=errest)
end
function phiv(t, Ks::KrylovSubspace{B, T}, k; cache=nothing, correct=false,
  errest=false) where {B, T}
  n = size(getV(Ks), 1)
  w = Matrix{T}(undef, n, k+1)
  phiv!(w, t, Ks, k; cache=cache, correct=correct, errest=errest)
end
"""
    phiv!(w,t,Ks,k[;cache,correct,errest]) -> w[,errest]

Non-allocating version of 'phiv' that uses precomputed Krylov subspace `Ks`.
"""
function phiv!(w::AbstractMatrix{T}, t::Number, Ks::KrylovSubspace{B, T}, k::Integer;
  cache=nothing, correct=false, errest=false) where {B, T <: Number}
  m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
  @assert size(w, 1) == size(V, 1) "Dimension mismatch"
  @assert size(w, 2) == k + 1 "Dimension mismatch"
  if cache == nothing
    cache = PhivCache{T}(m, k)
  elseif !isa(cache, PhivCache)
    throw(ArgumentError("Cache must be a PhivCache"))
  end
  e, Hcopy, C1, C2 = get_caches(cache, m, k)
  lmul!(t, copyto!(Hcopy, @view(H[1:m, :])))
  fill!(e, zero(T)); e[1] = one(T) # e is the [1,0,...,0] basis vector
  phiv_dense!(C2, Hcopy, e, k; cache=C1) # C2 = [ϕ0(H)e ϕ1(H)e ... ϕk(H)e]
  lmul!(beta, mul!(w, @view(V[:, 1:m]), C2)) # f(A) ≈ norm(b) * V * f(H)e
  if correct
    # Use the last Arnoldi vector for correction with little additional cost
    # correct_p = beta * h_{m+1,m} * (em^T phi_p+1(H) e1) * v_m+1
    betah = beta * H[end,end] * t
    vlast = @view(V[:,end])
    @inbounds for i = 1:k
      axpy!(betah * C2[end, i+1], vlast, @view(w[:, i]))
    end
  end
  if errest
    err = abs(beta * H[end, end] * t * C2[end, end])
    return w, err
  else
    return w
  end
end

###########################################
# Krylov phiv with internal time-stepping
"""
    exp_timestep(ts,A,b[;adaptive,tol,kwargs...]) -> U

Evaluates the matrix exponentiation-vector product using time stepping

```math
u = \\exp(tA)b
```

`ts`` is an array of time snapshots for u, with `U[:,j] ≈ u(ts[j])`. `ts` can
also be just one value, in which case only the end result is returned and `U`
is a vector.

The time stepping formula of Niesen & Wright is used [^1]. If the time step
`tau` is not specified, it is chosen according to (17) of Neisen & Wright. If
`adaptive==true`, the time step and Krylov subsapce size adaptation scheme of
Niesen & Wright is used, the relative tolerance of which can be set using the
keyword parameter `tol`. The delta and gamma parameter of the adaptation
scheme can also be adjusted.

Set `verbose=true` to print out the internal steps (for debugging). For the
other keyword arguments, consult `arnoldi` and `phiv`, which are used
internally.

Note that this function is just a special case of `phiv_timestep` with a more
intuitive interface (vector `b` instead of a n-by-1 matrix `B`).

[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for
evaluating the φ-functions in exponential integrators. arXiv preprint
arXiv:0907.4631.
"""
function expv_timestep(ts::Vector{tType}, A, b; kwargs...) where {tType <: Real}
  U = Matrix{eltype(A)}(undef, size(A, 1), length(ts))
  expv_timestep!(U, ts, A, b; kwargs...)
end
function expv_timestep(t::tType, A, b; kwargs...) where {tType <: Real}
  u = Vector{eltype(A)}(undef, size(A, 1))
  expv_timestep!(u, t, A, b; kwargs...)
end
"""
    expv_timestep!(u,t,A,b[;kwargs]) -> u

Non-allocating version of `expv_timestep`.
"""
function expv_timestep!(u::AbstractVector{T}, t::tType, A, b::AbstractVector{T};
  kwargs...) where {T <: Number, tType <: Real}
  expv_timestep!(reshape(u, length(u), 1), [t], A, b; kwargs...)
  return u
end
function expv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, b::AbstractVector{T};
  kwargs...) where {T <: Number, tType <: Real}
  B = reshape(b, length(b), 1)
  phiv_timestep!(U, ts, A, B; kwargs...)
end
"""
    phiv_timestep(ts,A,B[;adaptive,tol,kwargs...]) -> U

Evaluates the linear combination of phi-vector products using time stepping

```math
u = \\varphi_0(tA)b_0 + t\\varphi_1(tA)b_1 + \\cdots + t^p\\varphi_p(tA)b_p
```

`ts`` is an array of time snapshots for u, with `U[:,j] ≈ u(ts[j])`. `ts` can
also be just one value, in which case only the end result is returned and `U`
is a vector.

The time stepping formula of Niesen & Wright is used [^1]. If the time step
`tau` is not specified, it is chosen according to (17) of Neisen & Wright. If
`adaptive==true`, the time step and Krylov subsapce size adaptation scheme of
Niesen & Wright is used, the relative tolerance of which can be set using the
keyword parameter `tol`. The delta and gamma parameter of the adaptation
scheme can also be adjusted.

Set `verbose=true` to print out the internal steps (for debugging). For the
other keyword arguments, consult `arnoldi` and `phiv`, which are used
internally.

[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for
evaluating the φ-functions in exponential integrators. arXiv preprint
arXiv:0907.4631.
"""
function phiv_timestep(ts::Vector{tType}, A, B; kwargs...) where {tType <: Real}
  U = Matrix{eltype(A)}(undef, size(A, 1), length(ts))
  phiv_timestep!(U, ts, A, B; kwargs...)
end
function phiv_timestep(t::tType, A, B; kwargs...) where {tType <: Real}
  u = Vector{eltype(A)}(undef, size(A, 1))
  phiv_timestep!(u, t, A, B; kwargs...)
end
"""
    phiv_timestep!(U,ts,A,B[;kwargs]) -> U

Non-allocating version of `phiv_timestep`.
"""
function phiv_timestep!(u::AbstractVector{T}, t::tType, A, B::AbstractMatrix{T};
  kwargs...) where {T <: Number, tType <: Real}
  phiv_timestep!(reshape(u, length(u), 1), [t], A, B; kwargs...)
  return u
end
function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractMatrix{T}; tau::Real=0.0,
  m::Int=min(10, size(A, 1)), tol::Real=1e-7, opnorm=LinearAlgebra.opnorm, iop::Int=0,
  correct::Bool=false, caches=nothing, adaptive=false, delta::Real=1.2,
  gamma::Real=0.8, NA::Int=0, verbose=false) where {T <: Number, tType <: Real}
  # Choose initial timestep
  abstol = tol * opnorm(A, Inf)
  verbose && println("Absolute tolerance: $abstol")
  if iszero(tau)
    Anorm = opnorm(A, Inf)
    b0norm = norm(@view(B[:, 1]), Inf)
    tau = 10/Anorm * (abstol * ((m+1)/ℯ)^(m+1) * sqrt(2*pi*(m+1)) /
      (4*Anorm*b0norm))^(1/m)
    verbose && println("Initial time step unspecified, chosen to be $tau")
  end
  # Initialization
  n = size(U, 1)
  sort!(ts); tend = ts[end]
  p = size(B, 2) - 1
  @assert length(ts) == size(U, 2) "Dimension mismatch"
  @assert n == size(A, 1) == size(A, 2) == size(B, 1) "Dimension mismatch"
  if caches == nothing
    u = Vector{T}(undef, n)              # stores the current state
    W = Matrix{T}(undef, n, p+1)         # stores the w vectors
    P = Matrix{T}(undef, n, p+2)         # stores output from phiv!
    Ks = KrylovSubspace{T}(n, m)  # stores output from arnoldi!
    phiv_cache = nothing         # cache used by phiv!
  else
    u, W, P, Ks, phiv_cache = caches
    @assert length(u) == n && size(W, 1) == n && size(P, 1) == n "Dimension mismatch"
    # W and P may be bigger than actually needed
    W = @view(W[:, 1:p+1])
    P = @view(P[:, 1:p+2])
  end
  copyto!(u, @view(B[:, 1])) # u(0) = b0
  coeffs = ones(tType, p);
  if adaptive # initialization step for the adaptive scheme
    if ishermitian(A)
      iop = 2 # does not have an effect on arnoldi!, just for flops estimation
    end
    if iszero(NA)
      _A = convert(AbstractMatrix, A)
      if isa(_A, SparseMatrixCSC)
        NA = nnz(_A)
      else
        NA = count(!iszero, _A) # not constant operation, should be best avoided
      end
    end
  end

  t = 0.0       # current time
  snapshot = 1  # which snapshot to compute next
  while t < tend # time stepping loop
    if t + tau > tend # last step
      tau = tend - t
    end
    # Part 1: compute w0...wp using the recurrence relation (16)
    copyto!(@view(W[:, 1]), u) # w0 = u(t)
    @inbounds for l = 1:p-1 # compute cl = t^l/l!
      coeffs[l+1] = coeffs[l] * t / l
    end
    @views @inbounds for j = 1:p
      mul!(W[:, j+1], A, W[:, j])
      for l = 0:p-j
        axpy!(coeffs[l+1], B[:, j+l+1], W[:, j+1])
      end
    end
    # Part 2: compute ϕp(tau*A)wp using Krylov, possibly with adaptation
    arnoldi!(Ks, A, @view(W[:, end]); tol=tol, m=m, opnorm=opnorm, iop=iop, cache=u)
    _, epsilon = phiv!(P, tau, Ks, p + 1; cache=phiv_cache, correct=correct, errest=true)
    verbose && println("t = $t, m = $m, tau = $tau, error estimate = $epsilon")
    if adaptive
      omega = (tend / tau) * (epsilon / abstol)
      epsilon_old = epsilon; m_old = m; tau_old = tau
      q = m/4; kappa = 2.0; maxtau = tend - t
      while omega > delta # inner loop of Algorithm 3
        m_new, tau_new, q, kappa = _phiv_timestep_adapt(
          m, tau, epsilon, m_old, tau_old, epsilon_old, q, kappa,
          gamma, omega, maxtau, n, p, NA, iop, opnorm(getH(Ks), 1), verbose)
        m, m_old = m_new, m
        tau, tau_old = tau_new, tau
        # Compute ϕp(tau*A)wp using the new parameters
        arnoldi!(Ks, A, @view(W[:, end]); tol=tol, m=m, opnorm=opnorm, iop=iop, cache=u)
        _, epsilon_new = phiv!(P, tau, Ks, p + 1; cache=phiv_cache, correct=correct, errest=true)
        epsilon, epsilon_old = epsilon_new, epsilon
        omega = (tend / tau) * (epsilon / abstol)
        verbose && println("  * m = $m, tau = $tau, error estimate = $epsilon")
      end
    end
    # Part 3: update u using (15)
    lmul!(tau^p, copyto!(u, @view(P[:, end - 1])))
    @inbounds for l = 1:p-1 # compute cl = tau^l/l!
      coeffs[l+1] = coeffs[l] * tau / l
    end
    @views @inbounds for j = 0:p-1
      axpy!(coeffs[j+1], W[:, j+1], u)
    end
    # Fill out all snapshots in between the current step
    while snapshot <= length(ts) && t + tau >= ts[snapshot]
      tau_snapshot = ts[snapshot] - t
      u_snapshot = @view(U[:, snapshot])
      phiv!(P, tau_snapshot, Ks, p + 1; cache=phiv_cache, correct=correct)
      lmul!(tau_snapshot^p, copyto!(u_snapshot, @view(P[:, end - 1])))
      @inbounds for l = 1:p-1 # compute cl = tau^l/l!
        coeffs[l+1] = coeffs[l] * tau_snapshot / l
      end
      @views @inbounds for j = 0:p-1
        axpy!(coeffs[j+1], W[:, j+1], u_snapshot)
      end
      snapshot += 1
    end

    t += tau
  end

  return U
end
# Helper functions for phiv_timestep!
function _phiv_timestep_adapt(m, tau, epsilon, m_old, tau_old, epsilon_old, q, kappa,
  gamma, omega, maxtau, n, p, NA, iop, Hnorm, verbose)
  # Compute new m and tau (Algorithm 4)
  if tau_old > tau
    q = log(tau/tau_old) / log(epsilon/epsilon_old) - 1
  end # else keep q the same
  tau_new = tau * (gamma / omega)^(1/(q + 1))
  tau_new = min(max(tau_new, tau/5), 2*tau, maxtau)
  if m_old < m
    kappa = (epsilon/epsilon_old)^(1/(m_old - m))
  end # else keep kappa the same
  m_new = m + ceil(Int, log(omega / gamma) / log(kappa))
  m_new = min(max(m_new, div(3*m, 4), 1), Int(ceil(4*m / 3)))
  verbose && println("  - Proposed new m: $m_new, new tau: $tau_new")
  # Compare costs of using new m vs new tau (23)
  cost_tau = _phiv_timestep_estimate_flops(m, tau_new, n, p, NA, iop, Hnorm, maxtau)
  cost_m = _phiv_timestep_estimate_flops(m_new, tau, n, p, NA, iop, Hnorm, maxtau)
  verbose && println("  - Cost to use new m: $cost_m flops, new tau: $cost_tau flops")
  if cost_tau < cost_m
    m_new = m
  else
    tau_new = tau
  end
  return m_new, tau_new, q, kappa
end
function _phiv_timestep_estimate_flops(m, tau, n, p, NA, iop, Hnorm, maxtau)
  # Estimate flops for the update of W and u
  flops_W = 2 * (p - 1) * (NA + n)
  flops_u = (2 * p + 1) * n
  # Estimate flops for arnoldi!
  if iop == 0
    iop = m
  end
  flops_matvec = 2 * m * NA
  flops_vecvec = 0
  for i = 1:m
    flops_vecvec += 3 * min(i, iop)
  end
  # Estimate flops for phiv! (7)
  MH = 44/3 + 2 * ceil(max(0.0, log2(Hnorm / 5.37)))
  flops_phiv = round(Int, MH * (m + p)^3)

  flops_onestep = flops_W + flops_u + flops_matvec + flops_vecvec + flops_phiv
  return flops_onestep * Int(ceil(maxtau / tau))
end
function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int)
  n = length(u_prototype); T = eltype(u_prototype)
  u = similar(u_prototype)                      # stores the current state
  W = Matrix{T}(undef, n, p+1)                         # stores the w vectors
  P = Matrix{T}(undef, n, p+2)                         # stores output from phiv!
  Ks = KrylovSubspace{T}(n, maxiter)            # stores output from arnoldi!
  phiv_cache = PhivCache{T}(maxiter, p+1)       # cache used by phiv! (need +1 for error estimation)
  return u, W, P, Ks, phiv_cache
end


text:
# Snooping functions
function _snoopi_bot(snoop_script, tmin)
    return quote
        using SnoopCompileCore

        __the_data__ = @snoopi tmin=$tmin begin
            $snoop_script
        end
    end
end

function _snoopc_bot(snoop_script)
    return quote
        using SnoopCompileCore

        @snoopc "compiles.log" begin
            $snoop_script
        end

        using SnoopCompile

        __the_data__ = SnoopCompile.read("compiles.log")[2]
        Base.rm("compiles.log", force = true)
    end
end

function _snoop_analysis_bot(snooping_code, package_name, precompile_folder, subst, exclusions, check_eval)
    return quote

        ################################################################
        @info "Processsing the generated precompile signatures"

        using SnoopCompile

        ### Parse the compiles and generate precompilation scripts
        let packageSym = Symbol($package_name),
            pc = SnoopCompile.parcel(__the_data__; subst = $subst, exclusions = $exclusions, check_eval = $check_eval)
            if !haskey(pc, packageSym)
                @warn "no precompile signature is found for $($package_name). Don't load the package before snooping. Restart your Julia session."
                if !isdir($precompile_folder)
                    mkpath($precompile_folder)
                end
                Base.write("$($precompile_folder)/precompile_$($package_name).jl", """
                function _precompile_()
                    # nothing
                end
                """)
            else # if any precompilation script is generated
                onlypackage = Dict( packageSym => Base.sort(pc[packageSym]) )
                SnoopCompile.write($precompile_folder, onlypackage)
                @info "precompile signatures were written to $($precompile_folder)"
            end
        end
    end
end

function _snoop_bot_expr(config::BotConfig, snoop_script, test_modul::Module; snoop_mode::Symbol)
    package_name = config.package_name
    exclusions = config.exclusions
    check_eval = config.check_eval
    os = config.os
    else_os = config.else_os
    version = config.version
    else_version = config.else_version
    package_path = config.package_path
    precompiles_rootpath = config.precompiles_rootpath
    subst = config.subst
    tmin = config.tmin

    if test_modul != Main && string(test_modul) == package_name
        @error "Your example/test shouldn't be in the same module!. Use `Main` instead."
    end

    ################################################################
    package_rootpath = dirname(dirname(package_path))

    new_includer_file(package_name, package_path, precompiles_rootpath, os, else_os, version, else_version) # create an precompile includer file
    add_includer(package_name, package_path) # add the code to packages source for including the includer

    # precompile folder for writing
    if  isnothing(os)
        if isnothing(version)
            precompile_folder = precompiles_rootpath
        else
            precompile_folder = "$precompiles_rootpath/$(VersionFloat(VERSION))"
        end
    else
        if isnothing(version)
            precompile_folder = "$precompiles_rootpath/$(detectOS()[1])"
        else
            precompile_folder = "$precompiles_rootpath/$(detectOS()[1])/$(VersionFloat(VERSION))"
        end
    end

    # automatic (based on Julia version)
    if snoop_mode == :auto
        if VERSION < v"1.2"
            snoop_mode = :snoopc
        else
            snoop_mode = :snoopi
        end
    end

    if snoop_mode == :snoopi
        snooping_code = _snoopi_bot(snoop_script, tmin)
    elseif snoop_mode == :snoopc
        snooping_code = _snoopc_bot(snoop_script)
    else
        error("snoop_mode $snoop_mode is unkown")
    end
    snooping_code = toplevel_string(snooping_code)

    analysis_code = _snoop_analysis_bot(snooping_code, package_name, precompile_folder, subst, exclusions, check_eval)
    analysis_code = toplevel_string(analysis_code)

    snooping_analysis_code = "$snooping_code; $analysis_code;"

    if in(get(ENV, "SnoopCompile_coverage_ENV", false), [true, "true"])
        code_coverage = "user"
    else
        code_coverage = "none"
    end

    julia_cmd = `$(Base.julia_cmd()) --code-coverage=$code_coverage --project=@. -e $snooping_analysis_code`

    addpkg_ifnotfound(:SnoopCompileCore, test_modul)
    addpkg_ifnotfound(:SnoopCompile, test_modul)
    out = quote
        ################################################################
        using CompileBot

        # Environment variable to detect SnoopCompile bot
        global SnoopCompile_ENV = true

        ################################################################
        CompileBot.precompile_deactivator($package_path);
        ################################################################

        ### Log the compiles and analyze the compiles
        run($julia_cmd)

        ################################################################
        CompileBot.precompile_activator($package_path)

        global SnoopCompile_ENV = false
    end
    return out
end

function _snoop_bot_expr(config::BotConfig, test_modul::Module; snoop_mode::Symbol)

    package_name = config.package_name
    package_rootpath = dirname(dirname(pathof_noload(package_name)))
    runtestpath = "$package_rootpath/test/runtests.jl"

    package = Symbol(package_name)
    snoop_script = quote
        using $(package)
        include($runtestpath)
    end
    return _snoop_bot_expr(config, snoop_script, test_modul; snoop_mode = snoop_mode)
end

################################################################

"""
    snoop_bot(config::BotConfig, path_to_example_script::String, test_modul=Main; snoop_mode=:auto)

Generate precompile statements using a precompile script.
`config` can be generated by [`BotConfig`](@ref).
`path_to_example_script` is preferred to be an absolute path.
The example script will be run in the module specified by `test_modul`.
`snoop_mode` can be `:auto`, `:snoopi` (to run [`SnoopCompileCore.@snoopi`](@ref)),
or `:snoopc` (to run [`SnoopCompileCore.@snoopc`](@ref)),
where `:auto` chooses `:snoopi` on supported versions of Julia.

See the [online documentation](https://timholy.github.io/SnoopCompile.jl/stable/bot/)
for a more complete overview.

# Extended help

## Example

In this case, the bot-running script is placed in the same directory as the
precompile script, so we can use `@__DIR__` to find it:

```julia
using CompileBot

snoop_bot(BotConfig("MatLang"), "\$(@__DIR__)/example_script.jl")
```
"""
function snoop_bot(config::BotConfig, path_to_example_script::String, test_modul::Module=Main; snoop_mode::Symbol=:auto)
    # search for the script! - needed because of confusing paths when referencing pattern_or_file in CI
    path_to_example_script = searchdirsboth([pwd(),dirname(dirname(config.package_path))], path_to_example_script)
    snoop_script = quote
        include($path_to_example_script)
    end
    out =  _snoop_bot_expr(config, snoop_script, test_modul; snoop_mode=snoop_mode)
    Core.eval( test_modul, out )
end

################################################################
"""
    snoop_bot(config::BotConfig, test_modul::Module = Main; snoop_mode::Symbol = :auto)

Generate precompile statements using the package's `runtests.jl` file.

During snooping, `snoop_bot` sets the global variable `SnoopCompile_ENV` to `true`.
If needed, your `runtests.jl` can check for the existence and value of this variable to
customize test behavior specifically for snooping.
"""
function snoop_bot(config::BotConfig, test_modul::Module = Main; snoop_mode::Symbol = :auto)
    out = _snoop_bot_expr(config, test_modul; snoop_mode = snoop_mode)
    Core.eval( test_modul, out )
end

################################################################

"""
    snoop_bot(config::BotConfig, expression::Expr, test_modul::Module = Main; snoop_mode::Symbol = :auto)

Generate precompile statements by evaluating an expression, for example `:(using MyPackage)`.
Interpolation and macros are not supported.
"""
function snoop_bot(config::BotConfig, snoop_script::Expr, test_modul::Module  = Main; snoop_mode::Symbol = :auto)
    out = _snoop_bot_expr(config, snoop_script, test_modul; snoop_mode = snoop_mode)
    Core.eval( test_modul, out )
end