List of usage examples for org.apache.commons.math3.linear RealMatrix getRowDimension
int getRowDimension();
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); }/* w ww . 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:msi.gama.util.matrix.GamaFloatMatrix.java
public IMatrix fromApacheMatrix(final IScope scope, final RealMatrix rm) { if (rm == null) { return null; }/*ww w . jav a 2 s. c o m*/ final GamaFloatMatrix matrix = new GamaFloatMatrix(rm.getColumnDimension(), rm.getRowDimension()); for (int i = 0; i < numCols; i++) { for (int j = 0; j < numRows; j++) { matrix.set(scope, i, j, rm.getEntry(j, i)); } } return matrix; }
From source file:edu.cudenver.bios.matrix.DesignEssenceMatrix.java
/** * Fills in a random column in the full design matrix * * @param randomColumn column index in random submatrix * @param fullColumn column index in full design matrix * @param fullDesign full design matrix//from w w w. j a va 2 s.c o m */ private void fillRandomColumn(int randomColumn, int fullColumn, RealMatrix fullDesign) { // if the column represents a random predictor, build a normal distribution // from which to pull random values NormalDistribution dist = null; // note, the jsc library takes a standard deviation, not a variance so // we take the square root dist = new NormalDistribution(randomColMetaData[randomColumn].getMean(), Math.sqrt(randomColMetaData[randomColumn].getVariance())); dist.reseedRandomGenerator(randomSeed); for (int row = 0; row < fullDesign.getRowDimension(); row++) { // fill in the data fullDesign.setEntry(row, fullColumn, dist.sample()); } }
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); }/*from w ww .java 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:com.bolatu.gezkoncsvlogger.GyroOrientation.RotationKalmanFilter.java
/** * Correct the current state estimate with an actual measurement. * * @param z/*from ww w.j a v a2 s.c o m*/ * the measurement vector * @throws NullArgumentException * if the measurement vector is {@code null} * @throws DimensionMismatchException * if the dimension of the measurement vector does not fit * @throws SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) throws NullArgumentException, DimensionMismatchException, SingularMatrixException { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance).multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 // instead of calculating the inverse of S we can rearrange the formula, // and then solve the linear equation A x X = B with A = S', X = K' and // B = (H * P)' // K(k) * S = P(k)- * H' // S' * K(k)' = H * P(k)-' RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver() .solve(measurementMatrix.multiply(errorCovariance.transpose())).transpose(); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); }
From source file:edu.cmu.tetrad.data.DataUtils.java
private static RealMatrix times(final RealMatrix m, final RealMatrix n) { if (m.getColumnDimension() != n.getRowDimension()) throw new IllegalArgumentException("Incompatible matrices."); final int rowDimension = m.getRowDimension(); final int columnDimension = n.getColumnDimension(); final RealMatrix out = new BlockRealMatrix(rowDimension, columnDimension); final int NTHREADS = Runtime.getRuntime().availableProcessors(); final int all = rowDimension; ForkJoinPool pool = ForkJoinPoolInstance.getInstance().getPool(); for (int t = 0; t < NTHREADS; t++) { final int _t = t; Runnable worker = new Runnable() { @Override//from w ww .java 2 s . c o m public void run() { int chunk = all / NTHREADS + 1; for (int row = _t * chunk; row < Math.min((_t + 1) * chunk, all); row++) { if ((row + 1) % 100 == 0) System.out.println(row + 1); for (int col = 0; col < columnDimension; ++col) { double sum = 0.0D; int commonDimension = m.getColumnDimension(); for (int i = 0; i < commonDimension; ++i) { sum += m.getEntry(row, i) * n.getEntry(i, col); } // double sum = m.getRowVector(row).dotProduct(n.getColumnVector(col)); out.setEntry(row, col, sum); } } } }; pool.submit(worker); } while (!pool.isQuiescent()) { } // for (int row = 0; row < rowDimension; ++row) { // if ((row + 1) % 100 == 0) System.out.println(row + 1); // // for (int col = 0; col < columnDimension; ++col) { // double sum = 0.0D; // // int commonDimension = m.getColumnDimension(); // // for (int i = 0; i < commonDimension; ++i) { // sum += m.getEntry(row, i) * n.getEntry(i, col); // } // // out.setEntry(row, col, sum); // } // } return out; }
From source file:edu.cmu.tetrad.data.DataUtils.java
public static TetradMatrix cov2(TetradMatrix data) { RealMatrix covarianceMatrix = new Covariance(data.getRealMatrix()).getCovarianceMatrix(); return new TetradMatrix(covarianceMatrix, covarianceMatrix.getRowDimension(), covarianceMatrix.getColumnDimension()); }
From source file:com.clust4j.algo.DBSCAN.java
/** {@inheritDoc} */ @Override/* ww w. j av a 2 s .c om*/ public int[] predict(RealMatrix newData) { final int[] fit_labels = getLabels(); // propagates errors final int n = newData.getColumnDimension(); // Make sure matches dimensionally if (n != this.data.getColumnDimension()) throw new DimensionMismatchException(n, data.getColumnDimension()); // Fit a radius model RadiusNeighbors radiusModel = new RadiusNeighborsParameters(eps) // no scale necessary; may already have been done .setMetric(dist_metric).setSeed(getSeed()).fitNewModel(data); final int[] newLabels = new int[newData.getRowDimension()]; Neighborhood theHood = radiusModel.getNeighbors(newData); int[][] indices = theHood.getIndices(); int[] idx_row; for (int i = 0; i < indices.length; i++) { idx_row = indices[i]; int current_class = NOISE_CLASS; if (idx_row.length == 0) { /* * If there are no indices in this point's radius, * we can just avoid the next step and exit early */ } else { // otherwise, we know there is something in the radius--noise or other int j = 0; while (j < idx_row.length) { current_class = fit_labels[idx_row[j]]; /* * The indices are ordered ascendingly by dist. * Even if the closest point is a noise point, it * could be within a border point's radius, so we * need to keep going. */ if (NOISE_CLASS == current_class) { j++; } else { break; } } } newLabels[i] = current_class; } return newLabels; }
From source file:edu.cmu.tetrad.data.DataUtils.java
public static void simpleTest() { double[][] d = new double[][] { { 1, 2 }, { 3, 4 }, { 5, 6 }, }; RealMatrix m = new BlockRealMatrix(d); System.out.println(m);/* w w w . j a v a2 s . c om*/ System.out.println(times(m.transpose(), m)); System.out.println(MatrixUtils.transposeWithoutCopy(m).multiply(m)); TetradMatrix n = new TetradMatrix(m, m.getRowDimension(), m.getColumnDimension()); System.out.println(n); RealMatrix q = n.getRealMatrix(); RealMatrix q1 = MatrixUtils.transposeWithoutCopy(q); RealMatrix q2 = times(q1, q); System.out.println(new TetradMatrix(q2, q.getColumnDimension(), q.getColumnDimension())); }
From source file:com.clust4j.algo.HierarchicalAgglomerative.java
protected HierarchicalAgglomerative(RealMatrix data, HierarchicalAgglomerativeParameters planner) { super(data, planner, planner.getNumClusters()); this.linkage = planner.getLinkage(); if (!isValidMetric(this.dist_metric)) { warn(this.dist_metric.getName() + " is invalid for " + this.linkage + ". Falling back to default Euclidean dist"); setSeparabilityMetric(DEF_DIST); }//from ww w . j a va2 s . co m this.m = data.getRowDimension(); this.num_clusters = super.k; logModelSummary(); }