use of cern.colt.matrix.DoubleFactory2D in project tetrad by cmu-phil.
the class Ricf method fitConGraph.
/**
* Fits a concentration graph. Coding algorithm #2 only.
*/
private FitConGraphResult fitConGraph(Graph graph, ICovarianceMatrix cov, int n, double tol) {
DoubleFactory2D factory = DoubleFactory2D.dense;
Algebra algebra = new Algebra();
List<Node> nodes = graph.getNodes();
String[] nodeNames = new String[nodes.size()];
for (int i = 0; i < nodes.size(); i++) {
Node node = nodes.get(i);
if (!cov.getVariableNames().contains(node.getName())) {
throw new IllegalArgumentException("Node in graph not in cov matrix: " + node);
}
nodeNames[i] = node.getName();
}
DoubleMatrix2D S = new DenseDoubleMatrix2D(cov.getSubmatrix(nodeNames).getMatrix().toArray());
graph = graph.subgraph(nodes);
List<List<Node>> cli = cliques(graph);
int nc = cli.size();
if (nc == 1) {
return new FitConGraphResult(S, 0, 0, 1);
}
int k = S.rows();
int it = 0;
// Only coding alg #2 here.
DoubleMatrix2D K = algebra.inverse(factory.diagonal(factory.diagonal(S)));
int[] all = range(0, k - 1);
while (true) {
DoubleMatrix2D KOld = K.copy();
it++;
for (List<Node> aCli : cli) {
int[] a = asIndices(aCli, nodes);
int[] b = complement(all, a);
DoubleMatrix2D a1 = S.viewSelection(a, a);
DoubleMatrix2D a2 = algebra.inverse(a1);
DoubleMatrix2D a3 = K.viewSelection(a, b);
DoubleMatrix2D a4 = K.viewSelection(b, b);
DoubleMatrix2D a5 = algebra.inverse(a4);
DoubleMatrix2D a6 = K.viewSelection(b, a).copy();
DoubleMatrix2D a7 = algebra.mult(a3, a5);
DoubleMatrix2D a8 = algebra.mult(a7, a6);
a2.assign(a8, PlusMult.plusMult(1));
DoubleMatrix2D a9 = K.viewSelection(a, a);
a9.assign(a2);
}
DoubleMatrix2D a32 = K.copy();
a32.assign(KOld, PlusMult.plusMult(-1));
double diff = algebra.norm1(a32);
if (diff < tol)
break;
}
DoubleMatrix2D V = algebra.inverse(K);
int numNodes = graph.getNumNodes();
int df = numNodes * (numNodes - 1) / 2 - graph.getNumEdges();
double dev = lik(algebra.inverse(V), S, n, k);
return new FitConGraphResult(V, dev, df, it);
}
use of cern.colt.matrix.DoubleFactory2D in project tetrad by cmu-phil.
the class Ricf method ricf.
// =============================PUBLIC METHODS=========================//
public RicfResult ricf(SemGraph 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);
}
use of cern.colt.matrix.DoubleFactory2D in project tdq-studio-se by Talend.
the class TestMatrix2D method doubleTest25.
/**
*/
public static void doubleTest25(int size) {
System.out.println("\n\n");
System.out.println("initializing...");
boolean dense = true;
DoubleMatrix2D A;
DoubleFactory2D factory;
if (dense)
factory = Factory2D.dense;
else
factory = Factory2D.sparse;
double value = 0.5;
A = factory.make(size, size, value);
Property.generateNonSingular(A);
cern.colt.Timer timer = new cern.colt.Timer().start();
System.out.println(A);
System.out.println(Algebra.ZERO.inverse(A));
timer.stop().display();
System.out.println("done.");
}
use of cern.colt.matrix.DoubleFactory2D in project tdq-studio-se by Talend.
the class QRTest method main.
public static void main(String[] args) {
// For COLT
DoubleMatrix2D xmatrix, ymatrix, zmatrix;
DoubleFactory2D myfactory;
myfactory = DoubleFactory2D.dense;
xmatrix = myfactory.make(8, 2);
ymatrix = myfactory.make(8, 1);
xmatrix.set(0, 0, 1);
xmatrix.set(1, 0, 1);
xmatrix.set(2, 0, 1);
xmatrix.set(3, 0, 1);
xmatrix.set(4, 0, 1);
xmatrix.set(5, 0, 1);
xmatrix.set(6, 0, 1);
xmatrix.set(7, 0, 1);
xmatrix.set(0, 1, 80);
xmatrix.set(1, 1, 220);
xmatrix.set(2, 1, 140);
xmatrix.set(3, 1, 120);
xmatrix.set(4, 1, 180);
xmatrix.set(5, 1, 100);
xmatrix.set(6, 1, 200);
xmatrix.set(7, 1, 160);
ymatrix.set(0, 0, 0.6);
ymatrix.set(1, 0, 6.70);
ymatrix.set(2, 0, 5.30);
ymatrix.set(3, 0, 4.00);
ymatrix.set(4, 0, 6.55);
ymatrix.set(5, 0, 2.15);
ymatrix.set(6, 0, 6.60);
ymatrix.set(7, 0, 5.75);
Algebra myAlgebra = new Algebra();
zmatrix = myAlgebra.solve(xmatrix, ymatrix);
System.err.println(xmatrix);
System.err.println(ymatrix);
System.err.println(zmatrix);
/*
// For JAMA
Matrix amatrix,bmatrix,cmatrix;
amatrix = new Matrix(8,2);
bmatrix = new Matrix(8,1);
amatrix.set(0,0,1);
amatrix.set(1,0,1);
amatrix.set(2,0,1);
amatrix.set(3,0,1);
amatrix.set(4,0,1);
amatrix.set(5,0,1);
amatrix.set(6,0,1);
amatrix.set(7,0,1);
amatrix.set(0,1,80);
amatrix.set(1,1,220);
amatrix.set(2,1,140);
amatrix.set(3,1,120);
amatrix.set(4,1,180);
amatrix.set(5,1,100);
amatrix.set(6,1,200);
amatrix.set(7,1,160);
bmatrix.set(0,0,0.6);
bmatrix.set(1,0,6.70);
bmatrix.set(2,0,5.30);
bmatrix.set(3,0,4.00);
bmatrix.set(4,0,6.55);
bmatrix.set(5,0,2.15);
bmatrix.set(6,0,6.60);
bmatrix.set(7,0,5.75);
cmatrix = amatrix.solve(bmatrix);
amatrix.print(8,5);
bmatrix.print(8,5);
cmatrix.print(8,5);
*/
}
use of cern.colt.matrix.DoubleFactory2D in project tdq-studio-se by Talend.
the class Statistic method demo3.
/**
* Demonstrates usage of this class.
*/
public static void demo3(VectorVectorFunction norm) {
double[][] values = { { -0.9611052, -0.25421095 }, { 0.4308269, -0.69932648 }, { -1.2071029, 0.62030596 }, { 1.5345166, 0.02135884 }, { -1.1341542, 0.20388430 } };
System.out.println("\n\ninitializing...");
DoubleFactory2D factory = DoubleFactory2D.dense;
DoubleMatrix2D A = factory.make(values).viewDice();
System.out.println("\nA=" + A.viewDice());
System.out.println("\ndist=" + distance(A, norm).viewDice());
}
Aggregations