Added logistic regression example

Author: Mark Poscablo

Diff for: .gitmodules

@@ -0,0 +1,3 @@
+[submodule "assets"]
+	path = assets
+	url = https://github.com/arrayfire/assets.git

Diff for: assets

@@ -0,0 +1,64 @@
+#!/usr/bin/python
+
+#######################################################
+# Copyright (c) 2019, ArrayFire
+# All rights reserved.
+#
+# This file is distributed under 3-clause BSD license.
+# The complete license agreement can be obtained at:
+# http://arrayfire.com/licenses/BSD-3-Clause
+########################################################
+
+def reverse_char(b):
+    b = (b & 0xF0) >> 4 | (b & 0x0F) << 4
+    b = (b & 0xCC) >> 2 | (b & 0x33) << 2
+    b = (b & 0xAA) >> 1 | (b & 0x55) << 1
+    return b
+
+
+# http://stackoverflow.com/a/9144870/2192361
+def reverse(x):
+    x = ((x >> 1)  & 0x55555555) | ((x & 0x55555555) << 1)
+    x = ((x >> 2)  & 0x33333333) | ((x & 0x33333333) << 2)
+    x = ((x >> 4)  & 0x0f0f0f0f) | ((x & 0x0f0f0f0f) << 4)
+    x = ((x >> 8)  & 0x00ff00ff) | ((x & 0x00ff00ff) << 8)
+    x = ((x >> 16) & 0xffff)     | ((x & 0xffff)     << 16);
+    return x
+
+
+def read_idx(name):
+    with open(name, 'rb') as f:
+        # In the C++ version, bytes the size of 4 chars are being read
+        # May not work properly in machines where a char is not 1 byte
+        bytes_read = f.read(4)
+        bytes_read = bytearray(bytes_read)
+
+        if bytes_read[2] != 8:
+            raise RuntimeError('Unsupported data type')
+
+        numdims = bytes_read[3]
+        elemsize = 1
+
+        # Read the dimensions
+        elem = 1
+        dims = [0] * numdims
+        for i in range(numdims):
+            bytes_read = bytearray(f.read(4))
+
+            # Big endian to little endian
+            for j in range(4):
+                bytes_read[j] = reverse_char(bytes_read[j])
+            bytes_read_int = int.from_bytes(bytes_read, 'little')
+            dim = reverse(bytes_read_int)
+
+            elem = elem * dim;
+            dims[i] = dim;
+
+        # Read the data
+        cdata = f.read(elem * elemsize)
+        cdata = list(cdata)
+        data = [float(cdata_elem) for cdata_elem in cdata]
+
+        return (dims, data)
+
+

Diff for: examples/common/idxio.py

