gt:         for (int j = 0; j < n; j++) {
context: ///////////////////////////////////////////////////////////////////////////////
  // For information as to what this class does, see the Javadoc, below.       //
  // Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006,       //
  // 2007, 2008, 2009, 2010, 2014, 2015 by Peter Spirtes, Richard Scheines, Joseph   //
  // Ramsey, and Clark Glymour.                                                //
  //                                                                           //
  // This program is free software; you can redistribute it and/or modify      //
  // it under the terms of the GNU General Public License as published by      //
  // the Free Software Foundation; either version 2 of the License, or         //
  // (at your option) any later version.                                       //
  //                                                                           //
  // This program is distributed in the hope that it will be useful,           //
  // but WITHOUT ANY WARRANTY; without even the implied warranty of            //
  // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
  // GNU General Public License for more details.                              //
  //                                                                           //
  // You should have received a copy of the GNU General Public License         //
  // along with this program; if not, write to the Free Software               //
  // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA //
  ///////////////////////////////////////////////////////////////////////////////
  
  package edu.cmu.tetrad.search;
  
  import cern.colt.matrix.DoubleMatrix1D;
  import cern.colt.matrix.DoubleMatrix2D;
  import cern.colt.matrix.impl.DenseDoubleMatrix1D;
  import cern.colt.matrix.impl.DenseDoubleMatrix2D;
  import cern.colt.matrix.linalg.Algebra;
  import cern.jet.math.Mult;
  import cern.jet.math.PlusMult;
  
  /**
   * A translation from Tibshirani's 2008 Fortran implementation of glasso.
   *
   * @author Joseph Ramsey
   */
  public class Glasso {
  
      /**
       * Dimension of matrix.
       */
      private int n = -1;
  
      /**
       * Data covariance matrix.
       */
      private DoubleMatrix2D ss;
  
      /**
       * Regularization strength parameters for each element (must be symmetric, rho(i, j) = rho(j, i).
       * False by default.
       */
      private Rho rho = new Rho() {
          public double get(int i, int j) {
              return 0;
          }
      };
  
  
      /**
       * Maximum number of iterations (no effect if ia = true).
       */
      private int maxit = 10000;
  
  
      /**
       * Approximation flag. False if exact solution, true if Meinhausen-Buhlman approximation.
       * False by default.
       */
      private boolean ia = false;
  
      /**
       * Initialization flag. false if cold start, initialize using ss. True if warm start, initialize with
       * previous solution stored in ww and wwi. False by default.
       */
      private boolean is = false;
  
      /**
       * Trace flag. True if trace information printed. False if trace information not printed.
       * False by default.
       */
      private boolean itr = false;
  
      /**
       * Diagonal penalty flag. True if diagonal is penalized. False if diagonal is not penalized.
       * False by default.
       */
      private boolean ipen = false;
  
      /**
       * Convergence threshold: Iterations stop when absolute average parameter change is less than
       * thr * avg(abs(offdiagonal(ss)). (Suggested default 1.0e-4.)
       */
      private double thr = 1.0e-4;
  
      /**
       * Return value of the algorithm.
       */
      public static class Result {
          /**
           * solution covariance matrix estimate (ia = 0)
           * (not used for ia != 0)
           */
          private DoubleMatrix2D ww;
  
          /**
           * solution inverse covariance matrix estimate (ia = 0)
           * = off-diagonal lasso coefficients (ia != 0)
           */
          private DoubleMatrix2D wwi;
  
          /**
           * number of iterations
           */
          private int niter;
  
          /**
           * average absolute parameter change at termination
           * (not used for ia != 0)
           */
          private double del;
  
          public Result(DoubleMatrix2D ww, DoubleMatrix2D wwi, int niter, double del) {
              this.ww = ww;
              this.wwi = wwi;
              this.niter = niter;
              this.del = del;
          }
  
          public DoubleMatrix2D getWw() {
              return ww;
          }
  
          public DoubleMatrix2D getWwi() {
              return wwi;
          }
  
          public int getNiter() {
              return niter;
          }
  
          public double getDel() {
              return del;
          }
      }
  
      public Glasso(DoubleMatrix2D cov) {
  
          this.n = cov.rows();
          this.ss = cov;
      }
  
      public Result search() {
          int niter = 0;
          double eps = 1.0e-7;
          int n = getN();
          DoubleMatrix2D ss = getSs();
  
          boolean approximateAlgorithm = isIa();
          boolean warmStart = isIs();
          boolean itr = isItr();
          boolean pen = isIpen();
          double thr = getThr();
  
  //        System.out.println(ss);
  
          Rho rho = getRho();
          DoubleMatrix2D ww = new DenseDoubleMatrix2D(n, n);
          DoubleMatrix2D wwi = new DenseDoubleMatrix2D(n, n);
  
          double dlx;
          double del;
  
          int nm1 = n - 1;
  
          DoubleMatrix2D vv = new DenseDoubleMatrix2D(nm1, nm1);
          DoubleMatrix2D xs = null;
          if (!approximateAlgorithm) {
              xs = new DenseDoubleMatrix2D(nm1, n);
          }
          DoubleMatrix1D s = new DenseDoubleMatrix1D(nm1);
          DoubleMatrix1D so = new DenseDoubleMatrix1D(nm1);
          DoubleMatrix1D x = new DenseDoubleMatrix1D(n - 1);
          DoubleMatrix1D ws;
          DoubleMatrix1D z = new DenseDoubleMatrix1D(nm1);
          int[] mm = new int[nm1];
          DoubleMatrix1D ro = new DenseDoubleMatrix1D(nm1);
  
          // shr warmStart sum(abs(offdiagonal(ss))).
          double shr = 0.0;
  
          for (int j = 0; j < n; j++) {
              for (int k = 0; k < n; k++) {
                  if (j == k) continue;
                  shr += Math.abs(ss.get(j, k));
              }
          }
  
  
          // TO TEST THE FOLLOWING CODE I NEED A DIAGONAL COVARIANCE MATRIX.
  
          // If ss is diagonal, just return the inverse of the covariance matrix diagonal
          // (penalized if necessary).
          if (shr == 0.0) {
              for (int j = 0; j < n; j++) {
                  if (Thread.currentThread().isInterrupted()) {
                      break;
                  }
  
                  if (!pen) {
                      ww.set(j, j, ss.get(j, j));
                  } else {
                      ww.set(j, j, ss.get(j, j) + rho.get(j, j));
                  }
                  wwi.set(j, j, 1.0 / Math.max(ww.get(j, j), eps));
              }
              return new Result(ww, wwi, niter, Double.NaN);
          }
  
  
          shr = getThr() * shr / nm1;
  
          if (approximateAlgorithm) {
              if (!warmStart) {
                  zero(wwi);
              }
  
              for (int m = 0; m < n; m++) {
                  System.out.println("m = " + m);
  
                  // This sets up vv, s, and r--i.e., W.11, s.12, and r.12.
                  setup(m, n, ss, rho, ss, vv, s, ro);
  
                  // This sets up x.12--i.e. theta.12.
                  int l = -1;
  
                  for (int j = 0; j < n; j++) {
                      if (Thread.currentThread().isInterrupted()) {
                          break;
                      }
  
                      if (j == m) continue;
                      l = l + 1;
                      x.set(l, wwi.get(j, m));
                  }
  
                  lasso(ro, nm1, vv, s, shr / n, x, z, mm);
  
                  l = -1;
                  for (int j = 0; j < n; j++) {
                      if (Thread.currentThread().isInterrupted()) {
                          break;
                      }
  
                      if (j == m) continue;
                      l = l + 1;
                      wwi.set(j, m, x.get(l));
                  }
              }
  
  
              niter = 1;
              return new Result(ww, wwi, niter, Double.NaN);
          }
  
  //        System.out.println("wwi = " + wwi);
  
          if (!warmStart) {
              ww.assign(ss);
  
              if (xs != null) {
                  zero(xs);
              }
          } else {
              for (int j = 0; j < n; j++) {
                  double xjj = -wwi.get(j, j);
  //                System.out.println("xjj = " + xjj);
                  int l = -1;
  
                  for (int k = 0; k < n; k++) {
                      if (Thread.currentThread().isInterrupted()) {
                          break;
                      }
  
                      if (k == j) continue;
                      l = l + 1;
                      xs.set(l, j, wwi.get(k, j) / xjj);
                  }
              }
          }
  
  