gistfile1.dart

import 'dart:math';
import 'dart:scalarlist';

main() {
  print("Generating 10 40-length vector");
  List<List<double>> dataSmall = getData(10,40);  
  perfGaussian(dataSmall);
  perfSimdGaussian(dataSmall);  
  print("Generating 100000 40-length vector");  
  List<List<double>> dataLarge = getData(100000,40);
  perfGaussian(dataLarge);    
  perfSimdGaussian(dataLarge);    
  print("Done.");
}

void perfGaussian(List<List<double>> data) {
    Random random = new Random();
    var d = new MultivariateDiagonalGaussian(data[0], data[1]);
    Stopwatch sw = new Stopwatch()..start();
    double tot = 0.0;    
    for (int i = 0; i<data.length; ++i) {
      tot= tot+ d.logLikelihood(data[i]);
    }
    print("Elapsed: ${sw.elapsedMilliseconds} log-likelihood=$tot");  
}

List<List<double>> getData(int size, int dimension) {
  List<List<double>> data = new List<List<double>>(size);
  for (int i = 0; i < data.length; i++) {
    data[i] = new List<double>(dimension);
    for (int j = 0; j < dimension; j++) {
      data[i][j] = _random();
    }
  }
  return data;
}

Random random = new Random(0xcafebabe);
_random() {
    return random.nextInt(10) / 10 + 0.1;
}

void perfSimdGaussian(List<List<double>> _d) {
    List<Float32x4List> data = new List<Float32x4List>(_d.length);
    for (int i = 0; i < data.length; i++) {
      data[i] = toFloat32x4List(_d[i]);
    }    
    var d = new SimdMultivariateDiagonalGaussian(data[0], data[1]);
    Stopwatch sw = new Stopwatch()..start();
    double tot = 0.0;    
    for (int i = 0; i<data.length; ++i) {
      tot= tot+ d.logLikelihood(data[i]);
    }
    print("SIMD elapsed: ${sw.elapsedMilliseconds} log-likelihood=$tot");  
}

List<double> toDoubleList(Float32x4List lst) {
  List<double> result = new List(lst.length*4);
  for(int i =0; i<lst.length;i++){
    result[i*4] = lst[i].x;
    result[i*4+1] = lst[i].y;
    result[i*4+2] = lst[i].w;
    result[i*4+3] = lst[i].z;      
  }
  return result;
}

Float32x4List toFloat32x4List(List<double> list) {
  int dsize = list.length~/4;
  Float32x4List res = new Float32x4List(dsize);
  for (int j = 0; j < dsize; j++) {        
    res[j] = new Float32x4(
        list[j*4],
        list[j*4+1],
        list[j*4+2],
        list[j*4+3]);        
  }  
  return res;
}

/**
 * Represents a Multivariate Gaussian function with diagonal covariance matrix.
 * Normally a multivariate Gauss function with diagonal covariance matrix is:
 * f(x) = MUL[d=1..D] (1/(sqrt(2 *PI * var[d])) exp(-0.5 * (x[d] - mean[d])^2 / var[d]))
 * 
 * But we are interested in log likelihoods. Calculating linear likelihood values is expensive because of sqrt and exp
 * operations. And because this operation will likely to be done millions of times, it needs to be very
 * effective. log likelihood is much more efficient to calculate. Therefore it becomes
 * log f(x) = -0.5*SUM[d=1..D] ( log(2*PI) + log(var[d]) +  (x[d] - mean[d])^2 / var[d]) )
 * 
 * This can be effectively calculated if right side is split as follows
 * log f(x) = -0.5*SUM[d=1..D] ( log(2*PI) + log(var[d]) ) - 0.5* SUM[1..D]((xd - mean[d])^2 / var[d])
 * 
 * Here first part can be pre-computed as:
 * C = -0.5*log(2*PI)*D -0.5 SUM[d=1..D](log(var[d]))
 * 
 * For making remaining part faster, we pre-compute -0.5*1/var[d] values as well and store them as negative half precision
 * negativeHalfPrecisions[d]. So result log likelihood computation becomes:
 * log f(x) = C - SUM[d=1..D]((xd - mean[d])^2 * negativeHalfPrecisions[d])
 */
class MultivariateDiagonalGaussian {

  List<double> means;
  List<double> variances;
  List<double> negativeHalfPrecisions;
  double logPrecomputedDistance;

    /**
     * Constructs a Diagonal Covariance Matrix Multivariate Gaussian with given mean and variance vector.
     */
    MultivariateDiagonalGaussian(this.means, this.variances) {

        // instead of using [-0.5 * 1/var[d]] during likelihood calculation we pre-compute the values.
        // This saves 1 mul 1 div operation.
        negativeHalfPrecisions = new List<double>(variances.length);
        for (int i = 0; i < negativeHalfPrecisions.length; i++) {
            negativeHalfPrecisions[i] = -0.5 / variances[i];
        }

        // calculate the precomputed distance.
        // -0.5*SUM[d=1..D] ( log(2*PI) + log(var[d]) ) = -0.5*log(2*PI)*D -0.5 SUM[d=1..D](log(var[d]))
        double val = -0.5 * log(2 * PI) * variances.length;
        for (double variance in variances) {
            val -= (0.5 * log(variance));
        }
        logPrecomputedDistance = val;     
    }

    /// Calculates log likelihood of the given vector [data].
    double logLikelihood(List<double> data) {
        double result = 0.0;
        for (int i = 0; i < means.length; i++) {
            final double dif = data[i] - means[i];
            result += (dif * dif * negativeHalfPrecisions[i]);
        }
        return logPrecomputedDistance + result;
    }

    /// dimension of this multiavriate Gaussian.
    int get dimension => means.length;

}

class SimdMultivariateDiagonalGaussian {
  Float32x4List _means;
  Float32x4List _variances;
  Float32x4List _negativeHalfPrecisions;
  double logPrecomputedDistance;

  SimdMultivariateDiagonalGaussian(this._means, this._variances) {
    // instead of using [-0.5 * 1/var[d]] during likelihood calculation we pre-compute the values.
    _negativeHalfPrecisions = new Float32x4List(_variances.length);
    for (int i = 0; i < _negativeHalfPrecisions.length; i++) {
      Float32x4 v = _variances[i];
      _negativeHalfPrecisions[i] = new Float32x4(-0.5/v.x,-0.5/v.y,-0.5/v.w,-0.5/v.z);
    }

    // calculate the precomputed distance.
    // -0.5*SUM[d=1..D] ( log(2*PI) + log(var[d]) ) = -0.5*log(2*PI)*D -0.5 SUM[d=1..D](log(var[d]))
    List<double> variances = toDoubleList(_variances);
    double val = -0.5 * log(2 * PI) * variances.length;
    for (double variance in variances) {
        val -= (0.5 * log(variance));
    }
    logPrecomputedDistance = val;   
  }
  
  /// Calculates log likelihood of the given vector [data].
  double logLikelihood(Float32x4List data) {
    Float32x4 res = new Float32x4.zero();
    for (int i = 0; i < _means.length; i++) {
        final Float32x4 dif = data[i] - _means[i];
        res += (dif * dif * _negativeHalfPrecisions[i]);     
    }    
    return logPrecomputedDistance + res.x + res.y + res.w + res.z;
  }

  /// dimension of this multiavriate Gaussian.
  int get dimension => _means.length;

}