text:
"""
Main module for `SCARGC.jl`, the Julia implementation of SCARGC algorithm.

From this module, two functions are exported for public use:

- [`scargc_1NN`](@ref). The 1-Nearest Neighbor SCARGC implementation.
- [`extractValuesFromFile`](@ref). A functio to extract the data values from a dataset file.
"""
module SCARGC

using Distances, Statistics
using PyCall

export scargc_1NN, extractValuesFromFile

"""
    scargc_1NN(
        dataset         -> "dataset used in the algorighm"
        percentTraining -> "amount of data that is goung to be used as training data"
        maxPoolSize     -> "maximum instances that the pool size can store"
        K               -> "number of clusters"
    )

SCARGC implementation with the Nearest Neighbor classifier for label predicting.
The function prints the final accuracy and returns the vector with all predicted labels and the final accuracy.

The function starts getting the labeled and unlabeled and, with them, it creates the initial centroids.
Then, a loop starts over the unlabeled data, storing the instance and the predicted label (predicted with the
classification model).

When the pool reaches the maximum, represented by `maxPoolSize`, a clustering step is
made on the data stored in the pool to get the centroids from the current iteration (represented as `tempCurrentCentroids`
before receiving the labels and `currentCentroids` after it). With the centroids from the past iteration (represented
as `centroids`), we find the current iteration's centroids' labels and create the intermediary centroids, represented
as `intermed`, that stores the median between the past and current centroids and the current iretation's
centroid's labels (the reason to store that is to store the drift between the past and the current iteration).

After getting the labels, the past iteration's centroids receive the values stored in `intermed` and a new labels are
found using both centroids from current and past iterations. These labels are going to be useful to get the concordance
between them and the labels stored in the pool to know if the model is going to be updated.
The last part, after calculating the concordance, is update the classification model if the concordance is different of
100%.

For example, let's predict some data in the **1CHT** dataset.

We start loading the dataset and reading the values from it.

```@example
using SCARGC

# loading dataset
dataset = "datasets/synthetic/1CHT.txt"

# extracting values from the dataset file and storing into `data`
data = extractValuesFromFile(dataset, 16000, 3)

# printing the read dataset
prinln(data)

# predicting labels and storing them in `predictedLabels` and the accuracy in `accuracy`
predictedLabels, accuracy = scargc_1NN(data, 0.3125, 300, 2)
```

If we print the results, we're gonna have

```julia
# predicted labels
11000-element Array{Int64,1}:
 1
 2
 1
 1
 1
 2
 1
 1
 ⋮
 2
 1
 2
 1
 2
 2
 2

# accuracy ~
99.73040752351096
```

The accuracy happens if you're using a dataset with the actual labels of all instances. Then
the function can compare the prediction with the actual label to get the accuracy.
"""
function scargc_1NN(dataset::Array{T, 2} where {T<:Number}, percentTraining::Float64, maxPoolSize::Int64, K::Int64)
    labeledData, labeledDataLabels, streamData, streamLabels, features = fitData(dataset, percentTraining)
    centroids = findCentroids(features, labeledData, labeledDataLabels, K)

    predictedLabels = zeros(length(streamLabels))

    accuracyVector = zeros(size(streamData, 1))
    poolData = zeros(maxPoolSize, size(streamData, 2) + 1)
    poolDataIndex = 1

    for stream = 1:size(streamData, 1)
        testInstance = streamData[stream, :]
        currentLabel = streamLabels[stream]

        predictedLabel, _ = knnClassification(labeledData, labeledDataLabels, testInstance)

        poolData[poolDataIndex, 1:end-1] = testInstance
        poolData[poolDataIndex, end] = predictedLabel
        poolDataIndex += 1

        if poolDataIndex > maxPoolSize
            py"""
            from sklearn.cluster import KMeans
            import numpy as np
            import sys, os
            import warnings

            warnings.filterwarnings("ignore")

            def kmeans(X, k, init, centroids=[]):
                if (init == 0):
                    return KMeans(n_clusters=k).fit(X).cluster_centers_
                else:
                    centroids = np.array(centroids)
                    sys.stdout = open(os.devnull, 'w')
                    return KMeans(n_clusters=k, init=centroids).fit(X).cluster_centers_"""

            tempCurrentCentroids = py"kmeans"(poolData[:, 1:end-1], K, 1, centroids[end-K+1:end, 1:end-1])

            intermed, centroidLabels = findLabelForCurrentCentroids(centroids, tempCurrentCentroids)

            currentCentroids = zeros(size(tempCurrentCentroids, 1), size(tempCurrentCentroids, 2) + 1)
            currentCentroids[:, 1:end-1] = tempCurrentCentroids
            currentCentroids[:, end] = centroidLabels

            centroids = intermed

            newLabeledData = calculateNewLabeledData(poolData, currentCentroids, centroids)

            concordantLabelCount = count(label -> label == 1, poolData[:, end] .== newLabeledData)

            centroids, labeledData, labeledDataLabels = updateInformation(poolData, centroids, labeledData, labeledDataLabels,
                                                                        newLabeledData, currentCentroids, intermed,
                                                                        concordantLabelCount, maxPoolSize)

            poolData = zeros(maxPoolSize, size(streamData, 2) + 1)
            poolDataIndex = 1
        end

        predictedLabels[stream] = predictedLabel

        if predictedLabel == currentLabel
            global accuracyVector[stream] = 1
        end
    end

    finalAccuracy = (sum(accuracyVector)/size(streamData, 1)) * 100

    #println(finalAccuracy)

    return predictedLabels, finalAccuracy
end

"""
    extractValuesFromFile(
        fileName -> "name of the file that you're trying to read"
        rows     -> "number of rows in the file"
        columns  -> "number of columns in the file"
    )

The function reads a file and creates a matrix with the file's values.
Using the function `readuntil` and `parse` we can converte each value, separeted by comma.

The function returns a matrix with the values in the file.

For example, let's suppose we have a text file, called "example.txt", with **10** lines and **7** columns.

```
2.91120, 2.94941, 2.92468, 0.88324, 6.00512, 5.88021, 3
4.97443, 2.31345, 5.21826, 8.13938, 7.47390, 8.24444, 1
3.13772, 7.87192, 8.35722, 7.41002, 6.95894, 5.20002, 2
3.74301, 6.39082, 3.89166, 8.76716, 1.66374, 6.75246, 3
8.98766, 4.42780, 4.94278, 7.05217, 5.21106, 3.69790, 1
4.38337, 4.91497, 6.16593, 8.73148, 6.17557, 0.90185, 2
0.97263, 2.70122, 0.00343, 7.20105, 1.63296, 8.26112, 1
0.50329, 3.54209, 8.47787, 2.59826, 3.93825, 2.58200, 3
6.75223, 4.59601, 2.02472, 8.10523, 3.65602, 6.91874, 1
8.98924, 2.33177, 8.34892, 4.03544, 8.57646, 7.93690, 1
```

We're gonna run this function with the following values:

- `fileName` is going to receive the path to the file. In this case, it's just "example.txt";
- `rows` is going to receive 10 and
- `columns` is going to receive 7

Then, the function starts and reads each value separated by "," and stores these values into
a matrix. This matrix is returned in the end of the function, in the format:

```julia
10×7 Array{Float64,2}:
2.91120, 2.94941, 2.92468, 0.88324, 6.00512, 5.88021, 3
4.97443, 2.31345, 5.21826, 8.13938, 7.47390, 8.24444, 1
3.13772, 7.87192, 8.35722, 7.41002, 6.95894, 5.20002, 2
3.74301, 6.39082, 3.89166, 8.76716, 1.66374, 6.75246, 3
8.98766, 4.42780, 4.94278, 7.05217, 5.21106, 3.69790, 1
4.38337, 4.91497, 6.16593, 8.73148, 6.17557, 0.90185, 2
0.97263, 2.70122, 0.00343, 7.20105, 1.63296, 8.26112, 1
0.50329, 3.54209, 8.47787, 2.59826, 3.93825, 2.58200, 3
6.75223, 4.59601, 2.02472, 8.10523, 3.65602, 6.91874, 1
8.98924, 2.33177, 8.34892, 4.03544, 8.57646, 7.93690, 1
```

This matrix is then ready to be used along the code.

"""
function extractValuesFromFile(fileName::String, rows::Int64, columns::Int64)
    file = open(fileName)

    fileValues = zeros(rows, columns)

    for row = 1:rows
        for column = 1:columns-1
            global fileValues[row, column] = parse(Float64, readuntil(file, ','))
        end

        global fileValues[row, columns] = parse(Float64, readuntil(file, '\n'))
    end

    close(file)

    return fileValues
end

