Example usage for org.apache.commons.math.linear RealMatrix transpose

List of usage examples for org.apache.commons.math.linear RealMatrix transpose

Introduction

In this page you can find the example usage for org.apache.commons.math.linear RealMatrix transpose.

Prototype

RealMatrix transpose();

Source Link

Document

Returns the transpose of this matrix.

Usage

From source file:math2605.gn_exp.java

private static RealMatrix findQR(RealMatrix J, RealMatrix r) {
    qr_fact_househ m = new qr_fact_househ(J);
    RealMatrix Q = m.getQ();
    RealMatrix qTranspose = Q.transpose();
    System.out.println("qtranspose");
    System.out.println(qTranspose);
    RealMatrix R = m.getR();/*ww  w .ja  v a 2s. c o m*/
    MatrixMethods temp = new MatrixMethods(R);
    RealMatrix rInverse = temp.inverseMatrix();
    System.out.println("rInverse");
    System.out.println(rInverse);
    RealMatrix sub = rInverse.multiply(qTranspose).multiply(r);
    return sub;
}

From source file:math2605.qr_fact_givens.java

public qr_fact_givens(RealMatrix A) {
    int m = A.getRowDimension();
    int n = A.getColumnDimension();
    R = A.copy();/*from  ww  w . j  ava 2 s. c  o m*/
    RealMatrix a = A.copy();
    Q = null;
    boolean firstQ = true;
    for (int j = 0; j < n; j++) {
        int row = j;
        for (int i = j; i < m - 1; i++) {
            if (a.getEntry(row, j) == 0) {
                //System.out.println("in break");
                break;
            }
            if (upperTriangular(a)) {
                //System.out.println("is ut");
                continue;
            }

            double[] cs = new double[2];

            cs = givensRoatation(a.getEntry(i, j), a.getEntry(i + 1, j));
            RealMatrix G = createG(i, i + 1, cs, n);
            R = G.multiply(R);
            if (firstQ) {
                Q = G.transpose();
                firstQ = false;
            } else {
                Q = Q.multiply(G.transpose());
            }
            a = G.multiply(a);
            row++;
        }
    }
}

From source file:com.opengamma.analytics.math.matrix.CommonsMatrixAlgebra.java

/**
 * {@inheritDoc}/*from  ww w  .j  a  va  2  s.  c om*/
 */
@Override
public DoubleMatrix2D getTranspose(final Matrix<?> m) {
    Validate.notNull(m, "m");
    if (m instanceof DoubleMatrix2D) {
        final RealMatrix temp = CommonsMathWrapper.wrap((DoubleMatrix2D) m);
        return CommonsMathWrapper.unwrap(temp.transpose());
    }
    throw new IllegalArgumentException("Can only find transpose of DoubleMatrix2D; have " + m.getClass());
}

From source file:jml.matlab.utils.SingularValueDecompositionImpl.java

/** {@inheritDoc} */
public RealMatrix getCovariance(final double minSingularValue) {
    // get the number of singular values to consider
    final int p = s.length;
    int dimension = 0;
    while ((dimension < p) && (s[dimension] >= minSingularValue)) {
        ++dimension;/*from   w ww . ja v a  2  s. co  m*/
    }

    if (dimension == 0) {
        throw new NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE, minSingularValue,
                s[0], true);
    }

    final double[][] data = new double[dimension][p];
    getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {

        /** {@inheritDoc} */
        @Override
        public void visit(final int row, final int column, final double value) {
            data[row][column] = value / s[row];
        }
    }, 0, dimension - 1, 0, p - 1);

    RealMatrix jv = new Array2DRowRealMatrix(data, false);
    return jv.transpose().multiply(jv);
}

From source file:lib.regressions.MultipleRegression.java

/**
 * Perform the regression computations/* w w w. j  ava2 s.co m*/
 */
