<commit_before>package hex.svd;

import Jama.CholeskyDecomposition;
import Jama.Matrix;
import Jama.QRDecomposition;
import Jama.SingularValueDecomposition;
import hex.*;
import hex.gram.Gram;
import hex.gram.Gram.GramTask;
import hex.schemas.ModelBuilderSchema;
import hex.schemas.SVDV99;
import hex.svd.SVDModel.SVDParameters;
import hex.util.LinearAlgebraUtils.*;
import water.*;
import water.fvec.Chunk;
import water.fvec.Frame;
import water.fvec.NewChunk;
import water.fvec.Vec;
import water.util.ArrayUtils;
import water.util.Log;

import java.util.Arrays;

/**
 * Singular Value Decomposition
 * <a href = "http://www.cs.yale.edu/homes/el327/datamining2013aFiles/07_singular_value_decomposition.pdf">SVD via Power Method Algorithm</a>
 * <a href = "https://www.cs.cmu.edu/~venkatg/teaching/CStheory-infoage/book-chapter-4.pdf">Proof of Convergence for Power Method</a>
 * <a href = "http://arxiv.org/pdf/0909.4061.pdf">Probabilistic Algorithms for Matrix Approximation</a>
 * @author anqi_fu
 */
public class SVD extends ModelBuilder<SVDModel,SVDModel.SVDParameters,SVDModel.SVDOutput> {
  // Convergence tolerance
  private final double TOLERANCE = 1e-6;    // Cutoff for estimation error of right singular vector

  // Maximum number of columns when categoricals expanded
  private final int MAX_COLS_EXPANDED = 5000;

  // Number of columns in training set (p)
  private transient int _ncolExp;    // With categoricals expanded into 0/1 indicator cols

  @Override public ModelBuilderSchema schema() {
    return new SVDV99();
  }

  @Override public Job<SVDModel> trainModelImpl(long work, boolean restartTimer) {
    return start(new SVDDriver(), work, restartTimer);
  }

  @Override
  public long progressUnits() {
    return _parms._nv+1;
  }


  @Override public ModelCategory[] can_build() {
    return new ModelCategory[]{ ModelCategory.DimReduction };
  }

  @Override public BuilderVisibility builderVisibility() { return BuilderVisibility.Experimental; };

  // Called from an http request
  public SVD(SVDModel.SVDParameters parms) {
    super("SVD", parms);
    init(false);
  }

  @Override public void init(boolean expensive) {
    super.init(expensive);
    // if (_parms._u_key == null) _parms._u_key = Key.make("SVDUMatrix_" + Key.rand());
    if (_parms._u_name == null || _parms._u_name.length() == 0)
      _parms._u_name = "SVDUMatrix_" + Key.rand();
    if (_parms._max_iterations < 1)
      error("_max_iterations", "max_iterations must be at least 1");

    if(_train == null) return;
    _ncolExp = _train.numColsExp(_parms._use_all_factor_levels, false);
    if (_ncolExp > MAX_COLS_EXPANDED)
      warn("_train", "_train has " + _ncolExp + " columns when categoricals are expanded. Algorithm may be slow.");

    if(_parms._nv < 1 || _parms._nv > _ncolExp)
      error("_nv", "Number of right singular values must be between 1 and " + _ncolExp);
  }

  // Compute ivv_sum - vec * vec' for symmetric array ivv_sum
  public static double[][] updateIVVSum(double[][] ivv_sum, double[] vec) {
    double diff;
    for(int i = 0; i < vec.length; i++) {
      for(int j = 0; j < i; j++) {
        diff = ivv_sum[i][j] - vec[i] * vec[j];
        ivv_sum[i][j] = ivv_sum[j][i] = diff;
      }
      ivv_sum[i][i] -= vec[i] * vec[i];
    }
    return ivv_sum;
  }

  class SVDDriver extends H2O.H2OCountedCompleter<SVDDriver> {