"""
    fitData(
        dataset         -> "the whole dataset loaded from file system"
        percentTraining -> "percentage of the amount of rows to be used as labeled data"
    )

Divides the dataset into labeled data and unlabeled data, training and testing data respectively.
These arrays are gonna be used in the whole program, so the fitting part is considered extremelly important because, if
something's wrong in the this step, the whole application is gonna fail.

The return of this function has the following variables:
    - labeledData       -> "`percentTraining`% of the dataset features"
    - labeledDataLabels -> "also `percentTraining`% of the dataset labels"
    - streamData        -> "remaining data of the dataset"
    - streamLabels      -> "remaining labels of the dataset"
    - features          -> "number of features present in the data"
"""
function fitData(dataset::Array{T, 2} where {T<:Number}, percentTraining::Float64)
    rows = size(dataset, 1)
    columns = size(dataset, 2)

    labeledRowCount = Int(ceil((percentTraining/100) * rows))

    labeledData = dataset[1:labeledRowCount, 1:end-1]
    labeledDataLabels = dataset[1:labeledRowCount, end]

    streamData = dataset[labeledRowCount + 1:end, 1:end-1]
    streamLabels = dataset[labeledRowCount + 1:end, end]

    features = columns - 1

    return labeledData, labeledDataLabels, streamData, streamLabels, features
end

"""
    knnClassification(
        labeledData  -> "labeled data used to apply the euclidean distance"
        labels       -> "all data labels"
        testInstance -> "the instance we want to predict the output based on the nearest neighbor."
    )

Applies the Nearest Neighbor to the `testInstance`.
The function calculates, for each `labeledData`, the Euclidean Distance between it and the test instance and.
If it's value is smaller then the smallest distance, the smallest distance recieves this distance value and the label
value, at this index, is stored.
The return of the function is the output and the data from the nearest neighbor.
"""
function knnClassification(labeledData::Array{T, 2} where {T<:Number}, labels::Array{T} where {T<:Number},
                           testInstance::Array{T} where {T<:Number})
    output = nothing
    resultData = nothing
    smallerDistance = Inf

    for row = 1:size(labeledData, 1)
        data = labeledData[row, :]

        distance = Distances.euclidean(testInstance, data)

        if (distance < smallerDistance)
            smallerDistance = distance
            output = labels[row]
            resultData = data
        end
    end

    return output, resultData
end

"""
    findCentroids(
        features          -> "number of features"
        labeledData       -> "labeled data used in classification"
        labeledDataLabels -> "labels of the labeled data used in classification"
        K                 -> "number of clusters"
    )

Finds the initial centroids used in program, using KMeans if necessary, with the data already fit.
The function has two possible ways:
    - The value of K is the same as the number of classes; and
    - The value of K is different as the number of classes.

If the values are the same, the function creates the centroids as the median of the data. In other words, if these
numbers are equal, the centroids have the value of the median of each feature.
Otherwise, if the values aren't the same, the KMeans method is used to define the initial centroids.

The function returns the centroids array.
"""
function findCentroids(features::Int64, labeledData::Array{T, 2} where {T<:Number},
                       labeledDataLabels::Array{T} where {T<:Number}, K::Int64)
    classes = unique(labeledDataLabels)
    totalClasses = length(classes)

    # Finding the first centroids
    # If the number of classes is the same as K, the initial centroid of each class is the mean of each feature.
    # If it's not, a clustering is necessary to find the centroids.
    if K == totalClasses
        centroids = zeros(totalClasses, features + 1)

        for class = 1:totalClasses
            tempCentroids = zeros(features)

            for feature = 1:features
                tempCentroids[feature] = median(labeledData[findall(label -> label == classes[class],
                                                labeledDataLabels), feature])
            end

            centroids[class, 1:end-1] = tempCentroids
        end

        centroids[:, end] = classes
    else
        py"""
            from sklearn.cluster import KMeans
            import numpy as np
            import sys, os
            import warnings

            warnings.filterwarnings("ignore")

            def kmeans(X, k, init, centroids=[]):
                if (init == 0):
                    return KMeans(n_clusters=k).fit(X).cluster_centers_
                else:
                    centroids = np.array(centroids)
                    return KMeans(n_clusters=k, init=centroids).fit(X).cluster_centers_"""

        tempCentroids = py"kmeans"(labeledData, K, 0)

        centroidLabels = zeros(size(tempCentroids, 1))

        for row = 1:size(tempCentroids, 1)
            output, _ = knnClassification(labeledData, labeledDataLabels, tempCentroids[row, :])
            centroidLabels[row] = output
        end

        centroids = zeros(size(tempCentroids, 1), size(tempCentroids, 2) + 1)
        centroids[:, 1:end-1] = tempCentroids
        centroids[:, end] = centroidLabels
    end

    return centroids
end

"""
    findLabelForCurrentCentroids(
        pastCentroids    -> "centroids from previous iteration"
        currentCentroids -> "centroids from current iteration"
    )

The function uses the last iteration centroids' labels to define the current iteration centroids' labels.
Given the current centroids from the most recent unlabeled clusters and the past centroids from the previously labeled
clusters, each centroid from the past iteration have a label yi and each centroid from current iteration needs a label
ˆyi. This label is obtained by aplying nearest neighbor algorithm. In other words, the label given to a centroid in the
current iteration is the same as the label given to it's nearest neighbor in the past iteration.

The function returns two arrays:
    intermed       -> matrix storing the median between the nearest neighbor's data and the current centroids and it's labels
    centroidLabels -> labels got for the current centroids
"""
function findLabelForCurrentCentroids(pastCentroids::Array{T, 2} where {T<:Number}, currentCentroids::Array{T, 2} where {T<:Number})
    intermed = zeros(size(currentCentroids, 1), size(pastCentroids, 2))
    centroidLabels = zeros(size(currentCentroids, 1))

    for row = 1:size(currentCentroids, 1)
        label, nearestData = knnClassification(pastCentroids[:, 1:end-1], pastCentroids[:, end], currentCentroids[row, :])

        medianMatrix = zeros(2, size(currentCentroids, 2))
        medianMatrix[1, :] = nearestData
        medianMatrix[2, :] = currentCentroids[row, :]

        intermed[row, 1:end-1] = median(medianMatrix, dims=1)
        intermed[row, end] = label

        centroidLabels[row] = label
    end

    return intermed, centroidLabels
end

"""
    calculateNewLabeledData(
        poolData          -> "data stored in the pool"
        currentCentroids  -> "centroids from current iteration"
        pastCentroids     -> "matrix with centroids and their labels"
    )

Given the updated centroid values and the pool data, the function calculates the new labeled data using nearest neighbor
over the vcat of current centroids and the intermed centroids and the pool data.
The result is going to have the labels got with new centroid values, from the updated classifier.
"""
function calculateNewLabeledData(poolData::Array{T, 2} where {T<:Number}, currentCentroids::Array{T, 2} where {T<:Number},
                                 pastCentroids::Array{T, 2} where {T<:Number})
    newLabeledData = zeros(size(poolData, 1))

    lblData = vcat(currentCentroids[:, 1:end-1], pastCentroids[:, 1:end-1])
    lblLabels = vcat(currentCentroids[:, end], pastCentroids[:, end])

    for row = 1:size(poolData, 1)
        output, _ = knnClassification(lblData, lblLabels, poolData[row, 1:end-1])

        newLabeledData[row] = output
    end

    return newLabeledData
end

"""
    updateInformation(
        poolData             -> "test instances with it's labels"
        centroids            -> "centroids from past iteration"
        labeledData          -> "labeled data used in classification"
        labeledDataLabels    -> "labels from the data in `labeledData`"
        newLabeledData       -> "data calculated with the new centroid values"
        currentCentroids     -> "centroids from the current iteration"
        intermed             -> "matrix with all median data between past and current iteration centroids
        concordantLabelCount -> "amout of labels in concordance between the labels stored in pool data and the calculated
                                with new centroid values"
        maxPoolSize          -> "the maximum size of the pool"
    )

Updates the information stored in `centroids`, `labeledData` and `labeledDataLabels` if the concordance isn't 100% with the
new labeled data calculated with the updated centroid values.
"""
function updateInformation(poolData::Array{T, 2} where {T<:Number}, centroids::Array{T, 2} where {T<:Number}, labeledData::Array{T, 2}
                           where {T<:Number}, labeledDataLabels::Array{T} where {T<:Number}, newLabeledData::Array{T} where {T<:Number},
                           currentCentroids::Array{T, 2} where {T<:Number}, intermed::Array{T, 2} where {T<:Number},
                           concordantLabelCount::Int64, maxPoolSize::Int64)
    # If there's still some difference (concordantLabelCount/maxPoolSize < 1) means that the classifier elements
    # (as centroids, data and labels) need to be updated.
    if concordantLabelCount/maxPoolSize < 1
        centroids = vcat(currentCentroids, intermed)

        labeledData = poolData[:, 1:end-1]
        labeledDataLabels = newLabeledData
    end

    return centroids, labeledData, labeledDataLabels
end

end


text:
# I guarantee the brittleness of this macro, it tends to break my terminal and it is an abomination. Have fun!