private void compute() {
    // Set everything to 0.
    for (int i = 0; i < (myNumVar + 1); i++) {
        myCoef[i] = 0.0;
        myStdErr[i] = 0.0;
        myTStat[i] = 0.0;
    }
    myChiSq = 0.0;
    myRSq = 0.0;
    myAdjustedRSq = 0.0;

    // Set coefficients, t-stat, etc. if there has been enough data added.
    if (myCount >= (myNumVar + 1)) {
        RealMatrix dataMatrix = new RealMatrixImpl(mySums.getSumXX());
        RealMatrix xxMatrix = dataMatrix.getSubMatrix(1, myNumVar + 1, 1, myNumVar + 1);
        RealMatrix xyMatrix = dataMatrix.getSubMatrix(1, myNumVar + 1, 0, 0);

        computeOkX(); // Determine which X components to use.
        int[] listX = getListX();
        int[] listY = { 0 };
        int numX = listX.length;
        RealMatrix xxSubMatrix = xxMatrix.getSubMatrix(listX, listX);
        RealMatrix xySubMatrix = xyMatrix.getSubMatrix(listX, listY);

        double sumY = mySums.getSumXX()[0][1];
        double sumYY = mySums.getSumXX()[0][0];

        if (!xxSubMatrix.isSingular()) {
            RealMatrix xxInverse = xxSubMatrix.inverse();
            RealMatrix coefMatrix = xxInverse.multiply(xySubMatrix);
            double[] coef = coefMatrix.getColumn(0);

            // Compute chi-squared
            myChiSq = sumYY - 2 * coefMatrix.transpose().multiply(xySubMatrix).getEntry(0, 0)
                    + +coefMatrix.transpose().multiply(xxSubMatrix).multiply(coefMatrix).getEntry(0, 0);

            // Compute R^2 and adjusted R^2
            int offset = getUseIntercept() ? 1 : 0;
            myRSq = 1.0 - myChiSq / (sumYY - sumY * sumY / myCount);
            myAdjustedRSq = 1 - ((1 - myRSq) * (myCount - 1) + 1 - offset) / (myCount - numX);

            // Compute standard errors and t-stats
            int j = 0;
            for (int i = 0; i < (myNumVar + 1); i++) {
                if (myOkX[i]) {
                    j++;
                    myCoef[i] = coef[j - 1];
                    myStdErr[i] = Math.sqrt(myChiSq * xxInverse.getEntry(j - 1, j - 1) / (myCount - numX));
                    myTStat[i] = myCoef[i] / myStdErr[i];
                }
            }
        }

    }
    myIsComputed = true;
}

From source file:jml.matlab.utils.SingularValueDecompositionImpl.java

/**
 * Calculates the compact Singular Value Decomposition of the given matrix.
 *
 * @param matrix Matrix to decompose.//from w  w w. ja  va2s  .  c  om
 */