@@ -0,0 +1,205 @@
+#!/usr/bin/python
+
+#######################################################
+# Copyright (c) 2019, ArrayFire
+# All rights reserved.
+#
+# This file is distributed under 3-clause BSD license.
+# The complete license agreement can be obtained at:
+# http://arrayfire.com/licenses/BSD-3-Clause
+########################################################
+
+from mnist_common import display_results, setup_mnist
+
+import sys
+import time
+
+import arrayfire as af
+from arrayfire.algorithm import max, imax, count, sum
+from arrayfire.arith import abs, sigmoid, log
+from arrayfire.array import transpose
+from arrayfire.blas import matmul, matmulTN
+from arrayfire.data import constant, join, lookup, moddims
+from arrayfire.device import set_device, sync, eval
+
+
+def accuracy(predicted, target):
+    _, tlabels = af.imax(target, 1)
+    _, plabels = af.imax(predicted, 1)
+    return 100 * af.count(plabels == tlabels) / tlabels.elements()
+
+
+def abserr(predicted, target):
+    return 100 * af.sum(af.abs(predicted - target)) / predicted.elements()
+
+
+# Predict (probability) based on given parameters
+def predict_proba(X, Weights):
+    Z = af.matmul(X, Weights)
+    return af.sigmoid(Z)
+
+
+# Predict (log probability) based on given parameters
+def predict_log_proba(X, Weights):
+    return af.log(predict_proba(X, Weights))
+
+
+# Give most likely class based on given parameters
+def predict(X, Weights):
+    probs = predict_proba(X, Weights)
+    _, classes = af.imax(probs, 1)
+    return classes
+
+
+def cost(Weights, X, Y, lambda_param=1.0):
+    # Number of samples
+    m = Y.dims()[0]
+
+    dim0 = Weights.dims()[0]
+    dim1 = Weights.dims()[1] if len(Weights.dims()) > 1 else None
+    dim2 = Weights.dims()[2] if len(Weights.dims()) > 2 else None
+    dim3 = Weights.dims()[3] if len(Weights.dims()) > 3 else None
+    # Make the lambda corresponding to Weights(0) == 0
+    lambdat = af.constant(lambda_param, dim0, dim1, dim2, dim3)
+
+    # No regularization for bias weights
+    lambdat[0, :] = 0
+
+    # Get the prediction
+    H = predict_proba(X, Weights)
+
+    # Cost of misprediction
+    Jerr = -1 * af.sum(Y * af.log(H) + (1 - Y) * af.log(1 - H), dim=0)
+
+    # Regularization cost
+    Jreg = 0.5 * af.sum(lambdat * Weights * Weights, dim=0)
+
+    # Total cost
+    J = (Jerr + Jreg) / m
+
+    # Find the gradient of cost
+    D = (H - Y)
+    dJ = (af.matmulTN(X, D) + lambdat * Weights) / m
+
+    return J, dJ
+
+
+def train(X, Y, alpha=0.1, lambda_param=1.0, maxerr=0.01, maxiter=1000, verbose=False):
+    # Initialize parameters to 0
+    Weights = af.constant(0, X.dims()[1], Y.dims()[1])
+
+    for i in range(maxiter):
+        # Get the cost and gradient
+        J, dJ = cost(Weights, X, Y, lambda_param)
+
+        err = af.max(af.abs(J))
+        if err < maxerr:
+            print('Iteration {0:4d} Err: {1:4f}'.format(i + 1, err))
+            print('Training converged')
+            return Weights
+
+        if verbose and ((i+1) % 10 == 0):
+            print('Iteration {0:4d} Err: {1:4f}'.format(i + 1, err))
+
+        # Update the parameters via gradient descent
+        Weights = Weights - alpha * dJ
+
+    if verbose:
+        print('Training stopped after {0:d} iterations'.format(maxiter))
+
+    return Weights
+
+
+def benchmark_logistic_regression(train_feats, train_targets, test_feats):
+    t0 = time.time()
+    Weights = train(train_feats, train_targets, 0.1, 1.0, 0.01, 1000)
+    af.eval(Weights)
+    sync()
+    t1 = time.time()
+    dt = t1 - t0
+    print('Training time: {0:4.4f} s'.format(dt))
+
+    t0 = time.time()
+    iters = 100
+    for i in range(iters):
+        test_outputs = predict(test_feats, Weights)
+        af.eval(test_outputs)
+    sync()
+    t1 = time.time()
+    dt = t1 - t0
+    print('Prediction time: {0:4.4f} s'.format(dt / iters))
+
+
+# Demo of one vs all logistic regression
+def logit_demo(console, perc):
+    # Load mnist data
+    frac = float(perc) / 100.0
+    mnist_data = setup_mnist(frac, True)
+    num_classes = mnist_data[0]
+    num_train = mnist_data[1]
+    num_test = mnist_data[2]
+    train_images = mnist_data[3]
+    test_images = mnist_data[4]
+    train_targets = mnist_data[5]
+    test_targets = mnist_data[6]
+
+    # Reshape images into feature vectors
+    feature_length = int(train_images.elements() / num_train);
+    train_feats = af.transpose(af.moddims(train_images, feature_length, num_train))
+    test_feats = af.transpose(af.moddims(test_images, feature_length, num_test))
+
+    train_targets = af.transpose(train_targets)
+    test_targets = af.transpose(test_targets)
+
+    num_train = train_feats.dims()[0]
+    num_test = test_feats.dims()[0]
+
+    # Add a bias that is always 1
+    train_bias = af.constant(1, num_train, 1)
+    test_bias = af.constant(1, num_test, 1)
+    train_feats = af.join(1, train_bias, train_feats)
+    test_feats = af.join(1, test_bias, test_feats)
+
+    # Train logistic regression parameters
+    Weights = train(train_feats, train_targets,
+                    0.1,  # learning rate
+                    1.0,  # regularization constant
+                    0.01, # max error
+                    1000, # max iters
+                    True  # verbose mode
+    )
+    af.eval(Weights)
+    af.sync()
+
+    # Predict the results
+    train_outputs = predict_proba(train_feats, Weights)
+    test_outputs = predict_proba(test_feats, Weights)
+
+    print('Accuracy on training data: {0:2.2f}'.format(accuracy(train_outputs, train_targets)))
+    print('Accuracy on testing data: {0:2.2f}'.format(accuracy(test_outputs, test_targets)))
+    print('Maximum error on testing data: {0:2.2f}'.format(abserr(test_outputs, test_targets)))
+
+    benchmark_logistic_regression(train_feats, train_targets, test_feats)
+
+    if not console:
+        test_outputs = af.transpose(test_outputs)
+        # Get 20 random test images
+        display_results(test_images, test_outputs, af.transpose(test_targets), 20, True)
+
+def main():
+    argc = len(sys.argv)
+
+    device  = int(sys.argv[1])      if argc > 1 else 0
+    console = sys.argv[2][0] == '-' if argc > 2 else False
+    perc    = int(sys.argv[3])      if argc > 3 else 60
+
+    try:
+        af.set_device(device)
+        af.info()
+        logit_demo(console, perc)
+    except Exception as e:
+        print('Error: ', str(e))
+
+
+if __name__ == '__main__':
+    main()

