Search in sources :

Example 11 with DenseDoubleMatrix1D

use of cern.colt.matrix.impl.DenseDoubleMatrix1D 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);
}
Also used : ICovarianceMatrix(edu.cmu.tetrad.data.ICovarianceMatrix) Node(edu.cmu.tetrad.graph.Node) DoubleFactory2D(cern.colt.matrix.DoubleFactory2D) Endpoint(edu.cmu.tetrad.graph.Endpoint) Algebra(cern.colt.matrix.linalg.Algebra) SemGraph(edu.cmu.tetrad.graph.SemGraph) Graph(edu.cmu.tetrad.graph.Graph) DoubleMatrix2D(cern.colt.matrix.DoubleMatrix2D) DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D) DoubleMatrix1D(cern.colt.matrix.DoubleMatrix1D) DenseDoubleMatrix1D(cern.colt.matrix.impl.DenseDoubleMatrix1D) DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D) DenseDoubleMatrix1D(cern.colt.matrix.impl.DenseDoubleMatrix1D)

Example 12 with DenseDoubleMatrix1D

use of cern.colt.matrix.impl.DenseDoubleMatrix1D in project tetrad by cmu-phil.

the class Glasso method search.

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));
        }
    }
    // (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);
    }
    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);
            }
        }
    }
    for (int j = 0; j < n; j++) {
        if (pen) {
            ww.set(j, j, ss.get(j, j) + rho.get(j, j));
        } else {
            ww.set(j, j, ss.get(j, j));
        // System.out.println(ww);
        }
    }
    niter = 0;
    while (true) {
        dlx = 0.0;
        for (int m = 0; m < n; m++) {
            if (itr) {
                System.out.println("Outer loop = " + m);
            }
            x = xs.viewColumn(m);
            // System.out.println(x);
            // System.out.println();
            ws = ww.viewColumn(m);
            // This sets up vv, s, and ro--i.e., W.11, s.12, and r.12.
            setup(m, n, ss, rho, ww, vv, s, ro);
            so.assign(s);
            // System.out.println("ww = " + ww);
            // This updates s and x--the estimated correlation matrix and the reduced form of the
            // estimated inverse covariance.
            lasso(ro, nm1, vv, s, shr / sum_abs(vv), x, z, mm);
            // lasso(ro,nm1,vv,s,thr/sum_abs(vv),x,z,mm);
            int l = -1;
            for (int j = 0; j < n; j++) {
                if (Thread.currentThread().isInterrupted()) {
                    break;
                }
                if (j == m)
                    continue;
                l = l + 1;
                ww.set(j, m, so.get(l) - s.get(l));
                ww.set(m, j, ww.get(j, m));
            }
            dlx = Math.max(dlx, sum_abs_diff(ww.viewColumn(m), ws));
            // xs(:,m)=x
            xs.viewColumn(m).assign(x);
        }
        niter = niter + 1;
        if (niter < getMaxit())
            break;
        if (dlx < shr)
            break;
    }
    del = dlx / nm1;
    inv(n, ww, xs, wwi);
    return new Result(ww, wwi, niter, del);
}
Also used : DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D) DoubleMatrix2D(cern.colt.matrix.DoubleMatrix2D) DoubleMatrix1D(cern.colt.matrix.DoubleMatrix1D) DenseDoubleMatrix1D(cern.colt.matrix.impl.DenseDoubleMatrix1D) DenseDoubleMatrix2D(cern.colt.matrix.impl.DenseDoubleMatrix2D) DenseDoubleMatrix1D(cern.colt.matrix.impl.DenseDoubleMatrix1D)

Aggregations

DoubleMatrix1D (cern.colt.matrix.DoubleMatrix1D)12 DenseDoubleMatrix1D (cern.colt.matrix.impl.DenseDoubleMatrix1D)12 DoubleMatrix2D (cern.colt.matrix.DoubleMatrix2D)6 DenseDoubleMatrix2D (cern.colt.matrix.impl.DenseDoubleMatrix2D)6 Algebra (cern.colt.matrix.linalg.Algebra)5 Node (edu.cmu.tetrad.graph.Node)3 DoubleArrayList (cern.colt.list.DoubleArrayList)2 DoubleFactory2D (cern.colt.matrix.DoubleFactory2D)2 ICovarianceMatrix (edu.cmu.tetrad.data.ICovarianceMatrix)2 Endpoint (edu.cmu.tetrad.graph.Endpoint)2 Graph (edu.cmu.tetrad.graph.Graph)2 SemGraph (edu.cmu.tetrad.graph.SemGraph)2 AbstractFormatter (cern.colt.matrix.impl.AbstractFormatter)1 TetradMatrix (edu.cmu.tetrad.util.TetradMatrix)1 TetradVector (edu.cmu.tetrad.util.TetradVector)1 ArrayList (java.util.ArrayList)1 List (java.util.List)1 StopWatch (org.apache.commons.lang3.time.StopWatch)1 ByteArrayConverter (ubic.basecode.io.ByteArrayConverter)1 ExpressionDataDoubleMatrix (ubic.gemma.core.datastructure.matrix.ExpressionDataDoubleMatrix)1