public SingularValueDecompositionImpl(RealMatrix matrix) {
    double[][] U, V;

    // Derived from LINPACK code.
    // Initialize.
    double[][] A;
    m = matrix.getRowDimension();
    n = matrix.getColumnDimension();

    if (matrix.getRowDimension() < matrix.getColumnDimension()) {
        transposed = true;
        A = matrix.transpose().getData();
        m = matrix.getColumnDimension();
        n = matrix.getRowDimension();
    } else {
        transposed = false;
        A = matrix.getData();
        m = matrix.getRowDimension();
        n = matrix.getColumnDimension();
    }

    int nu = Math.min(m, n);
    s = new double[Math.min(m + 1, n)];
    U = new double[m][nu];
    V = new double[n][n];
    double[] e = new double[n];
    double[] work = new double[m];
    boolean wantu = true;
    boolean wantv = true;

    // Reduce A to bidiagonal form, storing the diagonal elements
    // in s and the super-diagonal elements in e.
    int nct = Math.min(m - 1, n);
    int nrt = Math.max(0, Math.min(n - 2, m));
    for (int k = 0; k < Math.max(nct, nrt); k++) {
        if (k < nct) {

            // Compute the transformation for the k-th column and
            // place the k-th diagonal in s[k].
            // Compute 2-norm of k-th column without under/overflow.
            s[k] = 0;
            for (int i = k; i < m; i++) {
                s[k] = hypot(s[k], A[i][k]);
            }
            if (s[k] != 0.0) {
                if (A[k][k] < 0.0) {
                    s[k] = -s[k];
                }
                for (int i = k; i < m; i++) {
                    A[i][k] /= s[k];
                }
                A[k][k] += 1.0;
            }
            s[k] = -s[k];
        }
        for (int j = k + 1; j < n; j++) {
            if ((k < nct) & (s[k] != 0.0)) {

                // Apply the transformation.

                double t = 0;
                for (int i = k; i < m; i++) {
                    t += A[i][k] * A[i][j];
                }
                t = -t / A[k][k];
                for (int i = k; i < m; i++) {
                    A[i][j] += t * A[i][k];
                }
            }

            // Place the k-th row of A into e for the
            // subsequent calculation of the row transformation.
            e[j] = A[k][j];
        }
        if (wantu & (k < nct)) {

            // Place the transformation in U for subsequent back
            // multiplication.

            for (int i = k; i < m; i++) {
                U[i][k] = A[i][k];
            }
        }
        if (k < nrt) {

            // Compute the k-th row transformation and place the
            // k-th super-diagonal in e[k].
            // Compute 2-norm without under/overflow.
            e[k] = 0;
            for (int i = k + 1; i < n; i++) {
                e[k] = hypot(e[k], e[i]);
            }
            if (e[k] != 0.0) {
                if (e[k + 1] < 0.0) {
                    e[k] = -e[k];
                }
                for (int i = k + 1; i < n; i++) {
                    e[i] /= e[k];
                }
                e[k + 1] += 1.0;
            }
            e[k] = -e[k];
            if ((k + 1 < m) & (e[k] != 0.0)) {

                // Apply the transformation.

                for (int i = k + 1; i < m; i++) {
                    work[i] = 0.0;
                }
                for (int j = k + 1; j < n; j++) {
                    for (int i = k + 1; i < m; i++) {
                        work[i] += e[j] * A[i][j];
                    }
                }
                for (int j = k + 1; j < n; j++) {
                    double t = -e[j] / e[k + 1];
                    for (int i = k + 1; i < m; i++) {
                        A[i][j] += t * work[i];
                    }
                }
            }
            if (wantv) {

                // Place the transformation in V for subsequent
                // back multiplication.
                for (int i = k + 1; i < n; i++) {
                    V[i][k] = e[i];
                }
            }
        }
    }

    // Set up the final bidiagonal matrix or order p.
    int p = Math.min(n, m + 1);
    if (nct < n) {
        s[nct] = A[nct][nct];
    }
    if (m < p) {
        s[p - 1] = 0.0;
    }
    if (nrt + 1 < p) {
        e[nrt] = A[nrt][p - 1];
    }
    e[p - 1] = 0.0;

    // If required, generate U.
    if (wantu) {
        for (int j = nct; j < nu; j++) {
            for (int i = 0; i < m; i++) {
                U[i][j] = 0.0;
            }
            U[j][j] = 1.0;
        }
        for (int k = nct - 1; k >= 0; k--) {
            if (s[k] != 0.0) {
                for (int j = k + 1; j < nu; j++) {
                    double t = 0;
                    for (int i = k; i < m; i++) {
                        t += U[i][k] * U[i][j];
                    }
                    t = -t / U[k][k];
                    for (int i = k; i < m; i++) {
                        U[i][j] += t * U[i][k];
                    }
                }
                for (int i = k; i < m; i++) {
                    U[i][k] = -U[i][k];
                }
                U[k][k] = 1.0 + U[k][k];
                for (int i = 0; i < k - 1; i++) {
                    U[i][k] = 0.0;
                }
            } else {
                for (int i = 0; i < m; i++) {
                    U[i][k] = 0.0;
                }
                U[k][k] = 1.0;
            }
        }
    }

    // If required, generate V.
    if (wantv) {
        for (int k = n - 1; k >= 0; k--) {
            if ((k < nrt) & (e[k] != 0.0)) {
                for (int j = k + 1; j < nu; j++) {
                    double t = 0;
                    for (int i = k + 1; i < n; i++) {
                        t += V[i][k] * V[i][j];
                    }
                    t = -t / V[k + 1][k];
                    for (int i = k + 1; i < n; i++) {
                        V[i][j] += t * V[i][k];
                    }
                }
            }
            for (int i = 0; i < n; i++) {
                V[i][k] = 0.0;
            }
            V[k][k] = 1.0;
        }
    }

    // Main iteration loop for the singular values.
    int pp = p - 1;
    int iter = 0;
    double eps = Math.pow(2.0, -52.0);
    double tiny = Math.pow(2.0, -966.0);
    while (p > 0) {
        int k, kase;

        // Here is where a test for too many iterations would go.

        // This section of the program inspects for
        // negligible elements in the s and e arrays.  On
        // completion the variables kase and k are set as follows.

        // kase = 1     if s(p) and e[k-1] are negligible and k<p
        // kase = 2     if s(k) is negligible and k<p
        // kase = 3     if e[k-1] is negligible, k<p, and
        //              s(k), ..., s(p) are not negligible (qr step).
        // kase = 4     if e(p-1) is negligible (convergence).

        for (k = p - 2; k >= -1; k--) {
            if (k == -1) {
                break;
            }
            if (Math.abs(e[k]) <= tiny + eps * (Math.abs(s[k]) + Math.abs(s[k + 1]))) {
                e[k] = 0.0;
                break;
            }
        }
        if (k == p - 2) {
            kase = 4;
        } else {
            int ks;
            for (ks = p - 1; ks >= k; ks--) {
                if (ks == k) {
                    break;
                }
                double t = (ks != p ? Math.abs(e[ks]) : 0.) + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.);
                if (Math.abs(s[ks]) <= tiny + eps * t) {
                    s[ks] = 0.0;
                    break;
                }
            }
            if (ks == k) {
                kase = 3;
            } else if (ks == p - 1) {
                kase = 1;
            } else {
                kase = 2;
                k = ks;
            }
        }
        k++;

        // Perform the task indicated by kase.

        switch (kase) {

        // Deflate negligible s(p).

        case 1: {
            double f = e[p - 2];
            e[p - 2] = 0.0;
            for (int j = p - 2; j >= k; j--) {
                double t = hypot(s[j], f);
                double cs = s[j] / t;
                double sn = f / t;
                s[j] = t;
                if (j != k) {
                    f = -sn * e[j - 1];
                    e[j - 1] = cs * e[j - 1];
                }
                if (wantv) {
                    for (int i = 0; i < n; i++) {
                        t = cs * V[i][j] + sn * V[i][p - 1];
                        V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
                        V[i][j] = t;
                    }
                }
            }
        }
            break;

        // Split at negligible s(k).

        case 2: {
            double f = e[k - 1];
            e[k - 1] = 0.0;
            for (int j = k; j < p; j++) {
                double t = hypot(s[j], f);
                double cs = s[j] / t;
                double sn = f / t;
                s[j] = t;
                f = -sn * e[j];
                e[j] = cs * e[j];
                if (wantu) {
                    for (int i = 0; i < m; i++) {
                        t = cs * U[i][j] + sn * U[i][k - 1];
                        U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
                        U[i][j] = t;
                    }
                }
            }
        }
            break;

        // Perform one qr step.

        case 3: {

            // Calculate the shift.

            double scale = Math.max(
                    Math.max(Math.max(Math.max(Math.abs(s[p - 1]), Math.abs(s[p - 2])), Math.abs(e[p - 2])),
                            Math.abs(s[k])),
                    Math.abs(e[k]));
            double sp = s[p - 1] / scale;
            double spm1 = s[p - 2] / scale;
            double epm1 = e[p - 2] / scale;
            double sk = s[k] / scale;
            double ek = e[k] / scale;
            double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
            double c = (sp * epm1) * (sp * epm1);
            double shift = 0.0;
            if ((b != 0.0) | (c != 0.0)) {
                shift = Math.sqrt(b * b + c);
                if (b < 0.0) {
                    shift = -shift;
                }
                shift = c / (b + shift);
            }
            double f = (sk + sp) * (sk - sp) + shift;
            double g = sk * ek;

            // Chase zeros.

            for (int j = k; j < p - 1; j++) {
                double t = hypot(f, g);
                double cs = f / t;
                double sn = g / t;
                if (j != k) {
                    e[j - 1] = t;
                }
                f = cs * s[j] + sn * e[j];
                e[j] = cs * e[j] - sn * s[j];
                g = sn * s[j + 1];
                s[j + 1] = cs * s[j + 1];
                if (wantv) {
                    for (int i = 0; i < n; i++) {
                        t = cs * V[i][j] + sn * V[i][j + 1];
                        V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
                        V[i][j] = t;
                    }
                }
                t = hypot(f, g);
                cs = f / t;
                sn = g / t;
                s[j] = t;
                f = cs * e[j] + sn * s[j + 1];
                s[j + 1] = -sn * e[j] + cs * s[j + 1];
                g = sn * e[j + 1];
                e[j + 1] = cs * e[j + 1];
                if (wantu && (j < m - 1)) {
                    for (int i = 0; i < m; i++) {
                        t = cs * U[i][j] + sn * U[i][j + 1];
                        U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
                        U[i][j] = t;
                    }
                }
            }
            e[p - 2] = f;
            iter = iter + 1;
        }
            break;

        // Convergence.

        case 4: {

            // Make the singular values positive.

            if (s[k] <= 0.0) {
                s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
                if (wantv) {
                    for (int i = 0; i <= pp; i++) {
                        V[i][k] = -V[i][k];
                    }
                }
            }

            // Order the singular values.

            while (k < pp) {
                if (s[k] >= s[k + 1]) {
                    break;
                }
                double t = s[k];
                s[k] = s[k + 1];
                s[k + 1] = t;
                if (wantv && (k < n - 1)) {
                    for (int i = 0; i < n; i++) {
                        t = V[i][k + 1];
                        V[i][k + 1] = V[i][k];
                        V[i][k] = t;
                    }
                }
                if (wantu && (k < m - 1)) {
                    for (int i = 0; i < m; i++) {
                        t = U[i][k + 1];
                        U[i][k + 1] = U[i][k];
                        U[i][k] = t;
                    }
                }
                k++;
            }
            iter = 0;
            p--;
        }
            break;
        }
    }

    if (!transposed) {
        cachedU = MatrixUtils.createRealMatrix(U);
        cachedV = MatrixUtils.createRealMatrix(V);
    } else {
        cachedU = MatrixUtils.createRealMatrix(V);
        cachedV = MatrixUtils.createRealMatrix(U);

    }
}