"""
A more convenient syntax for try/catch clauses.
You know that you want it!
Instead of the boring old:
```julia
try
	foo()
catch e
	bar(e)
finally
	baz()
end
```
you can now write the far more legible:
```julia
💣 = foo
😥 = bar
🍌 = baz

@🐒 begin
    🙈
    	💣()
    🙊(💥)
    	😥(💥)
    🙉
    	🍌()
    end
```
"""
macro 🐒(monkeyexpression::Expr)
	monkeyexpression.head == :block || error("You have to wrap this in a begin...end block, sorry!")
	newexpr = Expr(:block)
	tryblock = Expr(:block)
	catchme = false
	catchblock = false
	finallyblock = nothing
	state = :start # where are we in the expression
	for sub in monkeyexpression.args
		if state == :start
			if sub isa Symbol && sub == :🙈
				state = :try
			elseif sub isa LineNumberNode
				push!(newexpr.args, sub)
			else
				error("Missing 🙈")
			end
		elseif state == :try
			if sub isa Symbol && sub == :🙊
				state = :catch
				catchblock = Expr(:block)
			elseif sub isa Expr && sub.args[1] == :🙊
				state = :catch
				catchblock = Expr(:block)
				if length(sub.args) == 2
					catchme = sub.args[2]
				else
					error("Can only catch a single error at once, duh!")
				end
			elseif sub isa Symbol && sub == :🙉
				state = :finally
				finallyblock = Expr(:block)
			else
				push!(tryblock.args, sub)
			end
		elseif state == :catch
			if sub isa Symbol && sub == :🙉
				state = :finally
				finallyblock = Expr(:block)
			else
				push!(catchblock.args, sub)
			end
		elseif state == :finally
			push!(finallyblock.args, sub)
		end
	end
	if state == :try
		error("Syntax: 🙈 without 🙊 or 🙉")
	end
	tryexpr = Expr(:try, tryblock, catchme, catchblock)
	if !(finallyblock === nothing)
		push!(tryexpr.args, finallyblock)
	end
	push!(newexpr.args, tryexpr)
	return esc(newexpr)
end

ext: jl


text:
module occa

#---[ Setup ]-----------------
macro libocca()
    return bytestring(ENV["OCCA_DIR"], "/lib/libocca.so")
end

#---[ Types ]-----------------
type device
    cDevice::Ptr{Void}
end

function device(infos::String)
    cDevice = ccall((:occaGetDevice, @libocca()),
                    Ptr{Void},
                    (Ptr{Uint8},),
                    bytestring(infos));

   return device(cDevice);
end

type stream
    cStream::Ptr{Void}

end

function stream(ptr::Ptr{Void})
    return stream(cStream);
end

type kernel
    cKernel::Ptr{Void}
end
function kernel(cKernel)
    return kernel(cKernel);
end

type kernelInfo
    cKernelInfo::Ptr{Void}
end

function kernelInfo(ptr::Ptr{Void})
    return kernlInfo(cKernelInfo);
end

type memory
    cMemory::Ptr{Void}
    cTypes
end

function memory(cMemory, cTypes)
    return memory(cMemory,cTypes);
end

#---[ Device ]----------------
function finalizer(d::device)
    ccall((:occaDeviceFree, @libocca()),
          Void,
          (Ptr{Void},),
          d.cDevice)
end

function mode(d::device)
    cMode = ccall((:occaDeviceMode, @libocca()),
                  Ptr{Uint8},
                  (Ptr{Void},), d.cDevice)

    return bytestring(cMode)
end

function setCompiler(d::device,
                     compiler::String)
    ccall((:occaDeviceSetCompiler, @libocca()),
          Void,
          (Ptr{Void}, Ptr{Uint8},),
          d.cDevice, bytestring(compiler))
end

function setCompilerFlags(d::device,
                          compilerFlags::String)
    ccall((:occaDeviceSetCompilerFlags, @libocca()),
          Void,
          (Ptr{Void}, Ptr{Uint8},),
          d.cDevice, bytestring(compilerFlags))
end

function buildKernelFromSource(d::device,
                               filename::String,
                               functionName::String,
                               info = C_NULL)
    if info == C_NULL
        cKernel = ccall((:occaBuildKernelFromSource, @libocca()),
                        Ptr{Void},
                        (Ptr{Void}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Void},),
                        d.cDevice,
                        bytestring(filename),
                        bytestring(functionName),
                        C_NULL)
    else
        cKernel = ccall((:occaBuildKernelFromSource, @libocca()),
                        Ptr{Void},
                        (Ptr{Void}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Void},),
                        d.cDevice,
                        bytestring(filename),
                        bytestring(functionName),
                        info.cKernelInfo)
    end

    return kernel(cKernel)
end

function buildKernelFromBinary(d::device,
                               filename::String,
                               functionName::String)
    cKernel = ccall((:occaBuildKernelFromBinary, @libocca()),
                    Ptr{Void},
                    (Ptr{Void}, Ptr{Uint8}, Ptr{Uint8},),
                    d.cDevice,
                    bytestring(filename),
                    bytestring(functionName))

    return kernel(cKernel)
end

function malloc(d::device, source::Array)
    cTypes = typeof(source[1])
    bytes  = length(source) * sizeof(cTypes)

    convert(Uint, bytes)

    cMemory = ccall((:occaDeviceMalloc, @libocca()),
                    Ptr{Void},
                    (Ptr{Void}, Uint, Ptr{Void},),
                    d.cDevice, bytes, pointer(source))

    return memory(cMemory, cTypes)
end

function malloc(d::device, entriesAndType)
    if length(entriesAndType) != 2
        error("malloc second argument must be a tuple of (bytes, type) or Array")
    end

    cTypes = entriesAndType[2]

    bytes  = entriesAndType[1] * sizeof(cTypes)

    convert(Uint, bytes)

    cMemory = ccall((:occaDeviceMalloc, @libocca()),
                    Ptr{Void},
                    (Ptr{Void}, Uint, Ptr{Void},),
                    d.cDevice, bytes, C_NULL)

    return memory(cMemory, cTypes)
end

function createStream(d::device)
    cStream = ccall((:occaGenStream, @libocca()),
                    Ptr{Void},
                    (Ptr{Void},),
                    d.cDevice)

    return stream(cStream)
end

function getStream(d::device)
    cStream = ccall((:occaGetStream, @libocca()),
                    Ptr{Void},
                    (Ptr{Void},),
                    d.cDevice)

    return stream(cStream)
end

function setStream(d::device, s::stream)
    ccall((:occaSetStream, @libocca()),
          Void,
          (Ptr{Void}, Ptr{Void},),
          d.cDevice, s.cStream)
end

#---[ Kernel ]----------------
function finalizer(k::kernel)
    ccall((:occaKernelFree, @libocca()),
          Void,
          (Ptr{Void},),
          k.cKernel)
end

function mode(k::kernel)
    cMode = ccall((:occaKernelMode, @libocca()),
                  Ptr{Uint8},
                  (Ptr{Void},),
                  k.cKernel)

    return bytestring(cMode)
end

function getPreferredDimSize(k::kernel)
    return ccall((:occaKernelPreferredDimSize, @libocca()),
                 Int32,
                 (Ptr{Void},),
                 k.cKernel)
end

function setWorkingDims(k::kernel,
                        dims, items, groups)
    convert(Int32, dims)

    items_  = ones(Uint, 3)
    groups_ = ones(Uint, 3)

    for i = 1:dims
        items_[i]  = items[i]
        groups_[i] = groups[i]
    end

    ccall((:occaKernelSetAllWorkingDims, @libocca()),
          Void,
          (Ptr{Void},
           Int32,
           Uint, Uint, Uint,
           Uint, Uint, Uint,),
          k.cKernel,
          dims,
          items_[1] , items_[2] , items_[3],
          groups_[1], groups_[2], groups_[3])
end

argType(arg::Int8)  = ccall((:occaChar , @libocca()), Ptr{Void}, (Int8,) , arg)
argType(arg::Uint8) = ccall((:occaUChar, @libocca()), Ptr{Void}, (Uint8,), arg)

argType(arg::Int16)  = ccall((:occaShort , @libocca()), Ptr{Void}, (Int16,) , arg)
argType(arg::Uint16) = ccall((:occaUShort, @libocca()), Ptr{Void}, (Uint16,), arg)

argType(arg::Int32)  = ccall((:occaInt , @libocca()), Ptr{Void}, (Int32,) , arg)
argType(arg::Uint32) = ccall((:occaUInt, @libocca()), Ptr{Void}, (Uint32,), arg)

argType(arg::Int64)  = ccall((:occaLong , @libocca()), Ptr{Void}, (Int64,) , arg)
argType(arg::Uint64) = ccall((:occaULong, @libocca()), Ptr{Void}, (Uint64,), arg)

argType(arg::Float32) = ccall((:occaFloat , @libocca()), Ptr{Void}, (Float32,) , arg)
argType(arg::Float64) = ccall((:occaDouble, @libocca()), Ptr{Void}, (Float64,) , arg)

