In [1]:
%matplotlib inline
import numpy
import matplotlib.pyplot as pyplot
from km3net.kernels import QuadraticDifferenceSparse, PurgingSparse
import km3net.util as util
print ("We're good to go!")
Let's start with reading some input data. We assume the data is stored in a ASCI file, using four columns. The first column stores the time in nano seconds, column two through four store the x,y,z-coordinates of the hits. We assume that the input data is sorted on time, with the earliest hit first. Finally, we assume a file contains data from only one timeslice. For this we are going to use the get_real_input_data utility function.
In [2]:
N,x,y,z,ct = util.get_real_input_data("sample.txt")
print ("Read", N, "hits from file")
Before we continue, we have to initialize the GPU kernels used by the software. This means we have to instantiate the Python objects that will serve as an interface to our compiled GPU kernels and the data currently stored in the memory of the GPU.
In [3]:
context, cc = util.init_pycuda()
qd_kernel = QuadraticDifferenceSparse(N, cc=cc)
purging = PurgingSparse(N, cc)
The next step is of course start doing some computations. We want to compute all correlations between the hits we've just retrieved from a file. For this we use the 'qd_kernel' name we've just created to refer to the GPU Kernel that computes correlations based on the Quadratic Difference criterion and stores the result as a sparse matrix. We use the kernel by means of its compute function.
In [4]:
d_col_idx, d_prefix_sums, d_degrees, total_hits = qd_kernel.compute(x, y, z, ct)
Now let's take a look at the data we got back from the GPU. We start by converting the sparse matrix stored on the GPU to a dense matrix stored on the host as a Numpy array. Then we make a quick plot of the correlation matrix.
In [5]:
dense_matrix = util.sparse_to_dense(d_prefix_sums, d_col_idx, N, total_hits)
pyplot.imshow(dense_matrix, cmap=pyplot.cm.bone, interpolation='nearest')
Out[5]:
The next step is of course to look for a group of hits that are all correlated with eachother. We approximate this clique using the purging algorithm. The implementation that we will use is the PurgingSparse GPU kernel that performs the purging algorithm directly on the sparsely stored correlation matrix in GPU memory.
We can call the GPU implementation of PurgingSparse using the name 'purging', which we created earlier. PurgingSparse implements a method called compute, which we will use to compute the clique.
In [7]:
clique = purging.compute(d_col_idx, d_prefix_sums, d_degrees)
print("found clique of size", len(clique))
print(clique)
And that's it! This example has shown how to read input data from a file and call the two main components of the pipeline.