From source file:dr.math.Procrustes.java

public Procrustes(RealMatrix X, RealMatrix Xstar, boolean allowTranslation, boolean allowDilation) {
    rowDimension = X.getRowDimension();// w  w  w . j  ava  2 s . c  om
    columnDimension = X.getColumnDimension();

    if (Xstar.getRowDimension() != rowDimension) {
        throw new IllegalArgumentException("X and Xstar do not have the same number of rows");
    }
    if (Xstar.getColumnDimension() != columnDimension) {
        throw new IllegalArgumentException("X and Xstar do not have the same number of columns");
    }

    RealMatrix J = new Array2DRowRealMatrix(rowDimension, rowDimension);

    if (allowTranslation) {
        //           J <- diag(n) - 1/n * matrix(1, n, n)
        //           for n = 3, J = {{1, -2/3, -2/3}, {-2/3, 1, -2/3}, {-2/3, -2/3, 1}}

        for (int i = 0; i < rowDimension; i++) {
            J.setEntry(i, i, 1.0 - (1.0 / rowDimension));

            for (int j = i + 1; j < rowDimension; j++) {
                J.setEntry(i, j, -1.0 / rowDimension);
                J.setEntry(j, i, -1.0 / rowDimension);
            }
        }
    } else {
        //           J <- diag(n)

        for (int i = 0; i < rowDimension; i++) {
            J.setEntry(i, i, 1);
        }

    }

    //       C <- t(Xstar) %*% J %*% X

    RealMatrix C = Xstar.transpose().multiply(J.multiply(X));

    //       svd.out <- svd(C)
    //       R <- svd.out$v %*% t(svd.out$u)
    //      NB: Apache math does a different SVD from R.  TODO Should use Colt library
    SingularValueDecomposition SVD = new SingularValueDecompositionImpl(C);
    R = SVD.getV().multiply(SVD.getUT());

    //       s <- 1
    double s = 1.0; // scale = 1 unless dilation is being used

    if (allowDilation) {
        //           mat1 <- t(Xstar) %*% J %*% X %*% R
        RealMatrix mat1 = Xstar.transpose().multiply(J.multiply(X.multiply(R)));

        //           mat2 <- t(X) %*% J %*% X
        RealMatrix mat2 = X.transpose().multiply(J.multiply(X));

        //           s.numer <- 0
        //           s.denom <- 0
        double numer = 0.0;
        double denom = 0.0;

        //           for (i in 1:m) {
        //               s.numer <- s.numer + mat1[i, i]
        //               s.denom <- s.denom + mat2[i, i]
        //           }
        for (int i = 0; i < columnDimension; i++) {
            numer = numer + mat1.getEntry(i, i);
            denom = denom + mat2.getEntry(i, i);
        }
        //           s <- s.numer/s.denom
        s = numer / denom;
    }
    this.s = s;

    //       tt <- matrix(0, m, 1)
    RealMatrix tmpT = new Array2DRowRealMatrix(columnDimension, 1); // a translation vector of zero unless translation is being used

    if (allowTranslation) {
        //           tt <- 1/n * t(Xstar - s * X %*% R) %*% matrix(1, n, 1)
        RealMatrix tmp = new Array2DRowRealMatrix(rowDimension, 1);
        for (int i = 0; i < rowDimension; i++) {
            tmp.setEntry(i, 0, 1);
        }
        tmpT = Xstar.subtract(X.multiply(R).scalarMultiply(s)).transpose().scalarMultiply(1.0 / rowDimension)
                .multiply(tmp);
    }

    T = tmpT;
}