function runKernel(k::kernel, args...)
    argList = ccall((:occaGenArgumentList, @libocca()),
                    Ptr{Void}, ())

    pos = convert(Int32, 0)

    for arg in args
        if isa(arg, memory)
            ccall((:occaArgumentListAddArg, @libocca()),
                  Void,
                  (Ptr{Void}, Int32, Ptr{Void},),
                  argList, pos, arg.cMemory)
        else
            if length(arg) != 2
                error("Kernel argument should be in the form of (value, type)")
            end

            arg_ = arg[1]
            convert(arg[2], arg_)

            cArg = argType(arg_)

            ccall((:occaArgumentListAddArg, @libocca()),
                  Void,
                  (Ptr{Void}, Int32, Ptr{Void},),
                  argList, pos, cArg)
        end

        pos += 1
    end

    ccall((:occaKernelRun_, @libocca()),
          Void,
          (Ptr{Void}, Ptr{Void},),
          k.cKernel, argList)

    ccall((:occaArgumentListFree, @libocca()),
          Void,
          (Ptr{Void},),
          argList)
end

function timeTaken(k::kernel)
    return ccall((:occaKernelTimeTaken, @libocca()),
                 Float64,
                 (Ptr{Void},),
                 k.cKernel)
end

function addDefine(info::kernelInfo, macro_::String, value::String)
    occaValue = ccall((:occaString, @libocca()),
                      Ptr{Void},
                      (Ptr{Uint8},),
                      bytestring(value))

    ccall((:occaKernelInfoAddDefine, @libocca()),
          Void,
          (Ptr{Void}, Ptr{Uint8}, Ptr{Void},),
          info.cKernelInfo, bytestring(macro_), occaValue)
end

function finalizer(info::kernelInfo)
    ccall((:occaKernelInfoFree, @libocca()),
          Void,
          (Ptr{Void},),
          info.cKernelInfo)
end

#---[ Memory ]----------------
function finalizer(m::memory)
    ccall((:occaMemoryFree, @libocca()),
          Void,
          (Ptr{Void},),
          m.cMemory)
end

function mode(m::memory)
    cMode = ccall((:occaMemoryMode, @libocca()),
                  Ptr{Uint8},
                  (Ptr{Void},),
                  m.cMemory)

    return bytestring(cMode)
end

function memcpy(destTuple, srcTuple, bytes::Number = 0)
    if isa(destTuple, memory)
        dest = destTuple.cMemory

        destOffset = 0
        convert(Uint, destOffset)

        destIsAMemory = true
    elseif isa(destTuple, Array)
        dest = pointer(destTuple)

        destOffset = 0
        convert(Uint, destOffset)

        destIsAMemory = false
    else
        dest = destTuple[1]

        if isa(dest, memory)
            dest = dest.cMemory
            destIsAMemory = true
        else
            dest = pointer(destTuple[1])
            destIsAMemory = false
        end

        destOffset = destTuple[2]
        convert(Uint, destOffset)
    end

    if isa(srcTuple, memory)
        src = srcTuple.cMemory

        srcOffset = 0
        convert(Uint, srcOffset)

        srcIsAMemory = true
    elseif isa(srcTuple, Array)
        src = pointer(srcTuple)

        srcOffset = 0
        convert(Uint, srcOffset)

        srcIsAMemory = false
    else
        src = srcTuple[1]

        if isa(src, memory)
            src = src.cMemory
            srcIsAMemory = true
        else
            src = pointer(srcTuple[1])
            srcIsAMemory = false
        end

        srcOffset = srcTuple[2]
        convert(Uint, srcOffset)
    end

    convert(Uint, bytes)

    if destIsAMemory
        if srcIsAMemory
            ccall((:occaCopyMemToMem, @libocca()),
                  Void,
                  (Ptr{Void}, Ptr{Void}, Uint, Uint, Uint,),
                  dest, src, bytes, destOffset, srcOffset)
        else
            ccall((:occaCopyPtrToMem, @libocca()),
                  Void,
                  (Ptr{Void}, Ptr{Void}, Uint, Uint,),
                  dest, src, bytes, destOffset)
        end
    else
        if srcIsAMemory
            ccall((:occaCopyMemToPtr, @libocca()),
                  Void,
                  (Ptr{Void}, Ptr{Void}, Uint, Uint,),
                  dest, src, bytes, srcOffset)
        else
            error("One of the arguments should be an OCCA memory type")
        end
    end
end

function swap(a::memory, b::memory)
    tmp       = a.cMemory
    a.cMemory = b.cMemory
    b.cMemory = tmp
end

end


text:
using Test
using GraphMakie
using GeometryBasics

using GraphMakie: BezierPath, MoveTo, LineTo, CurveTo, interpolate, discretize, tangent, waypoints, Path, Line

@testset "interpolation and tangents" begin
    path = BezierPath([
        MoveTo(Point2f(0,0)),
        LineTo(Point2f(0.5,0.5)),
        CurveTo(Point2f(1,0),
                Point2f(1,0),
                Point2f(1,1)),
        CurveTo(Point2f(1,2),
                Point2f(0,2),
                Point2f(0,1))])

    lines(discretize(path))

    ts = 0:0.1:1.0
    pos = map(t->interpolate(path,t), ts)
    tan = map(t->tangent(path,t), ts)

    scatter!(pos)
    arrows!(pos, tan; lengthscale=0.1)

    path = BezierPath([
        MoveTo(Point3f(0,0,0)),
        LineTo(Point3f(1,0,0)),
        LineTo(Point3f(1,1,1)),
        LineTo(Point3f(0,1,1)),
        LineTo(Point3f(0,0,1))])
    ts = 0:0.1:1.0
    pos = map(t->interpolate(path,t), ts)
    tan = map(t->tangent(path,t), ts)
    fig, ax, p = scatter(pos)
    arrows!(pos, tan; lengthscale=0.1)
end

@testset "natural spline constructor" begin
    path = Path(Point2f(0.0,0.0), Point2f(0.0,1.0))
    @test path isa GraphMakie.Line

    path = Path(Point2f(0.0,0.0),
                Point2f(0.0,1.0),
                Point2f(1.0,1.0))
    lines(discretize(path))
    scatter!(waypoints(path))

    path = Path(Point2f(0.0,0.0),
                Point2f(-0.5,1.0),
                Point2f(0.5,1.0),
                Point2f(0.0,0.0))
    lines(discretize(path))
    scatter!(waypoints(path))

    path = Path(Point2f(0.0,0.0),
                Point2f(-0.5,1.0),
                Point2f(1.5,1.0),
                Point2f(2.0,0.0))
    lines(discretize(path))
    scatter!(waypoints(path))

    path = Path(Point2f(0.0,0.0),
                Point2f(-0.5,1.0),
                Point2f(1.5,1.0),
                Point2f(2.0,0.0);
                tangents=(Point2f(-1,0),
                          Point2f(-1,0)))
    lines(discretize(path), linewidth=10)
    scatter!(waypoints(path))
end

@testset "two points and tangents" begin
    p1 = Point2f(0,0)
    t1 = Point2f(0,1)
    p2 = Point2f(1,1)
    t2 = Point2f(0,1)
    path = Path(p1, p2; tangents=(t1, t2))
    lines(discretize(path))
    scatter!(waypoints(path))

    p1 = Point2f(0,0)
    t1 = Point2f(1,0)*10
    p2 = Point2f(1,1)
    t2 = Point2f(0,1)*5
    path = Path(p1, p2; tangents=(t1, t2))
    lines(discretize(path))
    scatter!(waypoints(path))

    p1 = Point2f(0,0)
    t1 = Point2f(-1,0)
    p2 = Point2f(1,1)
    t2 = Point2f(0,1)
    path = Path(p1, p2; tangents=(t1, t2))
    lines(discretize(path))
    scatter!(waypoints(path))

    p1 = Point2f(0,0)
    t1 = Point2f(-.5,.5)
    p2 = Point2f(0,0)
    t2 = Point2f(-.5,-.5)
    path = Path(p1, p2; tangents=(t1, t2))
    lines(discretize(path))
    scatter!(waypoints(path))

    p1 = Point2f(0,0)
    p2 = Point2f(1,0)
    midp = Point2f[]
    path = Path(p1, midp..., p2)
    lines(discretize(path))
end

@testset "straight lines with radi" begin
    using GraphMakie: plot_controlpoints!
    p1 = Point2f(0,0)
    p2 = Point2f(1,-.5)
    p3 = Point2f(2,.5)
    p4 = Point2f(3,0)
    path = Path(0.5, p1, p2, p3, p4)
    fig, ax, p = lines(discretize(path))
    plot_controlpoints!(ax, path)
    scatter!(ax, [p1,p2,p3,p4])

    path = Path(0.0, p1, p2, p3, p4)
    fig, ax, p = lines(discretize(path))
end

@testset "test beziersegments recipe" begin
    using GraphMakie: beziersegments
    paths = [Path(rand(Point2f, 4)...) for _ in 1:4]
    fig, ax, p = beziersegments(paths; linewidth=[2,4,6,8], color=[1,2,3,4])
    p.attributes
end

@testset "selfloop test" begin
    using Graphs, GraphMakie
    g = star_graph(10)
    add_edge!(g, 1, 1)
    add_edge!(g, 2, 2)
    fig, ax, p = graphplot(g)
    ax.aspect = DataAspect()

    g = star_graph(10)
    add_edge!(g, 1, 1)
    @test_throws ErrorException graphplot(g; layout=_->rand(Point3f, 10))
