examples/machine_learning/pca.cpp

#include <stdio.h>
#include <arrayfire.h>

#include "wines_data.h" // Dataset

using namespace af;
using namespace std;


// data:  input,  [n d] - n samples, d dimensions
// pvec:  output, principal vectors
// eigv:  output, eigen values
// pcomp: output, principal components
// norm:  output, normalized input data
void pca(array& data, array& pvec, array& eigv, array& pcomp, array& norm) {

    // sizes
    int n = data.dims()[0];

    // data normalization
    array means = mean(data);
    array stdvs = stdev(data);
    array sdata = (data - tile(means, n, 1)) / tile(stdvs, n, 1);

    // decomposition
    array V, D, E;
    eigen(V, D, cov(sdata));
    E = diag(D);
    array val, idx;
    sort(val, idx, -1*E);
    array comps = matmul(sdata, V); // principal components

    // output
    pvec = V;                 // principal vectors
    eigv = E(idx);            // sorted eigen values
    pcomp = comps(span, idx); // sorted principal components
    norm = sdata;             // normalized input

}


// data:  input,  [n m] - n features, m samples
// pvec:  output, principal vectors
// eigv:  output, eigen values
// pcomp: output, principal components
void pca(array& data, array& pvec, array& eigv, array& pcomp) {
    array norm;
    pca(data, pvec, eigv, pcomp, norm);
}


// Computes PCA on a wine dataset of 14 dimensions,
// and reduces this down to the top 3 for visualization.
void pca_demo(bool console) {
    array wines = array(14, 178, wines_h);
    array data = wines.T();

    array pvec, eigv, pcmp, norm;
    pca(data, pvec, eigv, pcmp, norm);
    print(eigv);

    array trim = pvec(span, seq(4));
    array proj = matmul(pcmp, trim);

    if (!console) {
        fig("color","gray");

        fig("sub",2,1,1); fig("title","Input");
        plot3(data(span,1), data(span,2), data(span,3));

        fig("sub",2,1,2); fig("title","14D->3D PCA");
        plot3(proj(span,1), proj(span,2), proj(span,3));

        fig("draw");
    }
}


// Entry point.
int main(int argc, char** argv) {

    bool console = false;
    if (argc > 2 || (argc == 2 && argv[1][0] != '-')) {
        printf("usage: pca_demo [-]\n");
        return -1;
    } else if (argc == 2 && argv[1][0] == '-') {
        console = true;
    }

    try {
        printf("** ArrayFire PCA Demo **\n\n");
        pca_demo(console);

    } catch (af::exception& e) {
        fprintf(stderr, "%s\n", e.what());
    }

    if (!console) {
        printf("hit [enter]...");
        getchar();
    }
    return 0;
}