From source file:name.mjw.cytospade.fcsFile.java

/**
 * getCompensatedEventList ---/*from w  w w .j  a v a 2s  .  c om*/
 * <p>
 * Returns the event list compensated by the SPILL matrix.
 * <p>
 *
 * @return array of double arrays containing the events.
 */
public double[][] getCompensatedEventList() {
    double[][] events = this.getEventList();
    if (events.length != this.getNumChannels())
        return events; // Unable to extract the underlying events

    // Convert the SPILL string to a compensation matrix
    String compString = this.getSpillString();
    if (compString == null)
        return events; // No compensation, just return the events

    // Split the compensation string into its values
    //
    // The basic structure for SPILL* is:
    // $SPILLOVER/n,string1,string2,...,f1,f2,f3,f4,.../

    String[] compValues = compString.split(",");
    String[] compNames = null;
    String[] compData = null;
    int compDataStart = 0;

    int n = 0;
    try {
        // Try to parse the number of acquisition parameters
        n = Integer.parseInt(compValues[0]);
        if (n <= 0 || n > this.parameters)
            throw new NumberFormatException();
    } catch (NumberFormatException nfe) {
        CyLogger.getLogger().error("Failed to parse parameter count in spill string", nfe);
        return events;
    }

    compNames = Arrays.copyOfRange(compValues, 1, n + 1);

    // Match names in spill string to columns in parameter lists
    compDataStart = Arrays.asList(this.channelShortname).indexOf(compNames[0]);
    if (compDataStart < 0) {
        CyLogger.getLogger().error("Failed to match channel " + compNames[0] + " to parameter in file");
        return events; // Failure match spill string names to channels
    }
    for (int i = 0; i < n; i++) {
        if (!compNames[i].equals(this.channelShortname[compDataStart + i])) {
            CyLogger.getLogger().error("Spill channel are not continguous parameters in file");
            return events; // Spill string columns not in order
        }
    }

    // Extract actual compensation data
    compData = Arrays.copyOfRange(compValues, n + 1, compValues.length);
    if (compData.length != (n * n))
        return events;

    /**
     * Populate the compensation matrix --- The values are stored in
     * row-major order, i.e., the elements in the first row appear
     * first.
     */
    double[][] matrix = new double[n][n];

    // Loop through the array of compensation values
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            try {
                matrix[i][j] = Double.parseDouble(compData[i * n + j]);
            } catch (NumberFormatException nfe) {
                // Set default value If a NumberFormatException occurred
                matrix[i][j] = 0.0d;
            }
        }
    }

    // Compute the inverse of the compensation data and then apply
    // to data matrix (which is column major). Specifically compute
    // transpose(inverse(<SPILL MATRIX>)) * data
    RealMatrix comp = (new LUDecompositionImpl(new Array2DRowRealMatrix(matrix))).getSolver().getInverse();
    RealMatrix data = new BlockRealMatrix(events);
    data.setSubMatrix( // Update compensated portion of data matrix
            comp.transpose()
                    .multiply(data.getSubMatrix(compDataStart, compDataStart + n - 1, 0,
                            this.getEventCount() - 1))
                    .getData(),

            compDataStart, 0);
    return data.getData();
}