end


text:
group = get(ARGS,1,"default")

if group == "default"
    include("gmres.jl")
    include("gcr.jl")
    include("cg.jl")
elseif group == "gpu"
    include("gmres_gpu.jl")
end

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
using QuantileRegressions
using DataFrames, CSV
cd(joinpath(dirname(pathof(QuantileRegressions)), "..", "examples"))
df = CSV.read("engel.csv")

out_ip = qreg(foodexp~income, df, IP()) # or just qreg(foodexp~income, df)
out_irls = qreg(foodexp~income, df, IRLS())

println("...... Interior point ......")
println(out_ip)
println(".......... IRLS ............")
println(out_irls)


text:
Base.@deprecate_binding ChannelView channelview

export ColorView

struct ColorView{C<:Colorant,N,A<:AbstractArray} <: AbstractArray{C,N}
    parent::A

    function ColorView{C,N,A}(parent::AbstractArray{T}) where {C,N,A,T<:Number}
        n = length(colorview_size(C, parent))
        n == N || throw(DimensionMismatch("for an $N-dimensional ColorView with color type $C, input dimensionality should be $n instead of $(ndims(parent))"))
        checkdim1(C, axes(parent))
        Base.depwarn("ColorView{C}(A) is deprecated, use colorview(C, A)", :ColorView)
        colorview(C, A)
    end
end

function ColorView{C}(A::AbstractArray) where C<:Colorant
    Base.depwarn("ColorView{C}(A) is deprecated, use colorview(C, A)", :ColorView)
    colorview(C, A)
end

ColorView(parent::AbstractArray) = error("must specify the colortype, use colorview(C, A)")

Base.@deprecate_binding squeeze1 true


text:
using Dates, Formatting, Plots, ColorSchemes, DataFrames

include("load_covid19_data.jl")

states_of_interest = ["Virginia","North Carolina","West Virginia","Delaware", "Wisconsin", "Minnesota", "Idaho", "Tennessee", "Alabama", "New York", "New Jersey", "Massachusetts", "Texas","Florida","California","Michigan", "Ohio", "Washington", "Oregon", "Illinois", "Oklahoma", "Maryland", "District of Columbia", "Alaska", "Arizona","Georgia","South Carolina", "Mississippi", "Maine", "Pennsylvania", "Puerto Rico", "Colorado", "New Hampshire", "Iowa", "Vermont","Hawaii"]

strnow = string(Dates.now())
strnow30 = strnow[1:4] * strnow[6:7] * strnow[9:10] * "T" * strnow[12:13] * strnow[15:16] * strnow[18:19]

# extracting county indices with land area and population
cind = findall((countyarea .> 0) .& (cpop .> 0));

reasonable_resolution() = (1000, 800)

tind = 1:length(covid19us[1].confirmed);
t = covid19us[1].time[tind];


text:
"""
Only packs `A`. Primitively does column-major packing: it packs blocks of `A` into a column-major temporary.
"""
function matmul_st_only_pack_A!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, M, K, N, ::StaticFloat64{W₁}, ::StaticFloat64{W₂}, ::StaticFloat64{R₁}, ::StaticFloat64{R₂}
) where {T, W₁, W₂, R₁, R₂}

    mᵣ, nᵣ = matmul_params()
    ((Mblock, Mblock_Mrem, Mremfinal, Mrem, Miter), (Kblock, Kblock_Krem, Krem, Kiter)) =
        solve_McKc(T, M, K, N, StaticFloat64{W₁}(), StaticFloat64{W₂}(), StaticFloat64{R₁}(), StaticFloat64{R₂}(), mᵣ)
    for ko ∈ CloseOpen(Kiter)
        ksize = ifelse(ko < Krem, Kblock_Krem, Kblock)
        let A = A, C = C
            for mo in CloseOpen(Miter)
                msize = ifelse((mo+1) == Miter, Mremfinal, ifelse(mo < Mrem, Mblock_Mrem, Mblock))
                # if ko == 0
                #     loopmul!(C, A, B, α, β, msize, ksize, N)
                # else
                #     loopmul!(C, A, B, α, One(), msize, ksize, N)
                # end
                if ko == 0
                    packaloopmul!(C, A, B, α, β, msize, ksize, N)
                else
                    packaloopmul!(C, A, B, α, One(), msize, ksize, N)
                end
                A = gesp(A, (msize, Zero()))
                C = gesp(C, (msize, Zero()))
            end
        end
        A = gesp(A, (Zero(), ksize))
        B = gesp(B, (ksize, Zero()))
    end
    nothing
end

function matmul_st_pack_A_and_B!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer, α, β, M, K, N, W₁, W₂, R₁, R₂, tid
) where {T}
    mᵣ, nᵣ = matmul_params()
    # TODO: if this is nested in other threaded code, use only a piece of BCACHE and make R₂ (and thus L₂ₑ) smaller
    (Mblock, Mblock_Mrem, Mremfinal, Mrem, Miter), (Kblock, Kblock_Krem, Krem, Kiter), (Nblock, Nblock_Nrem, Nrem, Niter) =
        solve_block_sizes(T, M, K, N, W₁, W₂, R₁, R₂, mᵣ)

    bcache = _use_bcache(tid)
    L3ptr = Base.unsafe_convert(Ptr{T}, pointer(bcache))
    for n ∈ CloseOpen(Niter)
        nsize = ifelse(n < Nrem, Nblock_Nrem, Nblock)
        let A = A, B = B
            for k ∈ CloseOpen(Kiter)
                ksize = ifelse(k < Krem, Kblock_Krem, Kblock)
                _B = default_zerobased_stridedpointer(L3ptr, (One(),ksize))
                unsafe_copyto_avx!(_B, B, ksize, nsize)
                let A = A, C = C, B = _B
                    for m in CloseOpen(Miter)
                        msize = ifelse((m+1) == Miter, Mremfinal, ifelse(m < Mrem, Mblock_Mrem, Mblock))
                        if k == 0
                            packaloopmul!(C, A, B, α,     β, msize, ksize, nsize)
                        else
                            packaloopmul!(C, A, B, α, One(), msize, ksize, nsize)
                        end
                        A = gesp(A, (msize, Zero()))
                        C = gesp(C, (msize, Zero()))
                    end
                end
                A = gesp(A, (Zero(), ksize))
                B = gesp(B, (ksize, Zero()))
            end
        end
        B = gesp(B, (Zero(), nsize))
        C = gesp(C, (Zero(), nsize))
    end
    _free_bcache!(bcache)
    nothing
end

@inline contiguousstride1(A) = ArrayInterface.contiguous_axis(A) === One()
@inline contiguousstride1(A::AbstractStridedPointer{T,N,1}) where {T,N} = true
# @inline bytestride(A::AbstractArray, i) = VectorizationBase.bytestrides(A)[i]
@inline bytestride(A::AbstractStridedPointer, i) = strides(A)[i]
@inline firstbytestride(A::AbstractStridedPointer) = bytestride(A, One())

@inline function vectormultiple(bytex, ::Type{Tc}, ::Type{Ta}) where {Tc,Ta}
    Wc = pick_vector_width(Tc) * static_sizeof(Ta) - One()
    iszero(bytex & (VectorizationBase.register_size() - One()))
end
@inline function dontpack(pA::AbstractStridedPointer{Ta}, M, K, ::StaticInt{mc}, ::StaticInt{kc}, ::Type{Tc}) where {mc, kc, Tc, Ta}
    (contiguousstride1(pA) &&
         ((((MᵣW_mul_factor() + StaticInt(5)) * pick_vector_width(Tc)) ≥ M) ||
          (vectormultiple(bytestride(pA, StaticInt{2}()), Tc, Ta) && ((M * K) ≤ (mc * kc)) && iszero(reinterpret(Int, pointer(pA)) & (VectorizationBase.register_size() - One())))))
end

@inline function alloc_matmul_product(A::AbstractArray{TA}, B::AbstractArray{TB}) where {TA,TB}
    # TODO: if `M` and `N` are statically sized, shouldn't return a `Matrix`.
    M, KA = size(A)
    KB, N = size(B)
    @assert KA == KB "Size mismatch."
    Matrix{promote_type(TA,TB)}(undef, M, N), (M, KA, N)
end
@inline function matmul_serial(A::AbstractMatrix, B::AbstractMatrix)
    C, (M,K,N) = alloc_matmul_product(A, B)
    _matmul_serial!(C, A, B, One(), Zero(), (M,K,N))
    return C
end

# These methods must be compile time constant
maybeinline(::Any, ::Any, ::Any, ::Any) = false
function maybeinline(::StaticInt{M}, ::StaticInt{N}, ::Type{T}, ::Val{true}) where {M,N,T}
    mᵣ, nᵣ = matmul_params()
    static_sizeof(T) * StaticInt{M}() * StaticInt{N}() < StaticInt{176}() * mᵣ * nᵣ