    // private double[] powerLoop(Gram gram) { return powerLoop(gram, ArrayUtils.gaussianVector(gram.fullN())); }
    private double[] powerLoop(Gram gram, long seed) { return powerLoop(gram, ArrayUtils.gaussianVector(gram.fullN(), seed)); }
    private double[] powerLoop(Gram gram, double[] vinit) {
      // TODO: What happens if Gram matrix is essentially zero? Numerical inaccuracies in PUBDEV-1161.
      assert vinit.length == gram.fullN();

      // Set initial value v_0 to standard normal distribution
      int iters = 0;
      double err = 2 * TOLERANCE;
      double[] v = vinit.clone();
      double[] vnew = new double[v.length];

      // Update v_i <- (A'Av_{i-1})/||A'Av_{i-1}|| where A'A = Gram matrix of training frame
      while(iters < _parms._max_iterations && err > TOLERANCE) {
        // Compute x_i <- A'Av_{i-1} and ||x_i||
        gram.mul(v, vnew);
        double norm = ArrayUtils.l2norm(vnew);

        double diff; err = 0;
        for (int i = 0; i < v.length; i++) {
          vnew[i] /= norm;        // Compute singular vector v_i = x_i/||x_i||
          diff = v[i] - vnew[i];  // Save error ||v_i - v_{i-1}||
          err += diff * diff;
          v[i] = vnew[i];         // Update v_i for next iteration
        }
        err = Math.sqrt(err);
        iters++;    // TODO: Should output vector of final iterations for each k
      }
      return v;
    }

    private double computeSigmaU(DataInfo dinfo, SVDModel model, int k, double[][] ivv_sum, Vec[] uvecs) {
      double[] ivv_vk = ArrayUtils.multArrVec(ivv_sum, model._output._v[k]);
      CalcSigmaU ctsk = new CalcSigmaU(self(), dinfo, ivv_vk).doAll(1, dinfo._adaptedFrame);
      model._output._d[k] = ctsk._sval;
      assert ctsk._nobs == model._output._nobs : "Processed " + ctsk._nobs + " rows but expected " + model._output._nobs;    // Check same number of skipped rows as Gram
      Frame tmp = ctsk.outputFrame();
      uvecs[k] = tmp.vec(0);   // Save output column of U
      tmp.unlock(self());
      return model._output._d[k];
    }