From source file:gephi.spade.panel.fcsFile.java

/**
 * getCompensatedEventList ---/*from   w  ww .j av a 2 s  .  c  o m*/
 * <p>
 * Returns the event list compensated by the SPILL matrix.
 * <p>
 *
 * @return array of double arrays containing the events.
 */
public double[][] getCompensatedEventList() {
    double[][] events = this.getEventList();
    if (events.length != this.getNumChannels())
        return events; // Unable to extract the underlying events

    // Convert the SPILL string to a compensation matrix
    String compString = this.getSpillString();
    if (compString == null)
        return events; // No compensation, just return the events

    // Split the compensation string into its values
    //
    // The basic structure for SPILL* is:
    // $SPILLOVER/n,string1,string2,...,f1,f2,f3,f4,.../

    String[] compValues = compString.split(",");
    String[] compNames = null;
    String[] compData = null;
    int compDataStart = 0;

    int n = 0;
    try {
        // Try to parse the number of acquisition parameters
        n = Integer.parseInt(compValues[0]);
        if (n <= 0 || n > this.parameters)
            throw new NumberFormatException();
    } catch (NumberFormatException nfe) {
        //CyLogger.getLogger().error("Failed to parse parameter count in spill string",nfe);
        return events;
    }

    compNames = Arrays.copyOfRange(compValues, 1, n + 1);

    // Match names in spill string to columns in parameter lists
    compDataStart = Arrays.asList(this.channelShortname).indexOf(compNames[0]);
    if (compDataStart < 0) {
        //CyLogger.getLogger().error("Failed to match channel "+compNames[0]+" to parameter in file");
        return events; // Failure match spill string names to channels
    }
    for (int i = 0; i < n; i++) {
        if (!compNames[i].equals(this.channelShortname[compDataStart + i])) {
            //CyLogger.getLogger().error("Spill channel are not continguous parameters in file");
            return events; // Spill string columns not in order
        }
    }

    // Extract actual compensation data
    compData = Arrays.copyOfRange(compValues, n + 1, compValues.length);
    if (compData.length != (n * n))
        return events;

    /**
     * Populate the compensation matrix --- The values are stored in
     * row-major order, i.e., the elements in the first row appear
     * first.
     */
    double[][] matrix = new double[n][n];

    // Loop through the array of compensation values
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            try {
                matrix[i][j] = Double.parseDouble(compData[i * n + j]);
            } catch (NumberFormatException nfe) {
                // Set default value If a NumberFormatException occurred
                matrix[i][j] = 0.0d;
            }
        }
    }

    // Compute the inverse of the compensation data and then apply
    // to data matrix (which is column major). Specifically compute
    // transpose(inverse(<SPILL MATRIX>)) * data
    RealMatrix comp = (new LUDecompositionImpl(new Array2DRowRealMatrix(matrix))).getSolver().getInverse();
    RealMatrix data = new BlockRealMatrix(events);
    data.setSubMatrix( // Update compensated portion of data matrix
            comp.transpose()
                    .multiply(data.getSubMatrix(compDataStart, compDataStart + n - 1, 0,
                            this.getEventCount() - 1))
                    .getData(),

            compDataStart, 0);
    return data.getData();
}

From source file:org.plista.kornakapi.core.recommender.FoldingFactorization.java

private double[][] computeUserFoldInMatrix(double[][] itemFeatures) {

    /* if there are no items, we cannot fold in anything */
    if (itemFeatures.length == 0) {
        return new double[0][0];
    }/* w  ww.j a va  2s .c  om*/

    if (log.isInfoEnabled()) {
        log.info("Computing fold-in matrix from a {} x {} item features matrix", factorization.numItems(),
                factorization.numFeatures());
    }

    RealMatrix Y = new Array2DRowRealMatrix(itemFeatures);
    RealMatrix YTY = Y.transpose().multiply(Y);
    RealMatrix YTYInverse = new LUDecompositionImpl(YTY).getSolver().getInverse();

    return Y.multiply(YTYInverse).getData();
}