Superceded, please disregard. Chapel Kernel Matrix Benchmark (WIP)

kernelmatrix.chpl

use CPtr;
use Math;

record DotProduct {
  type T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var ret: T = 0: T;
    for i in 0..#p {
      ret += xrow[i] * yrow[i];
    }
    return ret;
  }
}

record Gaussian {
  type T;
  const theta: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      var tmp = xrow[i] - yrow[i];
      dist += tmp * tmp;
    }
    return exp(-sqrt(dist)/this.theta);
  }
}

record Polynomial {
  type T;
  const d: T;
  const offset: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      dist += xrow[i]*yrow[i];
    }
    return (dist + this.offset)**this.d;
  }
}

record Exponential {
  type T;
  const theta: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      dist += abs(xrow[i] - yrow[i]);
    }
    return exp(dist/this.theta);
  }
}

record Log {
  type T;
  const beta: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      dist += abs(xrow[i] - yrow[i])**this.beta;
    }
    dist = dist**(1/this.beta);
    return -log(1 + dist);
  }
}

record Cauchy {
  type T;
  const theta: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      var tmp = xrow[i] - yrow[i];
      dist += tmp * tmp;
    }
    dist = sqrt(dist)/this.theta;
    return 1/(1 + dist);
  }
}

record Power {
  type T;
  const beta: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      dist += abs(xrow[i] - yrow[i])**this.beta;
    }
    return -dist**(1/this.beta);
  }
}

record Wave {
  type T;
  const theta: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      dist += abs(xrow[i] - yrow[i]);
    }
    var tmp = this.theta/dist;
    return tmp*sin(1/tmp);
  }
}

record Sigmoid {
  type T;
  const beta0: T;
  const beta1: T;
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var dist: T = 0: T;
    for i in 0..#p {
      dist += xrow[i] * yrow[i];
    }
    return tanh(this.beta0 * dist + this.beta1);
  }
}

/***************************************************************************/
use DynamicIters;
proc calculateKernelMatrix(K, data: [?D] ?T) /* : [?E] T */
{
  var n = D.dim(0).last;
  var p = D.dim(1).last;
  var E: domain(2) = {D.dim(0), D.dim(0)};
  var mat: [E] T;
  // code below assumes data starts at 1,1
  var rowPointers: [1..n] c_ptr(T) =
    forall i in 1..n do c_ptrTo(data[i, 1]);

  forall j in guided(1..n by -1) {
    for i in j..n {
      mat[i, j] = K.kernel(rowPointers[i], rowPointers[j], p);
      mat[j, i] = mat[i, j];
    }
  }
  return mat;
}

script.chpl

use kernelmatrix;

use IO;
use Time;
use Random;

proc bench(type T, Kernel, n: [?D] int(64), verbose: bool = true)
{
  var nitems: int(64) = D.dim(0).last: int(64);
  var times: [0..nitems] real(64);
  for i in 0..nitems {
    var _times: [1..3] real(64);
    var data: [1..n[i], 1..784] T;
    fillRandom(data);
    for j in 1..3 {
      var sw = new Timer();
      sw.start();
      var mat = calculateKernelMatrix(Kernel, data);
      sw.stop();
      _times[j] = (sw.elapsed(TimeUnits.microseconds)/1000_000): real(64);
    }
    times[i] = (_times[1] + _times[2] + _times[3])/3;
    if verbose {
      writeln("Average time for n = ", n[i], ", ", times[i], " seconds.");
      writeln("Detailed times: ", _times);
    }
  }
  return (n, times);
}

proc runKernelBenchmarks(type T, kernels, n: [?D] int(64), verbose: bool = true)
{
  var tmp = bench(T, kernels(0), n, verbose);
  type R = tmp.type;
  var results: [0..kernels.size] R;
  results[0] = tmp;
  var i = 0;
  for kernel in kernels {
    if verbose {
      writeln("\n\nRunning benchmarks for ", kernel: string);
    }
    results[i] = bench(T, kernel, n, verbose);
    i += 1;
  }
  return results;
}

/**
  To compile:
  chpl script.chpl kernelmatrix.chpl --fast && ./script
  bench: 0.040345 0.978326 3.85939 16.7087 40.2499
*/
proc runAllKernelBenchmarks(type T, verbose: bool = true)
{
  var n = [100, 500, 100];
  //var n = [1000, 5000, 10_000, 20_000, 30_000];
  //writeln("bench: ", bench(T, K, n));
  var kernels = (new DotProduct(T),        new Gaussian(T, 1: T), new Polynomial(T, 2.5: T, 1: T),
                 new Exponential(T, 1: T), new Log(T, 3: T),      new Cauchy(T, 1: T),
                 new Power(T, 2.5: T),     new Wave(T, 1: T),     new Sigmoid(T, 1: T, 1: T));
  var kernelNames = ["DotProduct",  "Gaussian", "Polynomial",
                       "Exponential", "Log",      "Cauchy",
                       "Power",       "Wave",     "Sigmoid"];
  var results = runKernelBenchmarks(T, kernels, n);
  

  var last: int(64) = n.domain.dim(0).last: int(64);
  var tabLen = n.size * kernels.size - 1;
  var table: [0..tabLen, 0..3] string;
  table[0, ..] = ["language", "kernel", "nitems", "time"];
  while (true)
  {
    var k = 1;
    for i in 0..kernels.size {
      var tmp = ["Chapel", kernelNames[i], "", ""];
      for j in 0..n.size {
        tmp[2] = results[i](0)[j]: string;
        tmp[3] = results[i](1)[j]: string;
        table[k, ..] = tmp;
        k += 1;
      }
    }
    if k > tabLen
    {
      break;
    }
  }
  var file = open("../data/chapelBench.csv", iomode.cw);
  var _channel = file.writer(); //iostyle.text()
  _channel.write(table);
  _channel.close();
  file.close();
  return;
}

proc main()
{
  runAllKernelBenchmarks(real(32));
}