gauss3.dart

import 'dart:math';
import 'dart:typeddata';

main() {
  print("Generating 10 40-length vector");
  GaussData g = new GaussData.random(4);
  InputData dataSmall = new InputData.random(4,4);  
  perfList(g, dataSmall);
  perfTyped(g, dataSmall);
  perfSIMD(g, dataSmall);
  print("Generating 100000 40-length vector");  
  InputData dataLarge = new InputData.random(100000,40);
  g = new GaussData.random(40);  
  perfList(g,dataLarge);    
  perfTyped(g,dataLarge);
  perfSIMD(g,dataLarge);
  print("Done.");
}

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

class GaussData {
  List<double> means;
  List<double> variances;
  int dimension;
  
  GaussData(this.means, this.variances){
    this.dimension = means.length;
  }
  
  GaussData.random(this.dimension) {
    means = new List(dimension);
    variances = new List(dimension);    
    for(int i = 0 ; i<dimension; i++){
      means[i]=_random();
      variances[i] = _random();
    }
  }
}

class InputData {
  List<List<double>> data;
  int size;
  int dimension;
  
  InputData.random(this.size, this.dimension){
    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();
      }
    }    
  }
  
  List<Float32List> getTyped() {
    List<Float32List> res = new List<Float32List>(size);
    for (int i = 0; i < size; i++) {
      Float32List lst = new Float32List(dimension);
      for(int j = 0; j<dimension; j++) {
        lst[j] = data[i][j];
      }
      res[i]=lst;
    }
    return res;
  }
  
  List<Float32x4List> getSIMD() {
    List<Float32x4List> res = new List<Float32x4List>(size);
    for (int i = 0; i < size; i++) {
      res[i] = toFloat32x4List(data[i]);
    }
    return res;
  } 
}

void calculate(DiagonalGaussian gaussian, List data) {
  Stopwatch sw = new Stopwatch()..start();
  double tot = 0.0;    
  for (int i = 0; i<data.length; ++i) {
    tot= tot+ gaussian.logLikelihood(data[i]);
  }
  print("log-likelihood=$tot");      
} 

void perfList(GaussData gauss, InputData input) {
    var d = new ListDiagonalGaussian(gauss.means, gauss.variances);
    Stopwatch sw = new Stopwatch()..start();
    double tot = 0.0;    
    for (int i = 0; i<input.size; ++i) {
      tot= tot+ d.logLikelihood(input.data[i]);
    }
    print("log-likelihood=$tot"); 
    print("List Elapsed: ${sw.elapsedMilliseconds}");  
}

void perfTyped(GaussData gauss, InputData input) {
    var d = new TypedDataDiagonalGaussian(toFloat32List(gauss.means),toFloat32List(gauss.variances));
    List<Float32List> data = input.getTyped();
    Stopwatch sw = new Stopwatch()..start();    
    double tot = 0.0;  
    for (int i = 0; i<input.size; ++i) {
      tot= tot+ d.logLikelihood(data[i]);
    }
    print("log-likelihood=$tot"); 
    print("Typed Elapsed: ${sw.elapsedMilliseconds}");  
}

void perfSIMD(GaussData gauss, InputData input) {
    List<Float32x4List> data = input.getSIMD();
    var d = new SimdMultivariateDiagonalGaussian(toFloat32x4List(gauss.means), toFloat32x4List(gauss.variances));
    Stopwatch sw = new Stopwatch()..start();
    double tot = 0.0;    
    for (int i = 0; i<data.length; ++i) {
      tot= tot+ d.logLikelihood(data[i]);
    }
    print("log-likelihood=$tot"); 
    print("SIMD Elapsed: ${sw.elapsedMilliseconds}");  
}

Float32List toFloat32List(List<double> input) {
  Float32List lst = new Float32List(input.length);
  for(int i =0; i< input.length; i++) {
    lst[i]=input[i];
  }
  return lst;
}

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].z;
    result[i*4+3] = lst[i].w;      
  }
  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;
}

abstract class DiagonalGaussian {
  
  double logPrecomputedDistance;  
  
  /// Calculates log likelihood of the given vector [data].
  double logLikelihood(List data);  
  
  /// dimension of this multiavriate Gaussian.
  int get dimension;    
}

/**
 * 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 ListDiagonalGaussian extends DiagonalGaussian {
  List<double> means;
  List<double> variances;
  List<double> negativeHalfPrecisions;
  
  ListDiagonalGaussian(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 res = 0.0;
    for (int i = 0; i < means.length; i++) {
        final double dif = data[i] - means[i];
        res += ((dif * dif) * negativeHalfPrecisions[i]);     
    }    
    return logPrecomputedDistance + res;
  }  
  
  /// dimension of this multiavriate Gaussian.
  int get dimension => means.length;    
 
}

class TypedDataDiagonalGaussian extends DiagonalGaussian {
  Float32List means;
  Float32List variances;
  Float32List negativeHalfPrecisions;
  
  TypedDataDiagonalGaussian(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 Float32List(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(Float32List data) {
    double res = 0.0;
    for (int i = 0; i < means.length; i++) {
        final double dif = data[i] - means[i];
        res += ((dif * dif) * negativeHalfPrecisions[i]);     
    }    
    return logPrecomputedDistance + res;
  }  
  
  /// dimension of this multiavriate Gaussian.
  int get dimension => means.length;    
 
}

class SimdMultivariateDiagonalGaussian implements DiagonalGaussian {
  
  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(means.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.z,-0.5/v.w);
    }

    // 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(this.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;   

}