use of cern.colt.matrix.linalg.Algebra in project Gemma by PavlidisLab.
the class ExpressionDataSVD method removeHighestComponents.
/**
* Provide a reconstructed matrix removing the first N components (the most significant ones). If the matrix was
* normalized first, removing the first component replicates the normalization approach taken by Nielsen et al.
* (Lancet 359, 2002) and Alter et al. (PNAS 2000). Correction by ANOVA would yield similar results if the nuisance
* variable is known.
*
* @param numComponentsToRemove The number of components to remove, starting from the largest eigenvalue.
* @return the reconstructed matrix; values that were missing before are re-masked.
*/
public ExpressionDataDoubleMatrix removeHighestComponents(int numComponentsToRemove) {
DoubleMatrix<Integer, Integer> copy = svd.getS().copy();
for (int i = 0; i < numComponentsToRemove; i++) {
copy.set(i, i, 0.0);
}
double[][] rawU = svd.getU().getRawMatrix();
double[][] rawS = copy.getRawMatrix();
double[][] rawV = svd.getV().getRawMatrix();
DoubleMatrix2D u = new DenseDoubleMatrix2D(rawU);
DoubleMatrix2D s = new DenseDoubleMatrix2D(rawS);
DoubleMatrix2D v = new DenseDoubleMatrix2D(rawV);
Algebra a = new Algebra();
DoubleMatrix<CompositeSequence, BioMaterial> reconstructed = new DenseDoubleMatrix<>(a.mult(a.mult(u, s), a.transpose(v)).toArray());
reconstructed.setRowNames(this.expressionData.getMatrix().getRowNames());
reconstructed.setColumnNames(this.expressionData.getMatrix().getColNames());
// re-mask the missing values.
for (int i = 0; i < reconstructed.rows(); i++) {
for (int j = 0; j < reconstructed.columns(); j++) {
if (Double.isNaN(this.missingValueInfo.get(i, j))) {
reconstructed.set(i, j, Double.NaN);
}
}
}
return new ExpressionDataDoubleMatrix(this.expressionData, reconstructed);
}
use of cern.colt.matrix.linalg.Algebra in project tetrad by cmu-phil.
the class StabilityUtils method totalInstabilityUndir.
// needs a symmetric matrix
// array of averages of instability matrix over [all, cc, cd, dd] edges
public static double[] totalInstabilityUndir(DoubleMatrix2D xi, List<Node> vars) {
if (vars.size() != xi.columns() || vars.size() != xi.rows()) {
throw new IllegalArgumentException("stability mat must have same number of rows and columns as there are vars");
}
Algebra al = new Algebra();
// DoubleMatrix2D xiu = MGM.upperTri(xi.copy().assign(al.transpose(xi)),1);
DoubleMatrix2D xiu = xi.copy().assign(xi.copy().assign(Functions.minus(1.0)), Functions.mult).assign(Functions.mult(-2.0));
double[] D = new double[4];
int[] discInds = MixedUtils.getDiscreteInds(vars);
int[] contInds = MixedUtils.getContinuousInds(vars);
int p = contInds.length;
int q = discInds.length;
double temp = MGM.upperTri(xiu.copy(), 1).zSum();
D[0] = temp / ((p + q - 1.0) * (p + q) / 2.0);
temp = MGM.upperTri(xiu.viewSelection(contInds, contInds).copy(), 1).zSum();
D[1] = temp / (p * (p - 1.0) / 2.0);
temp = xiu.viewSelection(contInds, discInds).zSum();
D[2] = temp / (p * q);
temp = MGM.upperTri(xiu.viewSelection(discInds, discInds).copy(), 1).zSum();
D[3] = temp / (q * (q - 1.0) / 2.0);
return D;
}
use of cern.colt.matrix.linalg.Algebra in project tetrad by cmu-phil.
the class IndTestRegression method determines.
// ==========================PRIVATE METHODS============================//
// /**
// * Return the p-value of the last calculated independence fact.
// *
// * @return this p-value. When accessed through the IndependenceChecker
// * interface, this p-value should only be considered to be a
// * relative strength.
// */
// private double getRelativeStrength() {
//
// // precondition: pdf is the most recently used partial
// // correlation distribution function, and storedR is the most
// // recently calculated partial correlation.
// return 2.0 * Integrator.getArea(npdf, Math.abs(storedR), 9.0, 100);
// }
// /**
// * Computes that value x such that P(abs(N(0,1) > x) < alpha. Note that
// * this is a two sided test of the null hypothesis that the Fisher's Z
// * value, which is distributed as N(0,1) is not equal to 0.0.
// */
// private double cutoffGaussian() {
// npdf = new NormalPdf();
// final double upperBound = 9.0;
// final double delta = 0.001;
// // double alpha = this.alpha/2.0; //Two sided test
// return CutoffFinder.getCutoff(npdf, upperBound, alpha, delta);
// }
// private int sampleSize() {
// return data.rows();
// }
public boolean determines(List<Node> zList, Node xVar) {
if (zList == null) {
throw new NullPointerException();
}
for (Node node : zList) {
if (node == null) {
throw new NullPointerException();
}
}
int size = zList.size();
int[] zCols = new int[size];
int xIndex = getVariables().indexOf(xVar);
for (int i = 0; i < zList.size(); i++) {
zCols[i] = getVariables().indexOf(zList.get(i));
}
int[] zRows = new int[data.rows()];
for (int i = 0; i < data.rows(); i++) {
zRows[i] = i;
}
DoubleMatrix2D Z = data.viewSelection(zRows, zCols);
DoubleMatrix1D x = data.viewColumn(xIndex);
DoubleMatrix2D Zt = new Algebra().transpose(Z);
DoubleMatrix2D ZtZ = new Algebra().mult(Zt, Z);
DoubleMatrix2D G = new DenseDoubleMatrix2D(new TetradMatrix(ZtZ.toArray()).inverse().toArray());
// Bug in Colt? Need to make a copy before multiplying to avoid
// a ClassCastException.
DoubleMatrix2D Zt2 = Zt.like();
Zt2.assign(Zt);
DoubleMatrix2D GZt = new Algebra().mult(G, Zt2);
DoubleMatrix1D b_x = new Algebra().mult(GZt, x);
DoubleMatrix1D xPred = new Algebra().mult(Z, b_x);
DoubleMatrix1D xRes = xPred.copy().assign(x, Functions.minus);
double SSE = xRes.aggregate(Functions.plus, Functions.square);
boolean determined = SSE < 0.0001;
if (determined) {
StringBuilder sb = new StringBuilder();
sb.append("Determination found: ").append(xVar).append(" is determined by {");
for (int i = 0; i < zList.size(); i++) {
sb.append(zList.get(i));
if (i < zList.size() - 1) {
sb.append(", ");
}
}
sb.append("}");
TetradLogger.getInstance().log("independencies", sb.toString());
}
return determined;
}
use of cern.colt.matrix.linalg.Algebra in project tetrad by cmu-phil.
the class SemUpdater method getUpdatedSemIm.
/**
* See http://en.wikipedia.org/wiki/Multivariate_normal_distribution.
*/
public SemIm getUpdatedSemIm() {
Algebra algebra = new Algebra();
// First manipulate the old semIm.
SemIm manipulatedSemIm = getManipulatedSemIm();
// Get out the means and implied covariances.
double[] means = new double[manipulatedSemIm.getVariableNodes().size()];
for (int i = 0; i < means.length; i++) {
means[i] = manipulatedSemIm.getMean(manipulatedSemIm.getVariableNodes().get(i));
}
DoubleMatrix1D mu = new DenseDoubleMatrix1D(means);
DoubleMatrix2D sigma = new DenseDoubleMatrix2D(manipulatedSemIm.getImplCovar(true).toArray());
// Updating on x2 = a.
SemEvidence evidence = getEvidence();
List nodesInEvidence = evidence.getNodesInEvidence();
int[] x2 = new int[nodesInEvidence.size()];
DoubleMatrix1D a = new DenseDoubleMatrix1D(nodesInEvidence.size());
for (int i = 0; i < nodesInEvidence.size(); i++) {
Node _node = (Node) nodesInEvidence.get(i);
x2[i] = evidence.getNodeIndex(_node);
}
for (int i = 0; i < nodesInEvidence.size(); i++) {
int j = evidence.getNodeIndex((Node) nodesInEvidence.get(i));
a.set(i, evidence.getProposition().getValue(j));
}
// x1 is all the variables.
// int[] x1 = new int[sigma.rows() - x2.length];
int[] x1 = new int[sigma.rows()];
for (int i = 0; i < sigma.rows(); i++) {
x1[i] = i;
}
// int index = -1;
// for (int i = 0; i < sigma.rows(); i++) {
// if (Arrays.binarySearch(x2, i) == -1) {
// x1[++index] = i;
// }
// }
// Calculate sigmaBar. (Don't know how to use it yet.)
// DoubleMatrix2D sigma11 = sigma.viewSelection(x1, x1);
DoubleMatrix2D sigma12 = sigma.viewSelection(x1, x2);
DoubleMatrix2D sigma22 = sigma.viewSelection(x2, x2);
// DoubleMatrix2D sigma21 = sigma.viewSelection(x2, x1);
DoubleMatrix2D inv_sigma22 = algebra.inverse(sigma22);
DoubleMatrix2D temp1 = algebra.mult(sigma12, inv_sigma22);
// DoubleMatrix2D temp2 = algebra.times(temp1, sigma21.copy());
// DoubleMatrix2D sigmaBar = sigma11.copy().assign(temp2, Functions.minus);
// Calculate muBar.
DoubleMatrix1D mu1 = mu.viewSelection(x1);
DoubleMatrix1D mu2 = mu.viewSelection(x2);
DoubleMatrix1D temp4 = a.copy().assign(mu2, Functions.minus);
DoubleMatrix1D temp5 = algebra.mult(temp1, temp4);
DoubleMatrix1D muBar = mu1.copy().assign(temp5, Functions.plus);
// Estimate a SEM with this sigmaBar and muBar.
// List variableNodes = manipulatedSemIm.getVariableNodes();
// String[] varNames = new String[variableNodes.size()];
//
// for (int i = 0; i < variableNodes.size(); i++) {
// varNames[i] = ((Node) variableNodes.get(i)).getNode();
// }
// System.out.println(sigmaBar);
//
// CovarianceMatrix covMatrix = new CovarianceMatrix(varNames,
// sigmaBar, 100);
// SemPm semPm = manipulatedSemIm.getEstIm();
// SemEstimator estimator = new SemEstimator(covMatrix, semPm);
// estimator.estimate();
// SemIm semIm = estimator.getEstimatedSem();
// semIm.setMeanValues(muBar.toArray());
// return semIm;
DoubleMatrix2D sigma2 = new DenseDoubleMatrix2D(manipulatedSemIm.getErrCovar().toArray());
// }
return manipulatedSemIm.updatedIm(new TetradMatrix(sigma2.toArray()), new TetradVector(muBar.toArray()));
}
use of cern.colt.matrix.linalg.Algebra in project tetrad by cmu-phil.
the class Ricf method ricf2.
/**
* same as above but takes a Graph instead of a SemGraph *
*/
public RicfResult ricf2(Graph mag, ICovarianceMatrix covMatrix, double tolerance) {
// mag.setShowErrorTerms(false);
DoubleFactory2D factory = DoubleFactory2D.dense;
Algebra algebra = new Algebra();
DoubleMatrix2D S = new DenseDoubleMatrix2D(covMatrix.getMatrix().toArray());
int p = covMatrix.getDimension();
if (p == 1) {
return new RicfResult(S, S, null, null, 1, Double.NaN, covMatrix);
}
List<Node> nodes = new ArrayList<>();
for (String name : covMatrix.getVariableNames()) {
nodes.add(mag.getNode(name));
}
DoubleMatrix2D omega = factory.diagonal(factory.diagonal(S));
DoubleMatrix2D B = factory.identity(p);
int[] ug = ugNodes(mag, nodes);
int[] ugComp = complement(p, ug);
if (ug.length > 0) {
List<Node> _ugNodes = new LinkedList<>();
for (int i : ug) {
_ugNodes.add(nodes.get(i));
}
Graph ugGraph = mag.subgraph(_ugNodes);
ICovarianceMatrix ugCov = covMatrix.getSubmatrix(ug);
DoubleMatrix2D lambdaInv = fitConGraph(ugGraph, ugCov, p + 1, tolerance).shat;
omega.viewSelection(ug, ug).assign(lambdaInv);
}
// Prepare lists of parents and spouses.
int[][] pars = parentIndices(p, mag, nodes);
int[][] spo = spouseIndices(p, mag, nodes);
int i = 0;
double _diff;
while (true) {
i++;
DoubleMatrix2D omegaOld = omega.copy();
DoubleMatrix2D bOld = B.copy();
for (int _v = 0; _v < p; _v++) {
// Exclude the UG part.
if (Arrays.binarySearch(ug, _v) >= 0) {
continue;
}
int[] v = new int[] { _v };
int[] vcomp = complement(p, v);
int[] all = range(0, p - 1);
int[] parv = pars[_v];
int[] spov = spo[_v];
DoubleMatrix2D a6 = B.viewSelection(v, parv);
if (spov.length == 0) {
if (parv.length != 0) {
if (i == 1) {
DoubleMatrix2D a1 = S.viewSelection(parv, parv);
DoubleMatrix2D a2 = S.viewSelection(v, parv);
DoubleMatrix2D a3 = algebra.inverse(a1);
DoubleMatrix2D a4 = algebra.mult(a2, a3);
a4.assign(Mult.mult(-1));
a6.assign(a4);
DoubleMatrix2D a7 = S.viewSelection(parv, v);
DoubleMatrix2D a9 = algebra.mult(a6, a7);
DoubleMatrix2D a8 = S.viewSelection(v, v);
DoubleMatrix2D a8b = omega.viewSelection(v, v);
a8b.assign(a8);
omega.viewSelection(v, v).assign(a9, PlusMult.plusMult(1));
}
}
} else {
if (parv.length != 0) {
DoubleMatrix2D oInv = new DenseDoubleMatrix2D(p, p);
DoubleMatrix2D a2 = omega.viewSelection(vcomp, vcomp);
DoubleMatrix2D a3 = algebra.inverse(a2);
oInv.viewSelection(vcomp, vcomp).assign(a3);
DoubleMatrix2D Z = algebra.mult(oInv.viewSelection(spov, vcomp), B.viewSelection(vcomp, all));
int lpa = parv.length;
int lspo = spov.length;
// Build XX
DoubleMatrix2D XX = new DenseDoubleMatrix2D(lpa + lspo, lpa + lspo);
int[] range1 = range(0, lpa - 1);
int[] range2 = range(lpa, lpa + lspo - 1);
// Upper left quadrant
XX.viewSelection(range1, range1).assign(S.viewSelection(parv, parv));
// Upper right quadrant
DoubleMatrix2D a11 = algebra.mult(S.viewSelection(parv, all), algebra.transpose(Z));
XX.viewSelection(range1, range2).assign(a11);
// Lower left quadrant
DoubleMatrix2D a12 = XX.viewSelection(range2, range1);
DoubleMatrix2D a13 = algebra.transpose(XX.viewSelection(range1, range2));
a12.assign(a13);
// Lower right quadrant
DoubleMatrix2D a14 = XX.viewSelection(range2, range2);
DoubleMatrix2D a15 = algebra.mult(Z, S);
DoubleMatrix2D a16 = algebra.mult(a15, algebra.transpose(Z));
a14.assign(a16);
// Build XY
DoubleMatrix1D YX = new DenseDoubleMatrix1D(lpa + lspo);
DoubleMatrix1D a17 = YX.viewSelection(range1);
DoubleMatrix1D a18 = S.viewSelection(v, parv).viewRow(0);
a17.assign(a18);
DoubleMatrix1D a19 = YX.viewSelection(range2);
DoubleMatrix2D a20 = S.viewSelection(v, all);
DoubleMatrix1D a21 = algebra.mult(a20, algebra.transpose(Z)).viewRow(0);
a19.assign(a21);
// Temp
DoubleMatrix2D a22 = algebra.inverse(XX);
DoubleMatrix1D temp = algebra.mult(algebra.transpose(a22), YX);
// Assign to b.
DoubleMatrix1D a23 = a6.viewRow(0);
DoubleMatrix1D a24 = temp.viewSelection(range1);
a23.assign(a24);
a23.assign(Mult.mult(-1));
// Assign to omega.
omega.viewSelection(v, spov).viewRow(0).assign(temp.viewSelection(range2));
omega.viewSelection(spov, v).viewColumn(0).assign(temp.viewSelection(range2));
// Variance.
double tempVar = S.get(_v, _v) - algebra.mult(temp, YX);
DoubleMatrix2D a27 = omega.viewSelection(v, spov);
DoubleMatrix2D a28 = oInv.viewSelection(spov, spov);
DoubleMatrix2D a29 = omega.viewSelection(spov, v).copy();
DoubleMatrix2D a30 = algebra.mult(a27, a28);
DoubleMatrix2D a31 = algebra.mult(a30, a29);
omega.viewSelection(v, v).assign(tempVar);
omega.viewSelection(v, v).assign(a31, PlusMult.plusMult(1));
} else {
DoubleMatrix2D oInv = new DenseDoubleMatrix2D(p, p);
DoubleMatrix2D a2 = omega.viewSelection(vcomp, vcomp);
DoubleMatrix2D a3 = algebra.inverse(a2);
oInv.viewSelection(vcomp, vcomp).assign(a3);
// System.out.println("O.inv = " + oInv);
DoubleMatrix2D a4 = oInv.viewSelection(spov, vcomp);
DoubleMatrix2D a5 = B.viewSelection(vcomp, all);
DoubleMatrix2D Z = algebra.mult(a4, a5);
// System.out.println("Z = " + Z);
// Build XX
DoubleMatrix2D XX = algebra.mult(algebra.mult(Z, S), Z.viewDice());
// System.out.println("XX = " + XX);
// Build XY
DoubleMatrix2D a20 = S.viewSelection(v, all);
DoubleMatrix1D YX = algebra.mult(a20, Z.viewDice()).viewRow(0);
// System.out.println("YX = " + YX);
// Temp
DoubleMatrix2D a22 = algebra.inverse(XX);
DoubleMatrix1D a23 = algebra.mult(algebra.transpose(a22), YX);
// Assign to omega.
DoubleMatrix1D a24 = omega.viewSelection(v, spov).viewRow(0);
a24.assign(a23);
DoubleMatrix1D a25 = omega.viewSelection(spov, v).viewColumn(0);
a25.assign(a23);
// System.out.println("Omega 2 " + omega);
// Variance.
double tempVar = S.get(_v, _v) - algebra.mult(a24, YX);
// System.out.println("tempVar = " + tempVar);
DoubleMatrix2D a27 = omega.viewSelection(v, spov);
DoubleMatrix2D a28 = oInv.viewSelection(spov, spov);
DoubleMatrix2D a29 = omega.viewSelection(spov, v).copy();
DoubleMatrix2D a30 = algebra.mult(a27, a28);
DoubleMatrix2D a31 = algebra.mult(a30, a29);
omega.set(_v, _v, tempVar + a31.get(0, 0));
// System.out.println("Omega final " + omega);
}
}
}
DoubleMatrix2D a32 = omega.copy();
a32.assign(omegaOld, PlusMult.plusMult(-1));
double diff1 = algebra.norm1(a32);
DoubleMatrix2D a33 = B.copy();
a33.assign(bOld, PlusMult.plusMult(-1));
double diff2 = algebra.norm1(a32);
double diff = diff1 + diff2;
_diff = diff;
if (diff < tolerance)
break;
}
DoubleMatrix2D a34 = algebra.inverse(B);
DoubleMatrix2D a35 = algebra.inverse(B.viewDice());
DoubleMatrix2D sigmahat = algebra.mult(algebra.mult(a34, omega), a35);
DoubleMatrix2D lambdahat = omega.copy();
DoubleMatrix2D a36 = lambdahat.viewSelection(ugComp, ugComp);
a36.assign(factory.make(ugComp.length, ugComp.length, 0.0));
DoubleMatrix2D omegahat = omega.copy();
DoubleMatrix2D a37 = omegahat.viewSelection(ug, ug);
a37.assign(factory.make(ug.length, ug.length, 0.0));
DoubleMatrix2D bhat = B.copy();
return new RicfResult(sigmahat, lambdahat, bhat, omegahat, i, _diff, covMatrix);
}
Aggregations