List of usage examples for org.apache.commons.math3.linear SingularValueDecomposition getRank
public int getRank()
From source file:com.opengamma.strata.math.impl.linearalgebra.SVDecompositionCommonsResult.java
/** * @param svd The result of the SV decomposition, not null *//*from ww w .j ava 2 s . com*/ public SVDecompositionCommonsResult(SingularValueDecomposition svd) { ArgChecker.notNull(svd, "svd"); _condition = svd.getConditionNumber(); _norm = svd.getNorm(); _rank = svd.getRank(); _s = CommonsMathWrapper.unwrap(svd.getS()); _singularValues = svd.getSingularValues(); _u = CommonsMathWrapper.unwrap(svd.getU()); _uTranspose = CommonsMathWrapper.unwrap(svd.getUT()); _v = CommonsMathWrapper.unwrap(svd.getV()); _vTranspose = CommonsMathWrapper.unwrap(svd.getVT()); _solver = svd.getSolver(); }
From source file:de.andreasschoknecht.LS3.LSSMCalculator.java
/** * Calculation of the singular value decomposition and initialization of corresponding variables. * * @param tdMatrix The weighted Term-Document Matrix for decomposition *///w ww . j av a 2 s.c o m void calculateSVD(double[][] tdMatrix) { weightedTDMatrix = new Array2DRowRealMatrix(tdMatrix); // execute actual SVD SingularValueDecomposition svd = new SingularValueDecomposition(weightedTDMatrix); U = svd.getU(); S = svd.getS(); Vt = svd.getVT(); // set rank of SVD rank = svd.getRank(); }
From source file:com.joptimizer.optimizers.OptimizationRequestHandler.java
/** * Find a solution of the linear (equalities) system A.x = b. * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 682" */// ww w. ja va 2 s. c om protected double[] findEqFeasiblePoint(double[][] A, double[] b) throws Exception { RealMatrix AT = new Array2DRowRealMatrix(A).transpose(); int p = A.length; SingularValueDecomposition dFact1 = new SingularValueDecomposition(AT); int rangoAT = dFact1.getRank(); if (rangoAT != p) { throw new RuntimeException("Equalities matrix A must have full rank"); } QRDecomposition dFact = new QRDecomposition(AT); //A = QR RealMatrix Q1Q2 = dFact.getQ(); RealMatrix R0 = dFact.getR(); RealMatrix Q1 = Q1Q2.getSubMatrix(0, AT.getRowDimension() - 1, 0, p - 1); RealMatrix R = R0.getSubMatrix(0, p - 1, 0, p - 1); double[][] rData = R.copy().getData(); //inversion double[][] rInvData = Utils.upperTriangularMatrixUnverse(rData); RealMatrix RInv = new Array2DRowRealMatrix(rInvData); //w = Q1 * Inv([R]T) . b double[] w = Q1.operate(RInv.transpose().operate(new ArrayRealVector(b))).toArray(); return w; }
From source file:com.joptimizer.optimizers.LPPresolverTest.java
/** * This problem has a deterministic solution. *//* w ww. ja v a 2 s . c om*/ public void testFromFile2() throws Exception { log.debug("testFromFile2"); String problemId = "2"; double[] c = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "c" + problemId + ".txt"); double[][] A = Utils.loadDoubleMatrixFromFile( "lp" + File.separator + "presolving" + File.separator + "A" + problemId + ".csv", ",".charAt(0)); double[] b = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "b" + problemId + ".txt"); double[] lb = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "lb" + problemId + ".txt"); double[] ub = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "ub" + problemId + ".txt"); double s = 0; try { s = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "s" + problemId + ".txt")[0]; } catch (Exception e) { } double[] expectedSolution = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "sol" + problemId + ".txt"); double expectedValue = Utils.loadDoubleArrayFromFile( "lp" + File.separator + "presolving" + File.separator + "value" + problemId + ".txt")[0]; double expectedTolerance = MatrixUtils.createRealMatrix(A) .operate(MatrixUtils.createRealVector(expectedSolution)).subtract(MatrixUtils.createRealVector(b)) .getNorm(); //must be: A pXn with rank(A)=p < n RealMatrix AMatrix = MatrixUtils.createRealMatrix(A); SingularValueDecomposition dec = new SingularValueDecomposition(AMatrix); int rankA = dec.getRank(); log.debug("p: " + AMatrix.getRowDimension()); log.debug("n: " + AMatrix.getColumnDimension()); log.debug("rank: " + rankA); LPPresolver lpPresolver = new LPPresolver(); lpPresolver.setNOfSlackVariables((short) s); lpPresolver.setExpectedSolution(expectedSolution);//this is just for test! lpPresolver.presolve(c, A, b, lb, ub); int n = lpPresolver.getPresolvedN(); //deterministic solution assertEquals(0, n); assertTrue(lpPresolver.getPresolvedC() == null); assertTrue(lpPresolver.getPresolvedA() == null); assertTrue(lpPresolver.getPresolvedB() == null); assertTrue(lpPresolver.getPresolvedLB() == null); assertTrue(lpPresolver.getPresolvedUB() == null); assertTrue(lpPresolver.getPresolvedYlb() == null); assertTrue(lpPresolver.getPresolvedYub() == null); assertTrue(lpPresolver.getPresolvedZlb() == null); assertTrue(lpPresolver.getPresolvedZub() == null); double[] sol = lpPresolver.postsolve(new double[] {}); assertEquals(expectedSolution.length, sol.length); for (int i = 0; i < sol.length; i++) { //log.debug("i: " + i); assertEquals(expectedSolution[i], sol[i], 1.e-9); } }
From source file:edu.cmu.tetrad.util.TetradMatrix.java
public int rank() { // return new RRQRDecomposition(apacheData).getRank(10); SingularValueDecomposition singularValueDecomposition = new SingularValueDecomposition(apacheData); return singularValueDecomposition.getRank(); }
From source file:com.joptimizer.optimizers.LPPresolverTest.java
private void doPresolving(double[] c, double[][] A, double[] b, double[] lb, double[] ub, double s, double[] expectedSolution, double expectedValue, double expectedTolerance) throws Exception { RealMatrix AMatrix = MatrixUtils.createRealMatrix(A); SingularValueDecomposition dec = new SingularValueDecomposition(AMatrix); int rankA = dec.getRank(); log.debug("p: " + AMatrix.getRowDimension()); log.debug("n: " + AMatrix.getColumnDimension()); log.debug("rank: " + rankA); LPPresolver lpPresolver = new LPPresolver(); lpPresolver.setNOfSlackVariables((short) s); lpPresolver.setExpectedSolution(expectedSolution);//this is just for test! //lpPresolver.setExpectedTolerance(expectedTolerance);//this is just for test! lpPresolver.presolve(c, A, b, lb, ub); int n = lpPresolver.getPresolvedN(); double[] presolvedC = lpPresolver.getPresolvedC().toArray(); double[][] presolvedA = lpPresolver.getPresolvedA().toArray(); double[] presolvedB = lpPresolver.getPresolvedB().toArray(); double[] presolvedLb = lpPresolver.getPresolvedLB().toArray(); double[] presolvedUb = lpPresolver.getPresolvedUB().toArray(); double[] presolvedYlb = lpPresolver.getPresolvedYlb().toArray(); double[] presolvedYub = lpPresolver.getPresolvedYub().toArray(); double[] presolvedZlb = lpPresolver.getPresolvedZlb().toArray(); double[] presolvedZub = lpPresolver.getPresolvedZub().toArray(); log.debug("n : " + n); log.debug("presolvedC : " + ArrayUtils.toString(presolvedC)); log.debug("presolvedA : " + ArrayUtils.toString(presolvedA)); log.debug("presolvedB : " + ArrayUtils.toString(presolvedB)); log.debug("presolvedLb : " + ArrayUtils.toString(presolvedLb)); log.debug("presolvedUb : " + ArrayUtils.toString(presolvedUb)); log.debug("presolvedYlb: " + ArrayUtils.toString(presolvedYlb)); log.debug("presolvedYub: " + ArrayUtils.toString(presolvedYub)); log.debug("presolvedZlb: " + ArrayUtils.toString(presolvedZlb)); log.debug("presolvedZub: " + ArrayUtils.toString(presolvedZub)); //check objective function double delta = expectedTolerance; RealVector presolvedX = MatrixUtils.createRealVector(lpPresolver.presolve(expectedSolution)); log.debug("presolved value: " + MatrixUtils.createRealVector(presolvedC).dotProduct(presolvedX)); RealVector postsolvedX = MatrixUtils.createRealVector(lpPresolver.postsolve(presolvedX.toArray())); double value = MatrixUtils.createRealVector(c).dotProduct(postsolvedX); assertEquals(expectedValue, value, delta); //check postsolved constraints for (int i = 0; i < lb.length; i++) { double di = Double.isNaN(lb[i]) ? -Double.MAX_VALUE : lb[i]; assertTrue(di <= postsolvedX.getEntry(i) + delta); }/*w w w .ja v a 2 s. c om*/ for (int i = 0; i < ub.length; i++) { double di = Double.isNaN(ub[i]) ? Double.MAX_VALUE : ub[i]; assertTrue(di + delta >= postsolvedX.getEntry(i)); } RealVector Axmb = AMatrix.operate(postsolvedX).subtract(MatrixUtils.createRealVector(b)); assertEquals(0., Axmb.getNorm(), expectedTolerance); //check presolved constraints assertEquals(presolvedLb.length, presolvedX.getDimension()); assertEquals(presolvedUb.length, presolvedX.getDimension()); AMatrix = MatrixUtils.createRealMatrix(presolvedA); RealVector bvector = MatrixUtils.createRealVector(presolvedB); for (int i = 0; i < presolvedLb.length; i++) { double di = Double.isNaN(presolvedLb[i]) ? -Double.MAX_VALUE : presolvedLb[i]; assertTrue(di <= presolvedX.getEntry(i) + delta); } for (int i = 0; i < presolvedUb.length; i++) { double di = Double.isNaN(presolvedUb[i]) ? Double.MAX_VALUE : presolvedUb[i]; assertTrue(di + delta >= presolvedX.getEntry(i)); } Axmb = AMatrix.operate(presolvedX).subtract(bvector); assertEquals(0., Axmb.getNorm(), expectedTolerance); //check rank(A): must be A pXn with rank(A)=p < n AMatrix = MatrixUtils.createRealMatrix(presolvedA); dec = new SingularValueDecomposition(AMatrix); rankA = dec.getRank(); log.debug("p: " + AMatrix.getRowDimension()); log.debug("n: " + AMatrix.getColumnDimension()); log.debug("rank: " + rankA); assertEquals(AMatrix.getRowDimension(), rankA); assertTrue(rankA < AMatrix.getColumnDimension()); }
From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ.java
public boolean isIndependent(Node x, Node y, List<Node> z) { if (verbose) { System.out.println("\n" + x + " _||_ " + y + " | " + z); }//from www .j a va 2 s . c o m List<Node> aa = nodeMap.get(x); List<Node> bb = nodeMap.get(y); List<Node> cc = new ArrayList<Node>(); for (Node _z : z) { cc.addAll(nodeMap.get(_z)); } TetradMatrix submatrix = subMatrix(cov, aa, bb, cc); TetradMatrix inverse; int rank; try { inverse = submatrix.inverse(); rank = inverse.columns(); } catch (Exception e) { SingularValueDecomposition svd = new SingularValueDecomposition(submatrix.getRealMatrix()); RealMatrix _inverse = svd.getSolver().getInverse(); inverse = new TetradMatrix(_inverse, _inverse.getRowDimension(), _inverse.getColumnDimension()); rank = svd.getRank(); } final List<Double> pValues = new ArrayList<Double>(); List<Integer> _i = new ArrayList<Integer>(); List<Integer> _m = new ArrayList<Integer>(); TetradMatrix stdSubmatrix = TestHippocampusUtils.subset(stdData, aa, bb, cc); for (int i = 0; i < aa.size(); i++) { for (int m = 0; m < bb.size(); m++) { int j = aa.size() + m; double a = -1.0 * inverse.get(i, j); double v0 = inverse.get(i, i); double v1 = inverse.get(j, j); double b = sqrt(v0 * v1); double r = a / b; int dof = cov.getSampleSize() - 1 - rank; if (dof < 0) { // out.println("Negative dof: " + dof + " n = " + cov.getSampleSize() + " cols = " + inverse.columns()); dof = 0; } double t2 = moment22(stdSubmatrix, i, j); double z_ = sqrt(dof) * 0.5 * (Math.log(1.0 + r) - Math.log(1.0 - r)); double p = 2.0 * (1.0 - RandomUtil.getInstance().normalCdf(0, sqrt(t2), abs(z_))); // if (p < alpha) p = 0; // System.out.println("z = " + z_ + " p = " + p); pValues.add(p); _i.add(i); _m.add(m); } } List<Integer> indices = new ArrayList<Integer>(); for (int i = 0; i < pValues.size(); i++) { indices.add(i); } Collections.sort(indices, new Comparator<Integer>() { @Override public int compare(Integer o1, Integer o2) { return pValues.get(o1).compareTo(pValues.get(o2)); } }); List<Double> pValues2 = new ArrayList<Double>(); List<Integer> _iIndices = new ArrayList<Integer>(); List<Integer> _mIndices = new ArrayList<Integer>(); for (int _y = 0; _y < indices.size(); _y++) { pValues2.add(pValues.get(indices.get(_y))); _iIndices.add(_i.get(indices.get(_y))); _mIndices.add(_m.get(indices.get(_y))); } int k = StatUtils.fdr(alpha, pValues2, false); double cutoff = StatUtils.fdrCutoff(alpha, pValues2, false); System.out.println("** cutoff = " + cutoff); int nonzero = -1; for (int i = 0; i < pValues2.size(); i++) { if (pValues2.get(i) != 0) { nonzero = i; break; } } double q = StatUtils.fdrQ(pValues2, nonzero); System.out.println("Q = " + q); // boolean independent = q >= alpha; //// boolean independent = k > nonzero; // System.out.println("k = " + k); System.out.println("nonzero = " + nonzero); System.out.println("fisher = " + fisherMethodP(pValues2)); // // boolean independent = k <= nonzero; double ratio = ((double) k) / pValues2.size(); System.out.println("ratio = " + ratio); boolean independent = ratio < 0.1; // System.out.println(independent ? "Independent" : "Dependent"); return independent; }
From source file:edu.cmu.tetrad.search.IndTestMultiCci.java
public boolean isIndependent(Node x, Node y, List<Node> z) { if (verbose) { System.out.println("\n" + x + " _||_ " + y + " | " + z); out.println("\n" + x + " _||_ " + y + " | " + z); }/* ww w . ja v a 2 s.c o m*/ List<Node> aa = nodeMap.get(x); List<Node> bb = nodeMap.get(y); // int[][] twod = new int[aa.size()][bb.size()]; List<Node> cc = new ArrayList<Node>(); for (Node _z : z) { cc.addAll(nodeMap.get(_z)); } TetradMatrix submatrix = TestHippocampusUtils.subMatrix(cov, aa, bb, cc); TetradMatrix inverse; int rank; try { inverse = submatrix.inverse(); rank = inverse.columns(); } catch (Exception e) { SingularValueDecomposition svd = new SingularValueDecomposition(submatrix.getRealMatrix()); RealMatrix _inverse = svd.getSolver().getInverse(); inverse = new TetradMatrix(_inverse, _inverse.getRowDimension(), _inverse.getColumnDimension()); rank = svd.getRank(); } final List<Double> pValues = new ArrayList<Double>(); List<Integer> _i = new ArrayList<Integer>(); List<Integer> _m = new ArrayList<Integer>(); System.out.println("# voxels for " + x.getName() + " = " + aa.size()); System.out.println("# voxels for " + y.getName() + " = " + bb.size()); System.out.println("# p values = " + aa.size() * bb.size()); IndTestConditionalCorrelation cciTest = new IndTestConditionalCorrelation(dataSet, alpha); for (int i = 0; i < aa.size(); i++) { for (int m = 0; m < bb.size(); m++) { int j = aa.size() + m; double a = -1.0 * inverse.get(i, j); double v0 = inverse.get(i, i); double v1 = inverse.get(j, j); double b = Math.sqrt(v0 * v1); double r = a / b; int dof = cov.getSampleSize() - 1 - rank; if (dof < 0) { out.println("Negative dof: " + dof + " n = " + cov.getSampleSize() + " cols = " + inverse.columns()); dof = 0; } double z_ = Math.sqrt(dof) * 0.5 * (Math.log(1.0 + r) - Math.log(1.0 - r)); double p = 2.0 * (1.0 - RandomUtil.getInstance().normalCdf(0, 1, abs(z_))); cciTest.isIndependent(aa.get(i), bb.get(m), cc); pValues.add(p); if (m == 0) { System.out.println("i = " + i + " m = " + m + " p = " + p); } _i.add(i); _m.add(m); } } List<Integer> indices = new ArrayList<Integer>(); for (int i = 0; i < pValues.size(); i++) { indices.add(i); } Collections.sort(indices, new Comparator<Integer>() { @Override public int compare(Integer o1, Integer o2) { return pValues.get(o1).compareTo(pValues.get(o2)); } }); List<Double> pValues2 = new ArrayList<Double>(); List<Integer> _iIndices = new ArrayList<Integer>(); List<Integer> _mIndices = new ArrayList<Integer>(); for (int _y = 0; _y < indices.size(); _y++) { pValues2.add(pValues.get(indices.get(_y))); _iIndices.add(_i.get(indices.get(_y))); _mIndices.add(_m.get(indices.get(_y))); } int k = StatUtils.fdr(alpha, pValues2, false); double cutoff = StatUtils.fdrCutoff(alpha, pValues2, false); int nonzero = -1; for (int i = 0; i < pValues2.size(); i++) { if (pValues2.get(i) != 0) { nonzero = i; break; } } boolean dependent = k > nonzero; boolean independent = !dependent; TetradMatrix xCoords = coords.get(x); TetradMatrix yCoords = coords.get(y); int[][][] X = TestHippocampusUtils.threeDView(_iIndices, k, xCoords, independent, nonzero); int[][][] Y = TestHippocampusUtils.threeDView(_mIndices, k, yCoords, independent, nonzero); all3D.add(X); all3D.add(Y); String fact = SearchLogUtils.independenceFact(x, y, z); // Printing out stuff for Ruben. First print out dependent voxels. // 2.Second file, contain a list of all the dependencies between voxels. // // 10 -- 50 // 30 -- 2 int[][][] Cx = getC(X); int[][][] Cy = getC(Y); out.println("\n\n" + fact); for (int g = independent ? nonzero : 0; g < k; g++) { int i = getIndex(_iIndices, Cx, g, coords.get(x)); int j = getIndex(_mIndices, Cy, g, coords.get(y)); if (i == -1 || j == -1) throw new IllegalArgumentException(); out.println(i + " -- " + j); } out.println(); // 1. First file, containing info of both ROIs and all their voxels. // Example: // // ROI_LABEL voxel_LABEL COORDINATES #Dependencies // ENT 10 -80 50 38 6 // CA1 50 -70 15 90 2 printDependencies(x, fact, X, Cx); printDependencies(y, fact, Y, Cy); // OK back to work. int xCount = countAboveThreshold(X, threshold); int yCount = countAboveThreshold(Y, threshold); System.out.println("Total above threshold count = " + (xCount + yCount)); out.println("Total above threshold count = " + (xCount + yCount)); boolean thresholdIndep = !(xCount > 0 && yCount > 0); String projection; projection = "Axial"; TestHippocampusUtils.printChart(X, xCoords, 0, 1, x.getName(), projection, fact, false, false, out, threshold); TestHippocampusUtils.printChart(Y, yCoords, 0, 1, y.getName(), projection, fact, false, false, out, threshold); projection = "Coronal"; TestHippocampusUtils.printChart(X, xCoords, 0, 2, x.getName(), projection, fact, true, false, out, threshold); TestHippocampusUtils.printChart(Y, yCoords, 0, 2, y.getName(), projection, fact, true, false, out, threshold); projection = "Saggital"; TestHippocampusUtils.printChart(X, xCoords, 1, 2, x.getName(), projection, fact, true, false, out, threshold); TestHippocampusUtils.printChart(Y, yCoords, 1, 2, y.getName(), projection, fact, true, false, out, threshold); if (thresholdIndep) { // if (independent) { if (verbose) { System.out.println("Independent"); out.println("Independent"); } out.flush(); return true; } else { if (verbose) { System.out.println("Dependent\n"); out.println("Dependent\n"); } out.flush(); return false; } }
From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ4.java
public boolean isIndependent(Node x, Node y, List<Node> z) { if (verbose) { System.out.println("\n" + x + " _||_ " + y + " | " + z); out.println("\n" + x + " _||_ " + y + " | " + z); }// w w w .j av a 2 s.c o m List<Node> aa = nodeMap.get(x); List<Node> bb = nodeMap.get(y); // int[][] twod = new int[aa.size()][bb.size()]; List<Node> cc = new ArrayList<Node>(); for (Node _z : z) { cc.addAll(nodeMap.get(_z)); } TetradMatrix submatrix = TestHippocampusUtils.subMatrix(cov, aa, bb, cc); TetradMatrix stdSubmatrix = TestHippocampusUtils.subset(stdData, aa, bb, cc); TetradMatrix inverse; int rank; try { inverse = submatrix.inverse(); rank = inverse.columns(); } catch (Exception e) { SingularValueDecomposition svd = new SingularValueDecomposition(submatrix.getRealMatrix()); RealMatrix _inverse = svd.getSolver().getInverse(); inverse = new TetradMatrix(_inverse, _inverse.getRowDimension(), _inverse.getColumnDimension()); rank = svd.getRank(); } final List<Double> pValues = new ArrayList<Double>(); List<Integer> _i = new ArrayList<Integer>(); List<Integer> _m = new ArrayList<Integer>(); System.out.println("# voxels for " + x.getName() + " = " + aa.size()); System.out.println("# voxels for " + y.getName() + " = " + bb.size()); System.out.println("# p values = " + aa.size() * bb.size()); for (int i = 0; i < aa.size(); i++) { for (int m = 0; m < bb.size(); m++) { int j = aa.size() + m; double a = -1.0 * inverse.get(i, j); double v0 = inverse.get(i, i); double v1 = inverse.get(j, j); double b = sqrt(v0 * v1); double r = a / b; int dof = cov.getSampleSize() - 1 - rank; if (dof < 0) { out.println("Negative dof: " + dof + " n = " + cov.getSampleSize() + " cols = " + inverse.columns()); dof = 0; } double t2 = moment22(stdSubmatrix, i, j); double z_ = sqrt(dof) * 0.5 * (Math.log(1.0 + r) - Math.log(1.0 - r)); double p = 2.0 * (1.0 - RandomUtil.getInstance().normalCdf(0, sqrt(t2), abs(z_))); // System.out.println("z = " + z_ + " p = " + p); pValues.add(p); _i.add(i); _m.add(m); } } List<Integer> indices = new ArrayList<Integer>(); for (int i = 0; i < pValues.size(); i++) { indices.add(i); } Collections.sort(indices, new Comparator<Integer>() { @Override public int compare(Integer o1, Integer o2) { return pValues.get(o1).compareTo(pValues.get(o2)); } }); List<Double> pValues2 = new ArrayList<Double>(); List<Integer> _iIndices = new ArrayList<Integer>(); List<Integer> _mIndices = new ArrayList<Integer>(); for (int _y = 0; _y < indices.size(); _y++) { pValues2.add(pValues.get(indices.get(_y))); _iIndices.add(_i.get(indices.get(_y))); _mIndices.add(_m.get(indices.get(_y))); } int k = StatUtils.fdr(alpha, pValues2, false); // double cutoff = StatUtils.fdrCutoff(alpha, pValues2, false); double cutoff = pValues2.get(k); int kP = -1; for (int j = 0; j < pValues2.size(); j++) { if (pValues2.get(j) > 0) { kP = j; break; } } double q = StatUtils.fdrQ(pValues2, kP); System.out.println("Q = " + q); // for (int k3 = 0; k3 < 500; k3++) { // System.out.println("Q for " + k3 + " = " + StatUtils.fdrQ(pValues2, k3)); // } System.out.println("** cutoff = " + cutoff); int nonzero = -1; for (int i = 0; i < pValues2.size(); i++) { if (pValues2.get(i) != 0) { nonzero = i; break; } } boolean dependent = k > nonzero; boolean independent = !dependent; TetradMatrix xCoords = coords.get(x); TetradMatrix yCoords = coords.get(y); int[][][] X = TestHippocampusUtils.threeDView(_iIndices, k, xCoords, independent, nonzero); int[][][] Y = TestHippocampusUtils.threeDView(_mIndices, k, yCoords, independent, nonzero); all3D.add(X); all3D.add(Y); String fact = SearchLogUtils.independenceFact(x, y, z); // Printing out stuff for Ruben. First print out dependent voxels. // 2.Second file, contain a list of all the dependencies between voxels. // // 10 -- 50 // 30 -- 2 int[][][] Cx = getC(X); int[][][] Cy = getC(Y); out.println("\n\n" + fact); for (int g = independent ? nonzero : 0; g < k; g++) { int i = getIndex(_iIndices, Cx, g, coords.get(x)); int j = getIndex(_mIndices, Cy, g, coords.get(y)); if (i == -1 || j == -1) throw new IllegalArgumentException(); out.println(i + " -- " + j); } out.println(); // 1. First file, containing info of both ROIs and all their voxels. // Example: // // ROI_LABEL voxel_LABEL COORDINATES #Dependencies // ENT 10 -80 50 38 6 // CA1 50 -70 15 90 2 printDependencies(x, fact, X, Cx); printDependencies(y, fact, Y, Cy); // OK back to work. int xCount = countAboveThreshold(X, threshold); int yCount = countAboveThreshold(Y, threshold); System.out.println("Total above threshold count = " + (xCount + yCount)); out.println("Total above threshold count = " + (xCount + yCount)); boolean thresholdIndep = !(xCount > 0 && yCount > 0); String projection; projection = "Axial"; TestHippocampusUtils.printChart(X, xCoords, 0, 1, x.getName(), projection, fact, false, false, out, threshold); TestHippocampusUtils.printChart(Y, yCoords, 0, 1, y.getName(), projection, fact, false, false, out, threshold); projection = "Coronal"; TestHippocampusUtils.printChart(X, xCoords, 0, 2, x.getName(), projection, fact, true, false, out, threshold); TestHippocampusUtils.printChart(Y, yCoords, 0, 2, y.getName(), projection, fact, true, false, out, threshold); projection = "Saggital"; TestHippocampusUtils.printChart(X, xCoords, 1, 2, x.getName(), projection, fact, true, false, out, threshold); TestHippocampusUtils.printChart(Y, yCoords, 1, 2, y.getName(), projection, fact, true, false, out, threshold); if (thresholdIndep) { // if (independent) { if (verbose) { System.out.println("Independent"); out.println("Independent"); } out.flush(); return true; } else { if (verbose) { System.out.println("Dependent\n"); out.println("Dependent\n"); } out.flush(); return false; } }
From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ2.java
private boolean indepCollection2(List<Node> aa, List<Node> bb, List<Node> cc, int n, double alpha) { // List<Double> ret = getCutoff1(aa, bb, cc); // List<Double> ret = getCutoff2(aa, bb, cc); // double mean = ret.get(0); // double sd = ret.get(1); int numPerm = 10; TetradMatrix submatrix = subMatrix(cov, aa, bb, cc); TetradMatrix inverse;/*from w ww. ja v a2s . co m*/ int rank; try { inverse = submatrix.inverse(); rank = inverse.columns(); } catch (Exception e) { SingularValueDecomposition svd = new SingularValueDecomposition(submatrix.getRealMatrix()); RealMatrix _inverse = svd.getSolver().getInverse(); inverse = new TetradMatrix(_inverse, _inverse.getRowDimension(), _inverse.getColumnDimension()); rank = svd.getRank(); } List<Double> pValues = new ArrayList<Double>(); for (int i = 0; i < aa.size(); i++) { for (int m = 0; m < bb.size(); m++) { int j = aa.size() + m; double a = -1.0 * inverse.get(i, j); double v0 = inverse.get(i, i); double v1 = inverse.get(j, j); double b = Math.sqrt(v0 * v1); double r = a / b; int dof = n - 1 - rank; if (dof < 0) { System.out.println("Negative dof: " + dof + " n = " + n + " cols = " + inverse.columns()); dof = 0; } double z = Math.sqrt(dof) * 0.5 * (Math.log(1.0 + r) - Math.log(1.0 - r)); double p = 2.0 * (1.0 - RandomUtil.getInstance().normalCdf(0, 1, abs(z))); pValues.add(p); } } List<Double> zeroes = new ArrayList<Double>(); for (double d : pValues) { if (d == 0) zeroes.add(d); } int k = 0; for (double p : pValues) { if (p < alpha) { k++; } } numTests++; int count = countGreater(aa, bb, cc, numPerm, k); // int count = countGreater2(aa, bb, cc, numPerm, k); double p = count / (double) numPerm; boolean indep = p > alpha; // double p = (1.0 - RandomUtil.getInstance().normalCdf(0, 1, (k - mean) / sd)); // boolean indep = p > alpha; if (verbose) { // System.out.println("\n#accepted " + k + " cutoff = " + kCutoff + " indep = " + indep); // System.out.println("\n#accepted " + k + " meam = " + mean + " sd = " + sd + " indep = " + indep); // System.out.println("standard deviations " + (k - mean) / sd + " p = " + p + " alpha = " + alpha); System.out.println("\n#accepted " + k + " indep = " + indep); System.out.println("p = " + p + " alpha = " + alpha); } numTests++; if (!indep) numDependent++; return indep; }