Diff for: examples/machine_learning/logistic_regression.py

@@ -0,0 +1,101 @@
+#!/usr/bin/python
+
+#######################################################
+# Copyright (c) 2019, ArrayFire
+# All rights reserved.
+#
+# This file is distributed under 3-clause BSD license.
+# The complete license agreement can be obtained at:
+# http://arrayfire.com/licenses/BSD-3-Clause
+########################################################
+
+import sys
+sys.path.insert(0, '../common')
+from idxio import read_idx
+
+import arrayfire as af
+from arrayfire.algorithm import where
+from arrayfire.array import Array
+from arrayfire.data import constant, lookup, moddims
+from arrayfire.random import randu
+
+
+def classify(arr, k, expand_labels):
+    ret_str = ''
+    if expand_labels:
+        vec = arr[:, k].as_type(af.Dtype.f32)
+        h_vec = vec.to_list()
+        data = []
+
+        for i in range(vec.elements()):
+            data.append((h_vec[i], i))
+
+        data = sorted(data, key=lambda pair: pair[0], reverse=True)
+
+        ret_str = str(data[0][1])
+
+    else:
+        ret_str = str(int(arr[k].as_type(af.Dtype.f32).scalar()))
+
+    return ret_str
+
+
+def setup_mnist(frac, expand_labels):
+    idims, idata = read_idx('../../assets/examples/data/mnist/images-subset')
+    ldims, ldata = read_idx('../../assets/examples/data/mnist/labels-subset')
+
+    idims.reverse()
+    numdims = len(idims)
+    images = af.Array(idata, tuple(idims))
+
+    R = af.randu(10000, 1);
+    cond = R < min(frac, 0.8)
+    train_indices = af.where(cond)
+    test_indices = af.where(~cond)
+
+    train_images = af.lookup(images, train_indices, 2) / 255
+    test_images = af.lookup(images, test_indices, 2) / 255
+
+    num_classes = 10
+    num_train = train_images.dims()[2]
+    num_test = test_images.dims()[2]
+
+    if expand_labels:
+        train_labels = af.constant(0, num_classes, num_train)
+        test_labels = af.constant(0, num_classes, num_test)
+
+        h_train_idx = train_indices.to_list()
+        h_test_idx = test_indices.to_list()
+
+        for i in range(num_train):
+            train_labels[ldata[h_train_idx[i]], i] = 1
+
+        for i in range(num_test):
+            test_labels[ldata[h_test_idx[i]], i] = 1
+
+    else:
+        labels = af.Array(ldata, tuple(ldims))
+        train_labels = labels[train_indices]
+        test_labels = labels[test_indices]
+
+    return (num_classes,
+            num_train,
+            num_test,
+            train_images,
+            test_images,
+            train_labels,
+            test_labels)
+
+
+def display_results(test_images, test_output, test_actual, num_display, expand_labels):
+    for i in range(num_display):
+        print('Predicted: ', classify(test_output, i, expand_labels))
+        print('Actual: ', classify(test_actual, i, expand_labels))
+
+        img = (test_images[:, :, i] > 0.1).as_type(af.Dtype.u8)
+        img = af.moddims(img, img.elements()).to_list()
+        for j in range(28):
+            for k in range(28):
+                print('\u2588' if img[j * 28 + k] > 0 else ' ', end='')
+            print()
+        input()