import 'dart:io';
import 'dart:math';


/**
 * See LinearRegression2() for better singular matrix handler
 *
 */

List list(int len, [func]){
  if( func == null ) return new List(len);
  if( func is Function ) return new List.generate(len, func);
  return new List.filled(len, func);
}
List<List> matrix(int height, int width, [func]){
  if( func is Function ) return new List.generate( height, (i){
    return new List.generate(width, (j) => func(i,j));
  }, growable: false );
  return new List.generate( height, (_) => list(width, func), growable: false );
}


// O( 4/3 * N^3 )
// mat is manipulated in-place, clone if it is not desired.
QRDecomposition(List<List<double> > mat){
  var h = mat.length;
  var w = mat[0].length;
  var mi = min(h,w);

  for(var i=0; i<mi; i++){
    var p = mat[i][i];

    // norm of column-i
    var s = 0.0;
    for(var j=i; j<h; j++){
      var x = mat[j][i];
      s += x*x;
    }

    // theta
    var t = sqrt(s);
    if( p.sign == 1 ) t = -t;

    // update columns i+1 to end
    var b0 = (s - p*t);
    for(var j=i+1; j<w; j++){

      var b1 = -t * mat[i][j];
      for(var k=i; k<h; k++) b1 += mat[k][i] * mat[k][j];

      var b = b1/b0;
      mat[i][j] += b * t;
      for(var k=i; k<h; k++) mat[k][j] -= b * mat[k][i];
    }

    // override columns-i
    mat[i][i] = t;
    for(var j=i+1; j<h; j++) mat[j][i] = 0.0;
  }
}

// QRDecomposition2(mat) = QRDecomposition( transpose(mat) )
QRDecomposition2(List<List<double> > mat){
  var h = mat.length;
  var w = mat[0].length;

  assert( w > h );

  dotFrom(i,a,b){
    var s = 0.0;
    for(var j=i; j<w; j++) s += a[j] * b[j];
    return s;
  }

  // a[i..w] += c * b[i..w]
  rowAddFrom(i, a, c, b){
    for(var j=i; j<w; j++){
      a[j] += c * b[j];
    }
  }

  // Householder transform
  hhTransform(i, b0, x, y){
    var b1 = dotFrom(i,x,y);
    rowAddFrom(i, y, -b1*b0, x);
  }

  for(var i=0; i<h; i++){
    // diagonal element (pivot)
    var p = mat[i][i];

    // norm of row-i s >= 0
    var s = dotFrom(i+1, mat[i], mat[i]);

    // theta
    var t = sqrt(p*p+s);
    // increase p - t

    // avoid subtraction near equal number,
    if( p.sign == 1 ) t = -t;
    p = mat[i][i] -= t;


    var b0 = 2/(p*p+s);
    for(var j=i+1; j<h; j++){
      hhTransform(i, b0, mat[i], mat[j]);
    }

    // override columns-i
    mat[i][i] = t;
    // this is not necessary for least square problem
    mat[i].fillRange(i+1, w, 0.0);
  }
}