    // Algorithm 4.4: Randomized subspace iteration from Halk et al (http://arxiv.org/pdf/0909.4061.pdf)
    private Frame randSubIter(DataInfo dinfo, SVDModel model, int iters, long seed) {
      DataInfo yinfo = null;
      Key yinfo_key = Key.make(), ykey = Key.make();
      Frame yinit = null, ybig = null, qfrm = null, ayqfrm = null;
      final int ncolA = dinfo._adaptedFrame.numCols();

      try {
        // 1) Initialize Y = AG where G ~ N(0,1) and compute Y = QR factorization
        double[][] gt = ArrayUtils.gaussianArray(_parms._nv, _ncolExp, seed);
        RandSubInit rtsk = new RandSubInit(self(), dinfo, gt);
        rtsk.doAll(_parms._nv, dinfo._adaptedFrame);
        yinit = rtsk.outputFrame(ykey, null, null);

        // Calculate Cholesky of Gram to get R' = L matrix
        yinfo = new DataInfo(yinfo_key, yinit, null, true, DataInfo.TransformType.NONE, true, false, false);
        DKV.put(yinfo._key, yinfo);
        GramTask gtsk = new GramTask(self(), yinfo);  // Gram is Y'Y/n where n = nrow(Y)
        gtsk.doAll(yinfo._adaptedFrame);
        // Gram.Cholesky chol = gtsk._gram.cholesky(null);   // If Y'Y = LL' Cholesky, then R = L'
        Matrix ygram = new Matrix(gtsk._gram.getXX());
        CholeskyDecomposition chol = new CholeskyDecomposition(ygram);
        System.out.println(ArrayUtils.pprint(chol.getL().getArray()));

        // Make input frame [A,Q,Y] where A = read-only training data, Y = AQ, Q from Y = QR factorization
        ayqfrm = new Frame(dinfo._adaptedFrame);
        ayqfrm.add(yinit);
        for (int i = 0; i < _parms._nv; i++)
          ayqfrm.add("qcol_" + i, ayqfrm.anyVec().makeZero());

        // Get Q from Y = QR giving R' = L
        double[] normSubY = new double[_parms._nv];
        double[] normMulY = new double[_parms._nv];
        Arrays.fill(normMulY, 1.0);
        Frame yqfrm = ayqfrm.subframe(ncolA, ayqfrm.numCols());   // Pass in [Y,Q]
        QRfromChol qrtsk = new QRfromChol(chol, gtsk._nobs, _parms._nv, _parms._nv, normSubY, normMulY);
        qrtsk.doAll(yqfrm);

        for (int q = 0; q < iters; q++) {
          // 2) Form \tilde{Y}_j = A'Q_{j-1} and compute \tilde{Y}_j = \tilde{Q}_j \tilde{R}_j factorization
          Frame aqfrm = ayqfrm.subframe(0, ncolA);
          aqfrm.add(ayqfrm.subframe(ncolA+_parms._nv, ayqfrm.numCols()));
          SMulTask stsk = new SMulTask(ncolA, _ncolExp, model._output._ncats, _parms._nv, model._output._normSub, model._output._normMul, model._output._catOffsets, _parms._use_all_factor_levels);
          stsk.doAll(aqfrm);    // Pass in [A,Q]

          Matrix ysmall = new Matrix(stsk._atq);
          QRDecomposition ysmall_qr = new QRDecomposition(ysmall);
          double[][] ysmall_q = ysmall_qr.getQ().getArray();

          // 3) Form Y_j = A\tilde{Q}_j and compute Y_j = Q_jR_j factorization
          Frame ayfrm = ayqfrm.subframe(0,ncolA+_parms._nv);
          BMulInPlaceTask tsk = new BMulInPlaceTask(ncolA, ArrayUtils.transpose(ysmall_q));
          tsk.doAll(ayfrm);     // TODO: Can we reuse [A,Q] frame and write Y_j into old Q_{j-1}?
          ybig = ayfrm.subframe(ncolA, ayfrm.numCols());
          ybig = new Frame(ykey, ybig.names(), ybig.vecs());
          DKV.put(ybig);

          // Calculate Cholesky of Gram to get R' = L matrix
          yinfo = new DataInfo(yinfo_key, ybig, null, true, DataInfo.TransformType.NONE, true, false, false);
          DKV.put(yinfo._key, yinfo);
          gtsk = new GramTask(self(), yinfo);  // Gram is Y'Y/n where n = nrow(Y)
          gtsk.doAll(yinfo._adaptedFrame);
          // Gram.Cholesky chol = gtsk._gram.cholesky(null);   // If Y'Y = LL' Cholesky, then R = L'
          ygram = new Matrix(gtsk._gram.getXX());
          chol = new CholeskyDecomposition(ygram);

          yqfrm = ayqfrm.subframe(ncolA, ayqfrm.numCols());
          qrtsk = new QRfromChol(chol, gtsk._nobs, _parms._nv, _parms._nv, normSubY, normMulY);
          qrtsk.doAll(yqfrm);   // Pass in [Y,Q]
        }

        // 4) Extract and save final Q_j from [A,Q] frame
        qfrm = ayqfrm.extractFrame(ncolA+_parms._nv, ayqfrm.numCols());
        qfrm = new Frame(Key.make(), qfrm.names(), qfrm.vecs());
        DKV.put(qfrm);
      } catch( Throwable t ) {
        Job thisJob = DKV.getGet(_key);
        if (thisJob._state == JobState.CANCELLED) {
          Log.info("Job cancelled by user.");
        } else {
          t.printStackTrace();
          failed(t);
          throw t;
        }
      } finally {
        if( yinfo != null ) yinfo.remove();
        if( yinit != null ) yinit.delete();
        if( ybig != null ) ybig.delete();
        if( ayqfrm != null ) {
          for(int i = ncolA; i < ncolA+_parms._nv; i++)
            ayqfrm.vec(i).remove();
        }
      }
      System.out.println(qfrm.toString());
      return qfrm;
    }