end
function maybeinline(::StaticInt{M}, ::StaticInt{N}, ::Type{T}, ::Val{false}) where {M,N,T}
    StaticInt{M}() * static_sizeof(T) ≤ StaticInt{2}() * VectorizationBase.register_size()
end

@inline function matmul_serial!(C::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix) where {T}
    matmul_serial!(C, A, B, One(), Zero(), nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α)
    matmul_serial!(C, A, B, α, Zero(), nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β)
    matmul_serial!(C, A, B, α, β, nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, MKN, ::StaticInt{2})
    _matmul_serial!(C', B', A', α, β, nothing)
    return C
end
@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, MKN, ::StaticInt)
    _matmul_serial!(C, A, B, α, β, nothing)
    return C
end

"""
  matmul_serial!(C, A, B[, α = 1, β = 0])

Calculates `C = α * (A * B) + β * C` in place.

A single threaded matrix-matrix-multiply implementation.
Supports dynamically and statically sized arrays.

Organizationally, `matmul_serial!` checks the arrays properties to try and dispatch to an appropriate implementation.
If the arrays are small and statically sized, it will dispatch to an inlined multiply.

Otherwise, based on the array's size, whether they are transposed, and whether the columns are already aligned, it decides to not pack at all, to pack only `A`, or to pack both arrays `A` and `B`.
"""
@inline function _matmul_serial!(
    C::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix, α, β, MKN
) where {T}
    M, K, N = MKN === nothing ? matmul_sizes(C, A, B) : MKN
    if M * N == 0
        return
    elseif K == 0
        matmul_only_β!(C, β)
        return
    end
    pA = zstridedpointer(A); pB = zstridedpointer(B); pC = zstridedpointer(C);
    Cb = preserve_buffer(C); Ab = preserve_buffer(A); Bb = preserve_buffer(B);
    Mc, Kc, Nc = block_sizes(T); mᵣ, nᵣ = matmul_params();
    GC.@preserve Cb Ab Bb begin
        if maybeinline(M, N, T, ArrayInterface.is_column_major(A)) # check MUST be compile-time resolvable
            inlineloopmul!(pC, pA, pB, One(), Zero(), M, K, N)
            return
        elseif (nᵣ ≥ N) || dontpack(pA, M, K, Mc, Kc, T)
            loopmul!(pC, pA, pB, α, β, M, K, N)
            return
        else
            matmul_st_pack_dispatcher!(pC, pA, pB, α, β, M, K, N)
            return
        end
    end
end # function

function matmul_only_β!(C::AbstractMatrix{T}, β::StaticInt{0}) where T
    @avx for i=1:length(C)
        C[i] = zero(T)
    end
end

function matmul_only_β!(C::AbstractMatrix{T}, β) where T
    @avx for i=1:length(C)
        C[i] = β * C[i]
    end
end

function matmul_st_pack_dispatcher!(pC::AbstractStridedPointer{T}, pA, pB, α, β, M, K, N) where {T}
    Mc, Kc, Nc = block_sizes(T)
    if (contiguousstride1(pB) ? (Kc * Nc ≥ K * N) : (firstbytestride(pB) ≤ 1600))
        matmul_st_only_pack_A!(pC, pA, pB, α, β, M, K, N, W₁Default(), W₂Default(), R₁Default(), R₂Default())
    # elseif notnested !== nothing && notnested
    #     matmul_st_pack_A_and_B!(pC, pA, pB, α, β, M, K, N, W₁Default(), W₂Default(), R₁Default(), R₂Default(), nothing)
    else
        matmul_st_pack_A_and_B!(pC, pA, pB, α, β, M, K, N, W₁Default(), W₂Default(), R₁Default(), R₂Default()/Threads.nthreads(), Threads.threadid() - 1)
    end
    nothing
end

"""
    matmul(A, B)

Multiply matrices `A` and `B`.
"""
@inline function matmul(A::AbstractMatrix, B::AbstractMatrix)
    C, (M,K,N) = alloc_matmul_product(A, B)
    _matmul!(C, A, B, One(), Zero(), nothing, (M,K,N))
    return C
end

"""
    matmul!(C, A, B[, α, β, max_threads])

Calculates `C = α * A * B + β * C` in place, overwriting the contents of `A`.
It may use up to `max_threads` threads. It will not use threads when nested in other threaded code.
"""
@inline function matmul!(C::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix) where {T}
    matmul!(C, A, B, One(), Zero(), nothing, nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α)
    matmul!(C, A, B, α, Zero(), nothing, nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β)
    matmul!(C, A, B, α, β, nothing, nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, nthread)
    matmul!(C, A, B, α, β, nthread, nothing, ArrayInterface.contiguous_axis(C))
end
@inline function matmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, nthread, MKN, ::StaticInt{2})
    _matmul!(C', B', A', α, β, nthread, MKN)
    return C
end
@inline function matmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, nthread, MKN, ::StaticInt)
    _matmul!(C, A, B, α, β, nthread, MKN)
    return C
end

@inline function dontpack(pA::AbstractStridedPointer{Ta}, M, K, ::StaticInt{mc}, ::StaticInt{kc}, ::Type{Tc}, nspawn) where {mc, kc, Tc, Ta}
    # TODO: perhaps consider K vs kc by themselves?
    (contiguousstride1(pA) && ((M * K) ≤ (mc * kc) * nspawn >>> 1))
end

# passing MKN directly would let osmeone skip the size check.
@inline function _matmul!(C::AbstractMatrix{T}, A, B, α, β, nthread, MKN) where {T}#::Union{Nothing,Tuple{Vararg{Integer,3}}}) where {T}
    M, K, N = MKN === nothing ? matmul_sizes(C, A, B) : MKN
    if M * N == 0
        return
    elseif K == 0
        matmul_only_β!(C, β)
        return
    end
    W = pick_vector_width(T)
    pA = zstridedpointer(A); pB = zstridedpointer(B); pC = zstridedpointer(C);
    Cb = preserve_buffer(C); Ab = preserve_buffer(A); Bb = preserve_buffer(B);
    mᵣ, nᵣ = matmul_params()
    GC.@preserve Cb Ab Bb begin
        if maybeinline(M, N, T, ArrayInterface.is_column_major(A)) # check MUST be compile-time resolvable
            inlineloopmul!(pC, pA, pB, One(), Zero(), M, K, N)
            return
        else
            (nᵣ ≥ N) && @goto LOOPMUL
            if (Sys.ARCH === :x86_64) || (Sys.ARCH === :i686)
                (M*K*N < (StaticInt{4_096}() * W)) && @goto LOOPMUL
            else
                (M*K*N < (StaticInt{32_000}() * W)) && @goto LOOPMUL
            end
            __matmul!(pC, pA, pB, α, β, M, K, N, nthread)
            return
            @label LOOPMUL
            loopmul!(pC, pA, pB, α, β, M, K, N)
            return
        end
    end
end

# This funciton is sort of a `pun`. It splits aggressively (it does a lot of "splitin'"), which often means it will split-N.
function matmulsplitn!(C::AbstractStridedPointer{T}, A, B, α, β, ::StaticInt{Mc}, M, K, N, nspawn, ::Val{PACK}) where {T, Mc, PACK}
    Mᵣ, Nᵣ = matmul_params()
    W = pick_vector_width(T)
    MᵣW = Mᵣ*W
    _Mblocks, Nblocks = divide_blocks(M, cld_fast(N, Nᵣ), nspawn, W)
    Mbsize, Mrem, Mremfinal, Mblocks = split_m(M, _Mblocks, W)
    # Nblocks = min(N, _Nblocks)
    Nbsize, Nrem = divrem_fast(N, Nblocks)

    _nspawn = Mblocks * Nblocks
    Mbsize_Mrem, Mbsize_ = promote(Mbsize +     W, Mbsize)
    Nbsize_Nrem, Nbsize_ = promote(Nbsize + One(), Nbsize)

    let _A = A, _B = B, _C = C, n = 0, tnum = 0, Nrc = Nblocks - Nrem, Mrc = Mblocks - Mrem, __Mblocks = Mblocks - One()
        while true
            nsize = ifelse(Nblocks > Nrc, Nbsize_Nrem, Nbsize_); Nblocks -= 1
            let _A = _A, _C = _C, __Mblocks = __Mblocks
                while __Mblocks != 0
                    msize = ifelse(__Mblocks ≥ Mrc, Mbsize_Mrem, Mbsize_); __Mblocks -= 1
                    launch_thread_mul!(_C, _A, _B, α, β, msize, K, nsize, (tnum += 1), Val{PACK}())
                    _A = gesp(_A, (msize, Zero()))
                    _C = gesp(_C, (msize, Zero()))
                end
                if Nblocks != 0
                    launch_thread_mul!(_C, _A, _B, α, β, Mremfinal, K, nsize, (tnum += 1), Val{PACK}())
                else
                    call_loopmul!(_C, _A, _B, α, β, Mremfinal, K, nsize, Val{PACK}())
                    waitonmultasks(CloseOpen(One(), _nspawn))
                    return
                end
            end
            _B = gesp(_B, (Zero(), nsize))
            _C = gesp(_C, (Zero(), nsize))
        end
    end
