Add PCA benchmarks (#5)

* Add native, scikit-learn, daal4py PCA benchmark.

* Add NPY serializer

Author: bibikar

daal4py/pca.py

@@ -0,0 +1,145 @@
+# Copyright (C) 2017-2019 Intel Corporation
+#
+# SPDX-License-Identifier: MIT
+
+from __future__ import print_function
+import numpy as np
+import timeit
+from numpy.random import rand
+from daal4py import pca, pca_transform, normalization_zscore
+from daal4py.sklearn.utils import getFPType
+from sklearn.utils.extmath import svd_flip
+from args import getArguments, coreString
+from bench import prepare_benchmark
+
+import argparse
+parser = argparse.ArgumentParser(description='daal4py PCA benchmark')
+parser.add_argument('--svd-solver', type=str, choices=['daal', 'full'],
+                    default='daal', help='SVD solver to use')
+parser.add_argument('--n-components', type=int, default=None,
+                    help='Number of components to find')
+parser.add_argument('--whiten', action='store_true', default=False,
+                    help='Perform whitening')
+parser.add_argument('--write-results', action='store_true', default=False,
+                    help='Write results to disk for verification')
+args = getArguments(parser)
+REP = args.iteration if args.iteration != '?' else 10
+core_number, daal_version = prepare_benchmark(args)
+
+
+def st_time(func):
+    def st_func(*args, **keyArgs):
+        times = []
+        for n in range(REP):
+            t1 = timeit.default_timer()
+            r = func(*args, **keyArgs)
+            t2 = timeit.default_timer()
+            times.append(t2-t1)
+        print(min(times))
+        return r
+    return st_func
+
+p = args.size[0]
+n = args.size[1]
+X = rand(p,n)
+Xp = rand(p,n)
+
+if args.n_components:
+    n_components = args.n_components
+else:
+    n_components = min((n, (2 + min((n, p))) // 3))
+
+def pca_fit_daal(X, n_components):
+
+    if n_components < 1:
+        n_components = min(X.shape)
+
+    fptype = getFPType(X)
+
+    centering_algo = normalization_zscore(
+            fptype=fptype,
+            doScale=False
+    )
+
+    pca_algorithm = pca(
+            fptype=fptype,
+            method='svdDense',
+            normalization=centering_algo,
+            resultsToCompute='mean|variance|eigenvalue',
+            isDeterministic=True,
+            nComponents=n_components
+    )
+
+    pca_result = pca_algorithm.compute(X)
+    eigenvectors = pca_result.eigenvectors
+    eigenvalues = pca_result.eigenvalues.ravel()
+    singular_values = np.sqrt((X.shape[0] - 1) * eigenvalues)
+
+    return pca_result, eigenvalues, eigenvectors, singular_values
+
+
+def pca_transform_daal(pca_result, X, n_components, fit_n_samples,
+                       eigenvalues, eigenvectors,
+                       whiten=False, scale_eigenvalues=False):
+
+    fptype = getFPType(X)
+
+    tr_data = {}
+    tr_data['mean'] = pca_result.dataForTransform['mean']
+
+    if whiten:
+        if scale_eigenvalues:
+            tr_data['eigenvalue'] = (fit_n_samples - 1) * pca_result.eigenvalues
+        else:
+            tr_data['eigenvalue'] = pca_result.eigenvalues
+    elif scale_eigenvalues:
+        tr_data['eigenvalue'] = np.full((1, pca_result.eigenvalues.size),
+                                        fit_n_samples - 1, dtype=X.dtype)
+
+    transform_algorithm = pca_transform(fptype=fptype, nComponents=n_components)
+    transform_result = transform_algorithm.compute(X, pca_result.eigenvectors, tr_data)
+    return transform_result.transformedData
+
+
+def pca_fit_full_daal(X, n_components):
+
+    fit_result, eigenvalues, eigenvectors, S = pca_fit_daal(X, min(X.shape))
+    U = pca_transform_daal(fit_result, X, min(X.shape), X.shape[0],
+                           eigenvalues, eigenvectors,
+                           whiten=True, scale_eigenvalues=True)
+    V = fit_result.eigenvectors
+
+    U, V = svd_flip(U, V)
+
+    eigenvalues = fit_result.eigenvalues[:n_components].copy()
+    eigenvectors = fit_result.eigenvectors[:n_components].copy()
+
+    return fit_result, eigenvalues, eigenvectors, U, S, V
+
+
+@st_time
+def test_fit(X):
+    if args.svd_solver == 'full':
+        return pca_fit_full_daal(X, n_components)
+    else:
+        return pca_fit_daal(X, n_components)
+
+@st_time
+def test_transform(Xp, pca_result, eigenvalues, eigenvectors):
+    return pca_transform_daal(pca_result, Xp, n_components, X.shape[0],
+                              eigenvalues, eigenvectors,
+                              whiten=args.whiten)
+
+
+print (','.join([args.batchID, args.arch, args.prefix, "PCA.fit", coreString(args.num_threads), "Double", "%sx%s" % (p,n)]), end=',')
+res = test_fit(X)
+print (','.join([args.batchID, args.arch, args.prefix, "PCA.transform", coreString(args.num_threads), "Double", "%sx%s" % (p,n)]), end=',')
+tr = test_transform(Xp, res[0], res[1], res[2])
+
+
+if args.write_results:
+    np.save('pca_daal4py_X.npy', X)
+    np.save('pca_daal4py_Xp.npy', Xp)
+    np.save('pca_daal4py_eigvals.npy', res[1])
+    np.save('pca_daal4py_eigvecs.npy', res[2])
+    np.save('pca_daal4py_tr.npy', tr)

native/common.hpp

@@ -7,6 +7,7 @@
 
 #include "CLI11.hpp"
 #include "daal.h"
+#include "npyfile.h"
 
 namespace dm = daal::data_management;
 namespace ds = daal::services;
@@ -133,6 +134,76 @@ dm::NumericTablePtr make_table(T *data, size_t rows, size_t cols) {
 }
 
 
+/*
+ * Write a NumericTable to npy format.
+ */
+template <typename T>
+void write_table(dm::NumericTablePtr table,
+                 const std::string descr, const std::string fn) {
+
+    dm::BlockDescriptor<T> block;
+    table->getBlockOfRows(0, table->getNumberOfRows(), dm::readOnly, block);
+    T *data = block.getBlockPtr();
+
+    size_t shape[2] = {table->getNumberOfRows(), table->getNumberOfColumns()};
+
+    char *c_descr = new char[descr.size() + 1];
+    strcpy(c_descr, descr.c_str());
+    const struct npyarr arr {
+        .descr = c_descr,
+        .fortran_order = false,
+        .shape_len = 2,
+        .shape = shape,
+        .data = data
+    };
+
+    save_npy(&arr, fn.c_str(), sizeof(T));
+
+    table->releaseBlockOfRows(block);
+
+}
+
+
+/*
+ * Create a new HomogenNumericTable from a submatrix of an existing
+ * HomogenNumericTable.
+ *
+ * equivalent in python:
+ *   return src[row_start:row_end, col_start:col_end].copy()
+ */
+template <typename T>
+dm::NumericTablePtr
+copy_submatrix(dm::NumericTablePtr src,
+               size_t row_start, size_t row_end,
+               size_t col_start, size_t col_end) {
+
+    dm::BlockDescriptor<T> srcblock;
+    src->getBlockOfRows(0, src->getNumberOfRows(), dm::readOnly, srcblock);
+    T *srcdata = srcblock.getBlockPtr();
+
+    size_t cols = col_end - col_start;
+    size_t rows = row_end - row_start;
+
+    dm::NumericTablePtr dest = dm::HomogenNumericTable<T>::create(
+            cols, rows, dm::NumericTable::doAllocate);
+    dm::BlockDescriptor<T> destblock;
+    dest->getBlockOfRows(0, dest->getNumberOfRows(), dm::readWrite, destblock);
+    T *destdata = destblock.getBlockPtr();
+
+    for (int i = 0; i < rows; i++) {
+        for (int j = 0; j < cols; j++) {
+            destdata[i*cols + j] = srcdata[(i+col_start)*cols + j+row_start];
+        }
+    }
+
+    src->releaseBlockOfRows(srcblock);
+    dest->releaseBlockOfRows(destblock);
+
+    return dest;
+
+}
+
+
 /*
  * Generate an array of random numbers.
  */

native/npyfile.h

@@ -10,6 +10,7 @@
 #include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 
 #ifndef _NPYFILE_H_
 #define _NPYFILE_H_
@@ -281,6 +282,105 @@ struct npyarr *load_npy(const char *path) {
     free_npy(arr);
     return NULL;
 }
+
+
+/*
+ * Write an array to disk.
+ * N.B. elem_size must be the correct size of each element in the
+ * array.
+ */
+void save_npy(const struct npyarr *arr, const char *path, size_t elem_size) {
+
+    if (path == NULL || arr == NULL) {
+        return;
+    }
+
+    /* open file */
+    FILE *f;
+    if ((f = fopen(path, "wb")) == NULL) {
+        return;
+    }
+
+    /* Determine header_len */
+    unsigned int header_len = 0;
+    static const char header_1[] = "{'descr': '";
+    static const char header_2[] = "', 'fortran_order': ";
+    static const char header_3[] = ", 'shape': (";
+    static const char header_4[] = ")}";
+    header_len += strlen(header_1) + strlen(header_2);
+    header_len += strlen(header_3) + strlen(header_4);
+
+    header_len += strlen(arr->descr);
+    header_len += (arr->fortran_order) ? 4 : 5;
+    header_len += 1; /* the terminating newline */
+
+    /* Determine shape length */
+    unsigned int shape_len = 0;
+    for (int i = 0; i < arr->shape_len; i++) {
+        /* Number of digits in the number */
+        if (arr->shape[i] > 0) {
+            shape_len += floor(log10(arr->shape[i]) + 1);
+        } else {
+            shape_len += 1;
+        }
+        if (i < arr->shape_len - 1 && arr->shape_len != 1) {
+            /* Adding commas where needed */
+            shape_len += 2;
+        }
+    }
+
+    header_len += shape_len;
 
+    /* 16-byte alignment */
+    unsigned int bytes_left = 16 - (strlen(NPY_HEADER) + 4 + header_len) % 16;
+    header_len += bytes_left;
+
+    unsigned char major_version = 0x01;
+    if (header_len > 65535) {
+        major_version = 0x02;
+    }
+    unsigned char minor_version = 0x00;
+
+    fputs(NPY_HEADER, f);
+    fputc(major_version, f);
+    fputc(minor_version, f);
+    
+    /* writing header_len in little-endian format always */
+    for (int i = 0; i < ((major_version >= 0x02) ? 4 : 2); i++) {
+        fputc((header_len >> (i * 8)) & 0xff, f);
+    }
+
+    fputs(header_1, f);
+    fputs(arr->descr, f);
+    fputs(header_2, f);
+    fputs((arr->fortran_order) ? "True" : "False", f);
+    fputs(header_3, f);
+
+    size_t nelem = 1;
+    for (int i = 0; i < arr->shape_len; i++) {
+        fprintf(f, "%d", arr->shape[i]);
+        nelem *= arr->shape[i];
+        if (i < arr->shape_len - 1 && arr->shape_len != 1) {
+            fputs(", ", f);
+        }
+    }
+    fputs(header_4, f);
+
+    /* Pad with spaces */
+    for (int i = 0; i < bytes_left; i++) {
+        fputc(' ', f);
+    }
+
+    /* Add terminating newline */
+    fputc('\n', f);
+
+    /* Write array data */
+    fwrite(arr->data, elem_size, nelem, f);
+
+    /* done. */
+    fflush(f);
+    fclose(f);
+
+}
 
 #endif /* _NPYFILE_H_ */