    // Algorithm 5.1: Direct SVD from Halko et al (http://arxiv.org/pdf/0909.4061.pdf)
    private Frame directSVD(DataInfo dinfo, Frame qfrm, SVDModel model) {
      DataInfo qinfo = null;
      Frame u = null;
      final int ncolA = dinfo._adaptedFrame.numCols();

      try {
        // 0) Make input frame [A,Q], where A = read-only training data, Q = matrix from randomized subspace iteration
        Vec[] vecs = new Vec[ncolA + _parms._nv];
        for (int i = 0; i < ncolA; i++) vecs[i] = dinfo._adaptedFrame.vec(i);
        for (int i = 0; i < _parms._nv; i++) vecs[ncolA + i] = qfrm.vec(i);
        Frame aqfrm = new Frame(vecs);

        // 1) Form the matrix B = Q'A = (A'Q)'
        SMulTask stsk = new SMulTask(ncolA, _ncolExp, model._output._ncats, _parms._nv, model._output._normSub, model._output._normMul, model._output._catOffsets, _parms._use_all_factor_levels);
        stsk.doAll(aqfrm);
        double[][] qta = ArrayUtils.transpose(stsk._atq);

        // 2) Compute SVD of small matrix B = WDV'
        Matrix qtaJ = new Matrix(qta);
        SingularValueDecomposition svdJ = qtaJ.svd();

        // 3) Form orthonormal matrix U = QW
        if (!_parms._only_v && _parms._keep_u) {<commit_msg>Forgot to remove debug print statements<commit_after>package hex.svd;

import Jama.CholeskyDecomposition;
import Jama.Matrix;
import Jama.QRDecomposition;
import Jama.SingularValueDecomposition;
import hex.*;
import hex.gram.Gram;
import hex.gram.Gram.GramTask;
import hex.schemas.ModelBuilderSchema;
import hex.schemas.SVDV99;
import hex.svd.SVDModel.SVDParameters;
import hex.util.LinearAlgebraUtils.*;
import water.*;
import water.fvec.Chunk;
import water.fvec.Frame;
import water.fvec.NewChunk;
import water.fvec.Vec;
import water.util.ArrayUtils;
import water.util.Log;

import java.util.Arrays;

/**
 * Singular Value Decomposition
 * <a href = "http://www.cs.yale.edu/homes/el327/datamining2013aFiles/07_singular_value_decomposition.pdf">SVD via Power Method Algorithm</a>
 * <a href = "https://www.cs.cmu.edu/~venkatg/teaching/CStheory-infoage/book-chapter-4.pdf">Proof of Convergence for Power Method</a>
 * <a href = "http://arxiv.org/pdf/0909.4061.pdf">Probabilistic Algorithms for Matrix Approximation</a>
 * @author anqi_fu
 */
public class SVD extends ModelBuilder<SVDModel,SVDModel.SVDParameters,SVDModel.SVDOutput> {
  // Convergence tolerance
  private final double TOLERANCE = 1e-6;    // Cutoff for estimation error of right singular vector

  // Maximum number of columns when categoricals expanded
  private final int MAX_COLS_EXPANDED = 5000;

  // Number of columns in training set (p)
  private transient int _ncolExp;    // With categoricals expanded into 0/1 indicator cols

  @Override public ModelBuilderSchema schema() {
    return new SVDV99();
  }

  @Override public Job<SVDModel> trainModelImpl(long work, boolean restartTimer) {
    return start(new SVDDriver(), work, restartTimer);
  }

  @Override
  public long progressUnits() {
    return _parms._nv+1;
  }


  @Override public ModelCategory[] can_build() {
    return new ModelCategory[]{ ModelCategory.DimReduction };
  }

  @Override public BuilderVisibility builderVisibility() { return BuilderVisibility.Experimental; };

  // Called from an http request
  public SVD(SVDModel.SVDParameters parms) {
    super("SVD", parms);
    init(false);
  }

  @Override public void init(boolean expensive) {
    super.init(expensive);
    // if (_parms._u_key == null) _parms._u_key = Key.make("SVDUMatrix_" + Key.rand());
    if (_parms._u_name == null || _parms._u_name.length() == 0)
      _parms._u_name = "SVDUMatrix_" + Key.rand();
    if (_parms._max_iterations < 1)
      error("_max_iterations", "max_iterations must be at least 1");

    if(_train == null) return;
    _ncolExp = _train.numColsExp(_parms._use_all_factor_levels, false);
    if (_ncolExp > MAX_COLS_EXPANDED)
      warn("_train", "_train has " + _ncolExp + " columns when categoricals are expanded. Algorithm may be slow.");

    if(_parms._nv < 1 || _parms._nv > _ncolExp)
      error("_nv", "Number of right singular values must be between 1 and " + _ncolExp);
  }

  // Compute ivv_sum - vec * vec' for symmetric array ivv_sum
  public static double[][] updateIVVSum(double[][] ivv_sum, double[] vec) {
    double diff;
    for(int i = 0; i < vec.length; i++) {
      for(int j = 0; j < i; j++) {
        diff = ivv_sum[i][j] - vec[i] * vec[j];
        ivv_sum[i][j] = ivv_sum[j][i] = diff;
      }
      ivv_sum[i][i] -= vec[i] * vec[i];
    }
    return ivv_sum;
  }

  class SVDDriver extends H2O.H2OCountedCompleter<SVDDriver> {

    // private double[] powerLoop(Gram gram) { return powerLoop(gram, ArrayUtils.gaussianVector(gram.fullN())); }
    private double[] powerLoop(Gram gram, long seed) { return powerLoop(gram, ArrayUtils.gaussianVector(gram.fullN(), seed)); }
    private double[] powerLoop(Gram gram, double[] vinit) {
      // TODO: What happens if Gram matrix is essentially zero? Numerical inaccuracies in PUBDEV-1161.
      assert vinit.length == gram.fullN();

      // Set initial value v_0 to standard normal distribution
      int iters = 0;
      double err = 2 * TOLERANCE;
      double[] v = vinit.clone();
      double[] vnew = new double[v.length];

      // Update v_i <- (A'Av_{i-1})/||A'Av_{i-1}|| where A'A = Gram matrix of training frame
      while(iters < _parms._max_iterations && err > TOLERANCE) {
        // Compute x_i <- A'Av_{i-1} and ||x_i||
        gram.mul(v, vnew);
        double norm = ArrayUtils.l2norm(vnew);

        double diff; err = 0;
        for (int i = 0; i < v.length; i++) {
          vnew[i] /= norm;        // Compute singular vector v_i = x_i/||x_i||
          diff = v[i] - vnew[i];  // Save error ||v_i - v_{i-1}||
          err += diff * diff;
          v[i] = vnew[i];         // Update v_i for next iteration
        }
        err = Math.sqrt(err);
        iters++;    // TODO: Should output vector of final iterations for each k
      }
      return v;
    }

    private double computeSigmaU(DataInfo dinfo, SVDModel model, int k, double[][] ivv_sum, Vec[] uvecs) {
      double[] ivv_vk = ArrayUtils.multArrVec(ivv_sum, model._output._v[k]);
      CalcSigmaU ctsk = new CalcSigmaU(self(), dinfo, ivv_vk).doAll(1, dinfo._adaptedFrame);
      model._output._d[k] = ctsk._sval;
      assert ctsk._nobs == model._output._nobs : "Processed " + ctsk._nobs + " rows but expected " + model._output._nobs;    // Check same number of skipped rows as Gram
      Frame tmp = ctsk.outputFrame();
      uvecs[k] = tmp.vec(0);   // Save output column of U
      tmp.unlock(self());
      return model._output._d[k];
    }

    // Algorithm 4.4: Randomized subspace iteration from Halk et al (http://arxiv.org/pdf/0909.4061.pdf)
    private Frame randSubIter(DataInfo dinfo, SVDModel model, int iters, long seed) {
      DataInfo yinfo = null;
      Key yinfo_key = Key.make(), ykey = Key.make();
      Frame yinit = null, ybig = null, qfrm = null, ayqfrm = null;
      final int ncolA = dinfo._adaptedFrame.numCols();

      try {
        // 1) Initialize Y = AG where G ~ N(0,1) and compute Y = QR factorization
        double[][] gt = ArrayUtils.gaussianArray(_parms._nv, _ncolExp, seed);
        RandSubInit rtsk = new RandSubInit(self(), dinfo, gt);
        rtsk.doAll(_parms._nv, dinfo._adaptedFrame);
        yinit = rtsk.outputFrame(ykey, null, null);

        // Calculate Cholesky of Gram to get R' = L matrix
        yinfo = new DataInfo(yinfo_key, yinit, null, true, DataInfo.TransformType.NONE, true, false, false);
        DKV.put(yinfo._key, yinfo);
        GramTask gtsk = new GramTask(self(), yinfo);  // Gram is Y'Y/n where n = nrow(Y)
        gtsk.doAll(yinfo._adaptedFrame);
        // Gram.Cholesky chol = gtsk._gram.cholesky(null);   // If Y'Y = LL' Cholesky, then R = L'
        Matrix ygram = new Matrix(gtsk._gram.getXX());
        CholeskyDecomposition chol = new CholeskyDecomposition(ygram);

        // Make input frame [A,Q,Y] where A = read-only training data, Y = AQ, Q from Y = QR factorization
        ayqfrm = new Frame(dinfo._adaptedFrame);
        ayqfrm.add(yinit);
        for (int i = 0; i < _parms._nv; i++)
          ayqfrm.add("qcol_" + i, ayqfrm.anyVec().makeZero());

        // Get Q from Y = QR giving R' = L
        double[] normSubY = new double[_parms._nv];
        double[] normMulY = new double[_parms._nv];
        Arrays.fill(normMulY, 1.0);
        Frame yqfrm = ayqfrm.subframe(ncolA, ayqfrm.numCols());   // Pass in [Y,Q]
        QRfromChol qrtsk = new QRfromChol(chol, gtsk._nobs, _parms._nv, _parms._nv, normSubY, normMulY);
        qrtsk.doAll(yqfrm);

        for (int q = 0; q < iters; q++) {
          // 2) Form \tilde{Y}_j = A'Q_{j-1} and compute \tilde{Y}_j = \tilde{Q}_j \tilde{R}_j factorization
          Frame aqfrm = ayqfrm.subframe(0, ncolA);
          aqfrm.add(ayqfrm.subframe(ncolA+_parms._nv, ayqfrm.numCols()));
          SMulTask stsk = new SMulTask(ncolA, _ncolExp, model._output._ncats, _parms._nv, model._output._normSub, model._output._normMul, model._output._catOffsets, _parms._use_all_factor_levels);
          stsk.doAll(aqfrm);    // Pass in [A,Q]

          Matrix ysmall = new Matrix(stsk._atq);
          QRDecomposition ysmall_qr = new QRDecomposition(ysmall);
          double[][] ysmall_q = ysmall_qr.getQ().getArray();

          // 3) Form Y_j = A\tilde{Q}_j and compute Y_j = Q_jR_j factorization
          Frame ayfrm = ayqfrm.subframe(0,ncolA+_parms._nv);
          BMulInPlaceTask tsk = new BMulInPlaceTask(ncolA, ArrayUtils.transpose(ysmall_q));
          tsk.doAll(ayfrm);     // TODO: Can we reuse [A,Q] frame and write Y_j into old Q_{j-1}?
          ybig = ayfrm.subframe(ncolA, ayfrm.numCols());
          ybig = new Frame(ykey, ybig.names(), ybig.vecs());
          DKV.put(ybig);

          // Calculate Cholesky of Gram to get R' = L matrix
          yinfo = new DataInfo(yinfo_key, ybig, null, true, DataInfo.TransformType.NONE, true, false, false);
          DKV.put(yinfo._key, yinfo);
          gtsk = new GramTask(self(), yinfo);  // Gram is Y'Y/n where n = nrow(Y)
          gtsk.doAll(yinfo._adaptedFrame);
          // Gram.Cholesky chol = gtsk._gram.cholesky(null);   // If Y'Y = LL' Cholesky, then R = L'
          ygram = new Matrix(gtsk._gram.getXX());
          chol = new CholeskyDecomposition(ygram);

          yqfrm = ayqfrm.subframe(ncolA, ayqfrm.numCols());
          qrtsk = new QRfromChol(chol, gtsk._nobs, _parms._nv, _parms._nv, normSubY, normMulY);
          qrtsk.doAll(yqfrm);   // Pass in [Y,Q]
        }

        // 4) Extract and save final Q_j from [A,Q] frame
        qfrm = ayqfrm.extractFrame(ncolA+_parms._nv, ayqfrm.numCols());
        qfrm = new Frame(Key.make(), qfrm.names(), qfrm.vecs());
        DKV.put(qfrm);
      } catch( Throwable t ) {
        Job thisJob = DKV.getGet(_key);
        if (thisJob._state == JobState.CANCELLED) {
          Log.info("Job cancelled by user.");
        } else {
          t.printStackTrace();
          failed(t);
          throw t;
        }
      } finally {
        if( yinfo != null ) yinfo.remove();
        if( yinit != null ) yinit.delete();
        if( ybig != null ) ybig.delete();
        if( ayqfrm != null ) {
          for(int i = ncolA; i < ncolA+_parms._nv; i++)
            ayqfrm.vec(i).remove();
        }
      }
      return qfrm;
    }

    // Algorithm 5.1: Direct SVD from Halko et al (http://arxiv.org/pdf/0909.4061.pdf)
    private Frame directSVD(DataInfo dinfo, Frame qfrm, SVDModel model) {
      DataInfo qinfo = null;
      Frame u = null;
      final int ncolA = dinfo._adaptedFrame.numCols();

      try {
        // 0) Make input frame [A,Q], where A = read-only training data, Q = matrix from randomized subspace iteration
        Vec[] vecs = new Vec[ncolA + _parms._nv];
        for (int i = 0; i < ncolA; i++) vecs[i] = dinfo._adaptedFrame.vec(i);
        for (int i = 0; i < _parms._nv; i++) vecs[ncolA + i] = qfrm.vec(i);
        Frame aqfrm = new Frame(vecs);

        // 1) Form the matrix B = Q'A = (A'Q)'
        SMulTask stsk = new SMulTask(ncolA, _ncolExp, model._output._ncats, _parms._nv, model._output._normSub, model._output._normMul, model._output._catOffsets, _parms._use_all_factor_levels);
        stsk.doAll(aqfrm);
        double[][] qta = ArrayUtils.transpose(stsk._atq);

        // 2) Compute SVD of small matrix B = WDV'
        Matrix qtaJ = new Matrix(qta);
        SingularValueDecomposition svdJ = qtaJ.svd();

        // 3) Form orthonormal matrix U = QW
        if (!_parms._only_v && _parms._keep_u) {
{'fraction_non_alphanumeric': 0.1330938403574871, 'fraction_numerical': 0.008411460615088058, 'mean_word_length': 2.7631058358061327, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 27, 'https://': 2, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2752528', 'n_tokens_mistral': 8342, 'n_tokens_neox': 7945, 'n_words': 1981}