end

function __matmul!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer, α, β, M, K, N, nthread
) where {T}
    Mᵣ, Nᵣ = matmul_params()
    W = pick_vector_width(T)
    Mc, Kc, Nc = block_sizes(T)
    MᵣW = Mᵣ*W

    # Not taking the fast path
    # But maybe we don't want to thread anyway
    # Maybe this is nested, or we have ≤ 1 threads
    nt = _nthreads()
    _nthread = nthread === nothing ? nt : min(nt, nthread)
    if _nthread < 2
        matmul_st_pack_dispatcher!(C, A, B, α, β, M, K, N)
        return
    end
    # We are threading, but how many threads?
    nspawn = if (Sys.ARCH === :x86_64) || (Sys.ARCH === :i686)
        clamp(div_fast(M * N, StaticInt{128}() * W), 1, _nthread)
    else
        clamp(div_fast(M * N, StaticInt{256}() * W), 1, _nthread)
    end
    # nkern = cld_fast(M * N,  MᵣW * Nᵣ)

    # Approach:
    # Check if we don't want to pack A,
    #    if not, aggressively subdivide
    # if so, check if we don't want to pack B
    #    if not, check if we want to thread `N` loop anyway
    #       if so, divide `M` first, then use ratio of desired divisions / divisions along `M` to calc divisions along `N`
    #       if not, only thread along `M`. These don't need syncing, as we're not packing `B`
    #    if so, `matmul_pack_A_and_B!`
    #
    # MᵣW * (MᵣW_mul_factor - One()) # gives a smaller Mc, then
    # if 2M/nspawn is less than it, we don't don't `A`
    # First check is: do we just want to split aggressively?
    mᵣ, nᵣ = matmul_params()
    if dontpack(A, M, K, Mc, Kc, T, nspawn) || (W ≥ M) || (nᵣ*((num_cores() ≥ StaticInt(8)) ? max(nspawn,8) : 8) ≥ N)
        # `nᵣ*nspawn ≥ N` is needed at the moment to avoid accidentally splitting `N` to be `< nᵣ` while packing
        # Should probably handle that with a smarter splitting function...
        matmulsplitn!(C, A, B, α, β, Mc, M, K, N, nspawn, Val{false}())
    elseif ((nspawn*(W+W) > M) || (contiguousstride1(B) ? (roundtostaticint(Kc * Nc * R₂Default()) ≥ K * N) : (firstbytestride(B) ≤ 1600)))
        matmulsplitn!(C, A, B, α, β, Mc, M, K, N, nspawn, Val{true}())
    else # TODO: Allow splitting along `N` for `matmul_pack_A_and_B!`
        matmul_pack_A_and_B!(C, A, B, α, β, M, K, N, nspawn, W₁Default(), W₂Default(), R₁Default(), R₂Default())
    end
    nothing
end

# If tasks is [0,1,2,3] (e.g., `CloseOpen(0,4)`), it will wait on `MULTASKS[i]` for `i = [1,2,3]`.
function waitonmultasks(tasks)
    for tid ∈ tasks
        wait(tid)
    end
end

@inline allocref(::StaticInt{N}) where {N} = Ref{NTuple{N,UInt8}}()
function matmul_pack_A_and_B!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer, α, β, M, K, N,
    tospawn::Int, ::StaticFloat64{W₁}, ::StaticFloat64{W₂}, ::StaticFloat64{R₁}, ::StaticFloat64{R₂}#, ::Val{1}
) where {T,W₁,W₂,R₁,R₂}
    W = pick_vector_width(T)
    mᵣ, nᵣ = matmul_params()
    mᵣW = mᵣ * W
    # atomicsync = Ref{NTuple{16,UInt}}()
    Mbsize, Mrem, Mremfinal, _to_spawn = split_m(M, tospawn, W) # M is guaranteed to be > W because of `W ≥ M` condition for `jmultsplitn!`...
    atomicsync = allocref(StaticInt{2}()*num_cores()*cache_linesize())
    p = reinterpret(Ptr{UInt}, Base.unsafe_convert(Ptr{UInt8}, atomicsync))
    GC.@preserve atomicsync begin
        for i ∈ CloseOpen(2_to_spawn)
            _atomic_store!(p + i*cache_linesize(), zero(UInt))
        end
        Mblock_Mrem, Mblock_ = promote(Mbsize + W, Mbsize)
        u_to_spawn = _to_spawn % UInt
        tid = 0
        bc = _use_bcache()
        bc_ptr = Base.unsafe_convert(typeof(pointer(C)), pointer(bc))
        last_id = _to_spawn - One()
        for m ∈ CloseOpen(last_id) # ...thus the fact that `CloseOpen()` iterates at least once is okay.
            Mblock = ifelse(m < Mrem, Mblock_Mrem, Mblock_)
            launch_thread_mul!(C, A, B, α, β, Mblock, K, N, p, bc_ptr, m % UInt, u_to_spawn, StaticFloat64{W₁}(),StaticFloat64{W₂}(),StaticFloat64{R₁}(),StaticFloat64{R₂}())
            A = gesp(A, (Mblock, Zero()))
            C = gesp(C, (Mblock, Zero()))
        end
        sync_mul!(C, A, B, α, β, Mremfinal, K, N, p, bc_ptr, last_id % UInt, u_to_spawn, StaticFloat64{W₁}(), StaticFloat64{W₂}(), StaticFloat64{R₁}(), StaticFloat64{R₂}())
        waitonmultasks(CloseOpen(One(), _to_spawn))
    end
    _free_bcache!(bc)
    return
end

function sync_mul!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer, α, β, M, K, N, atomicp::Ptr{UInt}, bc::Ptr, id::UInt, total_ids::UInt,
    ::StaticFloat64{W₁}, ::StaticFloat64{W₂}, ::StaticFloat64{R₁}, ::StaticFloat64{R₂}
) where {T, W₁, W₂, R₁, R₂}

    (Mblock, Mblock_Mrem, Mremfinal, Mrem, Miter), (Kblock, Kblock_Krem, Krem, Kiter), (Nblock, Nblock_Nrem, Nrem, Niter) =
        solve_block_sizes(T, M, K, N, StaticFloat64{W₁}(), StaticFloat64{W₂}(), StaticFloat64{R₁}(), StaticFloat64{R₂}(), One())

    # atomics = atomicp + 8sizeof(UInt)
    sync_iters = zero(UInt)
    myp = atomicp +   id *cache_linesize()
    atomicp -= cache_linesize()
    atomics = atomicp + total_ids*cache_linesize()
    mys = myp + total_ids*(cache_linesize() % UInt)
    Npackb_r_div, Npackb_r_rem = divrem_fast(Nblock_Nrem, total_ids)
    Npackb_r_block_rem, Npackb_r_block_ = promote(Npackb_r_div + One(), Npackb_r_div)

    Npackb___div, Npackb___rem = divrem_fast(Nblock, total_ids)
    Npackb___block_rem, Npackb___block_ = promote(Npackb___div + One(), Npackb___div)

    pack_r_offset = Npackb_r_div * id + min(id, Npackb_r_rem)
    pack___offset = Npackb___div * id + min(id, Npackb___rem)

    pack_r_len = ifelse(id < Npackb_r_rem, Npackb_r_block_rem, Npackb_r_block_)
    pack___len = ifelse(id < Npackb___rem, Npackb___block_rem, Npackb___block_)

    for n in CloseOpen(Niter)
        # Krem
        # pack kc x nc block of B
        nfull = n < Nrem
        nsize = ifelse(nfull, Nblock_Nrem, Nblock)
        pack_offset = ifelse(nfull, pack_r_offset, pack___offset)
        pack_len = ifelse(nfull, pack_r_len, pack___len)
        let A = A, B = B
            for k ∈ CloseOpen(Kiter)
                ksize = ifelse(k < Krem, Kblock_Krem, Kblock)
                _B = default_zerobased_stridedpointer(bc, (One(), ksize))
                unsafe_copyto_avx!(gesp(_B, (Zero(), pack_offset)), gesp(B, (Zero(), pack_offset)), ksize, pack_len)
                # synchronize before starting the multiplication, to ensure `B` is packed
                _mv = _atomic_add!(myp, one(UInt))
                sync_iters += one(UInt)
                let atomp = atomicp
                    for _ ∈ CloseOpen(total_ids)
                        atomp += cache_linesize()
                        atomp == myp && continue
                        while _atomic_load(atomp) != sync_iters
                            pause()
                        end
                    end
                end
                # multiply
                let A = A, B = _B, C = C
                    for m in CloseOpen(Miter)
                        msize = ifelse((m+1) == Miter, Mremfinal, ifelse(m < Mrem, Mblock_Mrem, Mblock))
                        if k == 0
                            packaloopmul!(C, A, B, α,     β, msize, ksize, nsize)
                        else
                            packaloopmul!(C, A, B, α, One(), msize, ksize, nsize)
                        end
                        A = gesp(A, (msize, Zero()))
                        C = gesp(C, (msize, Zero()))
                    end
                end
                A = gesp(A, (Zero(), ksize))
                B = gesp(B, (ksize, Zero()))
                # synchronize on completion so we wait until every thread is done with `Bpacked` before beginning to overwrite it
                _mv = _atomic_add!(mys, one(UInt))
                let atoms = atomics
                    for _ ∈ CloseOpen(total_ids)
                        atoms += cache_linesize()
                        atoms == mys && continue
                        while _atomic_load(atoms) != sync_iters
                            pause()
                        end
                    end
                end
            end
        end
        B = gesp(B, (Zero(), nsize))
        C = gesp(C, (Zero(), nsize))
    end
    nothing
end


text:
function part1(filename)
    layers = collect(Iterators.partition(readchomp(filename), 25 * 6))
    _, i = findmin(layer -> count(==('0'), layer), layers)
    count(==('1'), layers[i]) * count(==('2'), layers[i])
end

function part2(filename)
    layers = Iterators.Stateful(Iterators.partition(readchomp(filename), 25 * 6))
    image = popfirst!(layers)
    while any(==('2'), image)
        image = [old == '2' ? new : old for (old, new) in zip(image, popfirst!(layers))]
    end
    open("image.pbm", "w") do f
        write(f, "P1\n25 6\n", join(join.(Iterators.partition(image, 25), ' '), '\n'))
    end
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module CompositeTypes

export iscomposite,
    components,
    ncomponents,
    component,
    setcomponent!

"""
    iscomposite(x)

Is `x` a composite object?

If an object is composite, then its components can be queried using the
`components` function.
"""
iscomposite(x) = length(components(x)) > 0

"Return the components of the composite object `x`."
components(x) = ()

"Return the number of components of composite object `x`."
ncomponents(x) = length(components(x))

"""
    component(x, I...)

Return the component of `x` that corresponds to the given index.
"""
component(x, I...) = getindex(components(x), I...)

"""
    setcomponent!(x, v, I...)

Set the component of `x` that corresponds to the given index to `v`.
"""
setcomponent!(x, v, I...) = setindex!(components(x), v, I...)
# the line above is a best-effort: it will not always work that way

include("Indexing.jl")
include("Display.jl")

end


text:
module MdSlicing
export fillfps, hyperslice

include("types.jl")
include("Meshes.jl")
include("HypersliceSets.jl")
include("Hyperslice.jl")

hyperslice = Hyperslice.hyperslice

#function fillfps(ps::ProblemSpec, fp::Array{Union{Missing,Float64},1}, n::Int = 50)
function fillfps(ps, fp, n = 50)
  if !any(map(ismissing, fp)) # special case
    return [fp]
  end
  fpvals =
    [ismissing(x[1]) ? (x[2][1] .+ (x[2][2]-x[2][1]) .* rand(n)) : fill(x[1], n) for x in zip(fp, ps.vals)]
  # transpose everything
  [[fpvals[c][r] for c=1:length(fp)] for r = 1:n]
end

end # module

max_stars_repo_name: gabysbrain/multid-slicing

max_issues_repo_name: gabysbrain/multid-slicing

max_forks_repo_name: gabysbrain/multid-slicing


text:
using Documenter
using DecisionProgramming

makedocs(
    sitename = "DecisionProgramming.jl",
    format = Documenter.HTML(
        assets = ["assets/favicon.ico"]
    ),
    modules = [DecisionProgramming],
    authors = "Jaan Tollander de Balsch",
    pages = [
        "Home" => "index.md",
        "Decision Programming" => Any[
            "decision-programming/influence-diagram.md",
            "decision-programming/paths.md",
            "decision-programming/decision-model.md",
            "decision-programming/analyzing-decision-strategies.md",
            "decision-programming/computational-complexity.md",
        ],
        "usage.md",
        "Examples" => Any[
            "examples/used-car-buyer.md",
            "examples/pig-breeding.md",
            "examples/n-monitoring.md",
            "examples/contingent-portfolio-programming.md",
            "examples/CHD_preventative_care.md",
        ],
        "api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/gamma-opt/DecisionProgramming.jl.git"
)


text:
abstract type AbstractLightHouse end
abstract type LossFunctionLightHouse <: AbstractLightHouse end

abstract type AbstractConformalTransform end

abstract type AbstractGameState end
abstract type AbstractPlayerState end
abstract type AbstractPlayerStateHistory end
abstract type AbstractWorldState end

abstract type AbstractConfig end

abstract type AbstractInputHandler end
@enum PLAYERINTENT moveforward=1 movebackward=2 moveleft=3 moveright=4 zoomout=5 zoomin=6 rotateleft=7 rotateright=8 plantlighthouse=9 stoprenderloop=10

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
foundflags = Calculations.documentation(QE, "noncollinear")
@test length(foundflags) == 2

@test eltype(Calculations.documentation(QE, :electron_maxstep)) == Int

@test length(Calculations.documentation(Elk, "magnetism")) == 2
@test Calculations.documentation(Elk, "magnetism")[1][2][1] == Calculations.documentation(Elk, :cmagz)
@test Calculations.documentation(Elk, :tasks) ==
      Calculations.getfirst(x -> x.name == :tasks, Calculations.ELK_CONTROLBLOCKS)
@test Calculations.documentation(Elk, "version")[1][1] ==
      Calculations.getfirst(x -> x.name == :tasks, Calculations.ELK_CONTROLBLOCKS)


text:
using Gtk.GLib
using Gtk
using Gtk.ShortNames
using Test

@testset "glist" begin

@testset "pointers" begin

w = Window("Window", 400, 300)
nb = Notebook()
w = push!(Window("Notebook"),nb)
l = ccall((:gtk_container_get_children,Gtk.libgtk),Ptr{Gtk._GList{Gtk.GtkWidget}},(Ptr{Gtk.GObject},),w)

@test eltype(l)==Gtk.GtkWidget

@test !isempty(l)
@test length(l)==1
@test l[1]==nb

for item in l
    @test item==nb
end

end

@testset "string" begin

g = GList(String)

@test ndims(g)==1

@test isempty(g)

for i=1:10
    push!(g,string(i))
end

@test !isempty(g)

@test g[1]==string(1)
@test g[2]==string(2)

@test length(g)==10
@test g[end]==string(10)

for (i,item) = enumerate(g)
    @test string(i)==item
end

g2=copy(g)

for (i,item) = enumerate(g2)
    @test string(i)==item
end

@test popfirst!(g2)==string(1)
@test length(g2)==9
@test pop!(g2)==string(10)
@test length(g2)==8

@test length(g)==10

insert!(g,8,string(25))
@test length(g)==11
@test g[8]==string(25)
@test g[9]==string(8)

pushfirst!(g,string(24))
@test length(g)==12
@test g[1]==string(24)
@test g[9]==string(25)

g[8]="hello"
@test g[8]=="hello"
@test g[9]==string(25)

reverse!(g)
@test g[end]==string(24)

g[2]="test"
@test g[2]=="test"

g3=copy(g)

append!(g2,g3)
@test length(g2)==20

empty!(g)
@test isempty(g)

# repeat above for GSList
g = GSList(String)

@test ndims(g)==1

@test isempty(g)

for i=1:10
    push!(g,string(i))
end

@test !isempty(g)

@test g[1]==string(1)
@test g[2]==string(2)

@test length(g)==10
@test g[end]==string(10)

for (i,item) = enumerate(g)
    @test string(i)==item
end

g2=copy(g)

for (i,item) = enumerate(g2)
    @test string(i)==item
end

@test popfirst!(g2)==string(1)
@test length(g2)==9
@test pop!(g2)==string(10)
@test length(g2)==8

@test length(g)==10

insert!(g,8,string(25))
@test length(g)==11
@test g[8]==string(25)
@test g[9]==string(8)

pushfirst!(g,string(24))
@test length(g)==12
@test g[1]==string(24)
@test g[9]==string(25)

g[8]="hello"
@test g[8]=="hello"
@test g[9]==string(25)

reverse!(g)
@test g[end]==string(24)

g[2]="test"
@test g[2]=="test"

g3=copy(g)

append!(g2,g3)
@test length(g2)==20

empty!(g)
@test isempty(g)

end

@testset "numbers" begin

g = GList(Int64)

@test ndims(g)==1

@test isempty(g)

for i=1:10
    push!(g,i)
end

@test !isempty(g)

@test g[1]==1
@test g[2]==2

@test length(g)==10
@test g[end]==10

for (i,item) = enumerate(g)
    @test i==item
end

g2=copy(g)

for (i,item) = enumerate(g2)
    @test i==item
end

empty!(g)
@test isempty(g)

# repeat above for SList
g = GSList(Int64)

@test ndims(g)==1

@test isempty(g)

for i=1:10
    push!(g,i)
end

@test !isempty(g)

@test g[1]==1
@test g[2]==2

@test length(g)==10
@test g[end]==10

for (i,item) = enumerate(g)
    @test i==item
end

g2=copy(g)

for (i,item) = enumerate(g2)
    @test i==item
end

empty!(g)
@test isempty(g)

end

GC.gc()

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner