Example usage for org.apache.commons.math3.linear RealMatrix getRowDimension

List of usage examples for org.apache.commons.math3.linear RealMatrix getRowDimension

  1. HOME
  2. Java
  3. org.apache.commons
  4. org.apache.commons.math3.linear.*
  5. RealMatrix
  6. getRowDimension

Introduction

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

Prototype

int getRowDimension();

Source Link

Document

Returns the number of rows in the matrix.

Usage

From source file:com.github.tteofili.looseen.yay.SGM.java

/**
 * predict network output given an input
 *
 * @param input the input/*from   w  w  w .j  ava  2 s  . co  m*/
 * @return the output
 * @throws Exception
 */
private double[] predictOutput(double[] input) throws Exception {

    RealMatrix hidden = rectifierFunction.applyMatrix(
            MatrixUtils.createRowRealMatrix(input).multiply(weights[0].transpose()).add(biases[0]));
    RealMatrix scores = hidden.multiply(weights[1].transpose()).add(biases[1]);

    RealMatrix probs = scores.copy();
    int len = scores.getColumnDimension() - 1;
    for (int d = 0; d < configuration.window - 1; d++) {
        int startColumn = d * len / (configuration.window - 1);
        RealMatrix subMatrix = scores.getSubMatrix(0, scores.getRowDimension() - 1, startColumn,
                startColumn + input.length);
        for (int sm = 0; sm < subMatrix.getRowDimension(); sm++) {
            probs.setSubMatrix(softmaxActivationFunction.applyMatrix(subMatrix.getRowMatrix(sm)).getData(), sm,
                    startColumn);
        }
    }

    RealVector d = probs.getRowVector(0);
    return d.toArray();
}

From source file:edu.stanford.cfuller.imageanalysistools.filter.LocalBackgroundEstimationFilter.java

/**
 * Applies the LocalBackgroundEstimationFilter to an Image.
 * @param im    The Image that will be replaced by the output Image.  This can be anything of the correct dimensions except a shallow copy of the reference Image.
 *///from w ww .  j  ava  2s.  c om
@Override
public void apply(WritableImage im) {

    if (this.referenceImage == null) {
        throw new ReferenceImageRequiredException(
                "LocalBackgroundEstimationFilter requires a reference image.");
    }

    edu.stanford.cfuller.imageanalysistools.image.Histogram h = new edu.stanford.cfuller.imageanalysistools.image.Histogram(
            this.referenceImage);

    int numPixelsInBox = boxSize * boxSize;

    ImageCoordinate ic = this.referenceImage.getDimensionSizes();

    ImageCoordinate icnew = ImageCoordinate.createCoordXYZCT(ic.get(ImageCoordinate.X) + 2 * boxSize,
            ic.get(ImageCoordinate.Y) + 2 * boxSize, ic.get(ImageCoordinate.Z), ic.get(ImageCoordinate.C),
            ic.get(ImageCoordinate.T));

    WritableImage padded = ImageFactory.createWritable(icnew, -1.0f);

    float maxValue = 0;

    for (ImageCoordinate i : this.referenceImage) {

        icnew.quickSet(ImageCoordinate.X, i.quickGet(ImageCoordinate.X) + boxSize);
        icnew.quickSet(ImageCoordinate.Y, i.quickGet(ImageCoordinate.Y) + boxSize);
        icnew.quickSet(ImageCoordinate.Z, i.quickGet(ImageCoordinate.Z));
        icnew.quickSet(ImageCoordinate.C, i.quickGet(ImageCoordinate.C));
        icnew.quickSet(ImageCoordinate.T, i.quickGet(ImageCoordinate.T));

        padded.setValue(icnew, this.referenceImage.getValue(i));

        if (this.referenceImage.getValue(i) > maxValue)
            maxValue = this.referenceImage.getValue(i);

    }

    RealVector overallCounts = new org.apache.commons.math3.linear.ArrayRealVector(h.getMaxValue() + 1);

    RealMatrix countsByRow = new org.apache.commons.math3.linear.Array2DRowRealMatrix(2 * boxSize + 1,
            h.getMaxValue() + 1);

    //loop over columns

    for (int i = boxSize; i < im.getDimensionSizes().quickGet(ImageCoordinate.X) + boxSize; i++) {

        overallCounts.mapMultiplyToSelf(0.0);
        double[] overallCounts_a = overallCounts.toArray();
        countsByRow = countsByRow.scalarMultiply(0.0);
        double[][] countsByRow_a = countsByRow.getData();

        int countsByRow_rowZero_pointer = 0;

        for (int m = i - boxSize; m <= i + boxSize; m++) {
            for (int n = 0; n < 2 * boxSize + 1; n++) {
                icnew.quickSet(ImageCoordinate.X, m);
                icnew.quickSet(ImageCoordinate.Y, n);
                int value = (int) padded.getValue(icnew);

                if (value == -1)
                    continue;

                overallCounts_a[value]++;
                countsByRow_a[(n + countsByRow_rowZero_pointer) % countsByRow.getRowDimension()][value]++;

            }
        }

        int currMedian = 0;
        int runningSum = 0;
        int k = 0;

        while (runningSum < (numPixelsInBox >> 1)) {
            runningSum += overallCounts_a[k++];
        }

        runningSum -= overallCounts_a[k - 1];

        currMedian = k - 1;

        icnew.quickSet(ImageCoordinate.X, i - boxSize);
        icnew.quickSet(ImageCoordinate.Y, 0);

        im.setValue(icnew, currMedian);

        int num_rows = countsByRow.getRowDimension();

        for (int j = boxSize + 1; j < im.getDimensionSizes().quickGet(ImageCoordinate.Y) + boxSize; j++) {

            double[] toRemove = countsByRow_a[(countsByRow_rowZero_pointer) % num_rows];

            for (int oc_counter = 0; oc_counter < overallCounts_a.length; oc_counter++) {
                overallCounts_a[oc_counter] -= toRemove[oc_counter];
            }

            for (int c = 0; c < toRemove.length; c++) {
                if (c < currMedian) {
                    runningSum -= toRemove[c];
                }

                countsByRow_a[countsByRow_rowZero_pointer % num_rows][c] *= 0.0;
            }

            countsByRow_rowZero_pointer++;

            for (int c = i - boxSize; c <= i + boxSize; c++) {
                icnew.quickSet(ImageCoordinate.X, c);
                icnew.quickSet(ImageCoordinate.Y, j + boxSize);
                int value = (int) padded.getValue(icnew);

                if (value == -1)
                    continue;

                countsByRow_a[(countsByRow_rowZero_pointer + num_rows - 1) % num_rows][value] += 1;

                overallCounts_a[value]++;

                if (value < currMedian) {
                    runningSum++;
                }

            }

            //case 1: runningSum > half of box

            if (runningSum > (numPixelsInBox >> 1)) {

                k = currMedian - 1;
                while (runningSum > (numPixelsInBox >> 1)) {

                    runningSum -= overallCounts_a[k--];
                }

                currMedian = k + 1;

            } else if (runningSum < (numPixelsInBox >> 1)) { // case 2: runningSum < half of box

                k = currMedian;

                while (runningSum < (numPixelsInBox >> 1)) {

                    runningSum += overallCounts_a[k++];
                }

                currMedian = k - 1;

                runningSum -= overallCounts_a[k - 1];

            }

            //cast 3: spot on, do nothing

            icnew.quickSet(ImageCoordinate.X, i - boxSize);
            icnew.quickSet(ImageCoordinate.Y, j - boxSize);

            im.setValue(icnew, currMedian);

        }

    }

    icnew.recycle();

}

From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ2.java

private boolean indepCollection2(List<Node> aa, List<Node> bb, List<Node> cc, int n, double alpha) {
    //        List<Double> ret = getCutoff1(aa, bb, cc);
    //        List<Double> ret = getCutoff2(aa, bb, cc);
    //        double mean = ret.get(0);
    //        double sd = ret.get(1);

    int numPerm = 10;

    TetradMatrix submatrix = subMatrix(cov, aa, bb, cc);

    TetradMatrix inverse;// w  ww . j av  a 2 s .  co  m
    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();
    }

    List<Double> pValues = new ArrayList<Double>();

    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 = n - 1 - rank;

            if (dof < 0) {
                System.out.println("Negative dof: " + dof + " n = " + n + " 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)));

            pValues.add(p);
        }
    }

    List<Double> zeroes = new ArrayList<Double>();

    for (double d : pValues) {
        if (d == 0)
            zeroes.add(d);
    }

    int k = 0;

    for (double p : pValues) {
        if (p < alpha) {
            k++;
        }
    }

    numTests++;

    int count = countGreater(aa, bb, cc, numPerm, k);
    //        int count = countGreater2(aa, bb, cc, numPerm, k);

    double p = count / (double) numPerm;

    boolean indep = p > alpha;

    //        double p = (1.0 - RandomUtil.getInstance().normalCdf(0, 1, (k - mean) / sd));

    //        boolean indep = p > alpha;

    if (verbose) {
        //            System.out.println("\n#accepted " + k + " cutoff = " + kCutoff + " indep = " + indep);
        //            System.out.println("\n#accepted " + k + " meam = " + mean + " sd = " + sd + " indep = " + indep);
        //            System.out.println("standard deviations " + (k - mean) / sd + " p = " + p + " alpha = " + alpha);
        System.out.println("\n#accepted " + k + " indep = " + indep);
        System.out.println("p = " + p + " alpha = " + alpha);
    }

    numTests++;
    if (!indep)
        numDependent++;

    return indep;

}

From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ2.java

private boolean indepCollection(List<Node> aa, List<Node> bb, List<Node> cc, int n, double alpha) {
    //        List<Double> ret = getCutoff1(aa, bb, cc);
    //        List<Double> ret = getCutoff2(aa, bb, cc);
    //        double mean = ret.get(0);
    //        double sd = ret.get(1);

    int numPerm = 10;

    TetradMatrix submatrix = subMatrix(cov, aa, bb, cc);

    TetradMatrix inverse;// w  ww .j  a va 2  s . c om
    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();
    }

    List<Double> pValues = new ArrayList<Double>();

    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 = n - 1 - rank;

            if (dof < 0) {
                System.out.println("Negative dof: " + dof + " n = " + n + " 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)));

            pValues.add(p);
        }
    }

    List<Double> zeroes = new ArrayList<Double>();

    for (double d : pValues) {
        if (d == 0)
            zeroes.add(d);
    }

    int k = 0;

    for (double p : pValues) {
        if (p < alpha) {
            k++;
        }
    }

    numTests++;

    int count = countGreater(aa, bb, cc, numPerm, k);
    //        int count = countGreater2(aa, bb, cc, numPerm, k);
    //        int count = countGreater3(aa, bb, cc, numPerm, k, submatrix, inverse);
    //        int count = countGreater4(aa, bb, cc, numPerm, k);

    double p = count / (double) numPerm;

    boolean indep = p > alpha;

    //        double p = (1.0 - RandomUtil.getInstance().normalCdf(0, 1, (k - mean) / sd));

    //        boolean indep = p > alpha;

    if (verbose) {
        //            System.out.println("\n#accepted " + k + " cutoff = " + kCutoff + " indep = " + indep);
        //            System.out.println("\n#accepted " + k + " meam = " + mean + " sd = " + sd + " indep = " + indep);
        //            System.out.println("standard deviations " + (k - mean) / sd + " p = " + p + " alpha = " + alpha);
        System.out.println("\n#accepted " + k + " indep = " + indep);
        System.out.println("p = " + p + " alpha = " + alpha);
    }

    numTests++;
    if (!indep)
        numDependent++;

    return indep;

}

From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ2.java

public boolean isIndependent(Node x, Node y, List<Node> z) {
    if (verbose) {
        System.out.println("\n" + x + " _||_ " + y + " | " + z);
    }//from w  ww  .j  a  va  2s  . c om

    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 = 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>();

    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 - inverse.columns();

            if (dof < 0) {
                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_)));

            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));
        }
    });

    // Sort pvalues, _i, and _m together.
    List<Double> pValues2 = new ArrayList<Double>();
    List<Integer> _i2 = new ArrayList<Integer>();
    List<Integer> _m2 = new ArrayList<Integer>();

    for (int _y = 0; _y < indices.size(); _y++) {
        pValues2.add(pValues.get(indices.get(_y)));
        _i2.add(_i.get(indices.get(_y)));
        _m2.add(_m.get(indices.get(_y)));
    }

    Collections.sort(pValues2);

    int k = StatUtils.fdr(alpha, pValues, false);

    int nonzero = -1;

    for (int i = 0; i < pValues2.size(); i++) {
        if (pValues.get(i) != 0) {
            nonzero = i;
            break;
        }
    }

    if (nonzero < k) {
        for (int g = 0; g < k; g++) {
            int x3 = _i2.get(g);
            int y3 = _m2.get(g);
            twod[x3][y3] = 1;
        }

        //            if (verbose) {
        //                System.out.println("Dependent");
        //            }

        return false;
    } else {
        if (verbose) {
            System.out.println("Independent");
        }

        return true;
    }
}

From source file:com.joptimizer.algebra.CholeskyRCFactorizationTest.java

public void testInvert1() throws Exception {
    log.debug("testInvert1");
    double[][] QData = new double[][] { { 1, .12, .13, .14, .15 }, { .12, 2, .23, .24, .25 },
            { .13, .23, 3, 0, 0 }, { .14, .24, 0, 4, 0 }, { .15, .25, 0, 0, 5 } };
    RealMatrix Q = MatrixUtils.createRealMatrix(QData);

    CholeskyRCFactorization myc = new CholeskyRCFactorization(DoubleFactory2D.dense.make(QData));
    myc.factorize();/*from   w  w  w .  j a  v a 2  s. c o m*/
    RealMatrix L = new Array2DRowRealMatrix(myc.getL().toArray());
    RealMatrix LT = new Array2DRowRealMatrix(myc.getLT().toArray());
    log.debug("L: " + ArrayUtils.toString(L.getData()));
    log.debug("LT: " + ArrayUtils.toString(LT.getData()));
    log.debug("L.LT: " + ArrayUtils.toString(L.multiply(LT).getData()));
    log.debug("LT.L: " + ArrayUtils.toString(LT.multiply(L).getData()));

    // check Q = L.LT
    double norm = L.multiply(LT).subtract(Q).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);

    RealMatrix LInv = new SingularValueDecomposition(L).getSolver().getInverse();
    log.debug("LInv: " + ArrayUtils.toString(LInv.getData()));
    RealMatrix LInvT = LInv.transpose();
    log.debug("LInvT: " + ArrayUtils.toString(LInvT.getData()));
    RealMatrix LTInv = new SingularValueDecomposition(LT).getSolver().getInverse();
    log.debug("LTInv: " + ArrayUtils.toString(LTInv.getData()));
    RealMatrix LTInvT = LTInv.transpose();
    log.debug("LTInvT: " + ArrayUtils.toString(LTInvT.getData()));
    log.debug("LInv.LInvT: " + ArrayUtils.toString(LInv.multiply(LInvT).getData()));
    log.debug("LTInv.LTInvT: " + ArrayUtils.toString(LTInv.multiply(LTInvT).getData()));

    RealMatrix Id = MatrixUtils.createRealIdentityMatrix(Q.getRowDimension());
    //check Q.(LTInv * LInv) = 1
    norm = Q.multiply(LTInv.multiply(LInv)).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 5.E-15);

    // check Q.QInv = 1
    RealMatrix QInv = MatrixUtils.createRealMatrix(myc.getInverse().toArray());
    norm = Q.multiply(QInv).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);
}

From source file:com.joptimizer.algebra.CholeskyRCTFactorizationTest.java

public void testInvert1() throws Exception {
    log.debug("testInvert1");
    double[][] QData = new double[][] { { 1, .12, .13, .14, .15 }, { .12, 2, .23, .24, .25 },
            { .13, .23, 3, 0, 0 }, { .14, .24, 0, 4, 0 }, { .15, .25, 0, 0, 5 } };
    RealMatrix Q = MatrixUtils.createRealMatrix(QData);

    CholeskyRCTFactorization myc = new CholeskyRCTFactorization(DoubleFactory2D.dense.make(QData));
    myc.factorize();/*from ww  w. j  a v a 2s. co m*/
    RealMatrix L = new Array2DRowRealMatrix(myc.getL().toArray());
    RealMatrix LT = new Array2DRowRealMatrix(myc.getLT().toArray());
    log.debug("L: " + ArrayUtils.toString(L.getData()));
    log.debug("LT: " + ArrayUtils.toString(LT.getData()));
    log.debug("L.LT: " + ArrayUtils.toString(L.multiply(LT).getData()));
    log.debug("LT.L: " + ArrayUtils.toString(LT.multiply(L).getData()));

    // check Q = L.LT
    double norm = L.multiply(LT).subtract(Q).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);

    RealMatrix LInv = new SingularValueDecomposition(L).getSolver().getInverse();
    log.debug("LInv: " + ArrayUtils.toString(LInv.getData()));
    RealMatrix LInvT = LInv.transpose();
    log.debug("LInvT: " + ArrayUtils.toString(LInvT.getData()));
    RealMatrix LTInv = new SingularValueDecomposition(LT).getSolver().getInverse();
    log.debug("LTInv: " + ArrayUtils.toString(LTInv.getData()));
    RealMatrix LTInvT = LTInv.transpose();
    log.debug("LTInvT: " + ArrayUtils.toString(LTInvT.getData()));
    log.debug("LInv.LInvT: " + ArrayUtils.toString(LInv.multiply(LInvT).getData()));
    log.debug("LTInv.LTInvT: " + ArrayUtils.toString(LTInv.multiply(LTInvT).getData()));

    RealMatrix Id = MatrixUtils.createRealIdentityMatrix(Q.getRowDimension());
    //check Q.(LTInv * LInv) = 1
    norm = Q.multiply(LTInv.multiply(LInv)).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 5.E-15);

    // check Q.QInv = 1
    RealMatrix QInv = MatrixUtils.createRealMatrix(myc.getInverse().toArray());
    norm = Q.multiply(QInv).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);
}

From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java

/**
 * Pre-calculate intermediate matrices, perform setup, etc.
 *//*from  ww w . ja  v  a 2 s.c om*/
private void initialize(Test test, RealMatrix FEssence, RealMatrix FtFinverse, int perGroupN,
        FixedRandomMatrix CFixedRand, RealMatrix U, RealMatrix thetaNull, RealMatrix beta,
        RealMatrix sigmaError, RealMatrix sigmaG, boolean exact) throws PowerException {
    debug("entering initialize");

    // reset member variables
    this.T1 = null;
    this.FT1 = null;
    this.S = null;
    this.mzSq = null;
    this.H0 = 0;
    this.sStar = 0;

    // cache inputs
    this.test = test;
    this.FEssence = FEssence;
    this.FtFinverse = FtFinverse;
    this.perGroupN = perGroupN;
    this.CFixedRand = CFixedRand;
    this.U = U;
    this.thetaNull = thetaNull;
    this.beta = beta;
    this.sigmaError = sigmaError;
    this.sigmaG = sigmaG;

    // calculate intermediate matrices
    //        RealMatrix FEssence = params.getDesignEssence().getFullDesignMatrixFixed();
    // TODO: do we ever get here with values that can cause integer overflow,
    //       and if so, does it matter?
    this.N = (double) FEssence.getRowDimension() * perGroupN;
    this.exact = exact;
    try {
        // TODO: need to calculate H0, need to adjust H1 for Unirep
        // get design matrix for fixed parameters only
        qF = FEssence.getColumnDimension();
        // a = CFixedRand.getCombinedMatrix().getRowDimension();

        // get fixed contrasts
        RealMatrix Cfixed = CFixedRand.getFixedMatrix();
        RealMatrix CGaussian = CFixedRand.getRandomMatrix();

        // build intermediate terms h1, S
        if (FtFinverse == null) {
            FtFinverse = new LUDecomposition(FEssence.transpose().multiply(FEssence)).getSolver().getInverse();
            debug("FEssence", FEssence);
            debug("FtFinverse = (FEssence transpose * FEssence) inverse", FtFinverse);
        } else {
            debug("FtFinverse", FtFinverse);
        }

        RealMatrix PPt = Cfixed.multiply(FtFinverse.scalarMultiply(1 / (double) perGroupN))
                .multiply(Cfixed.transpose());
        debug("Cfixed", Cfixed);
        debug("n = " + perGroupN);
        debug("PPt = Cfixed * FtF inverse * (1/n) * Cfixed transpose", PPt);

        T1 = forceSymmetric(new LUDecomposition(PPt).getSolver().getInverse());
        debug("T1 = PPt inverse", T1);

        FT1 = new CholeskyDecomposition(T1).getL();
        debug("FT1 = Cholesky decomposition (L) of T1", FT1);

        // calculate theta difference
        //            RealMatrix thetaNull = params.getTheta();
        RealMatrix C = CFixedRand.getCombinedMatrix();
        //            RealMatrix beta = params.getScaledBeta();
        //            RealMatrix U = params.getWithinSubjectContrast();

        // thetaHat = C * beta * U
        RealMatrix thetaHat = C.multiply(beta.multiply(U));
        debug("C", C);
        debug("beta", beta);
        debug("U", U);
        debug("thetaHat = C * beta * U", thetaHat);

        // thetaDiff = thetaHat - thetaNull
        RealMatrix thetaDiff = thetaHat.subtract(thetaNull);
        debug("thetaNull", thetaNull);
        debug("thetaDiff = thetaHat - thetaNull", thetaDiff);

        // TODO: specific to HLT or UNIREP
        RealMatrix sigmaStarInverse = getSigmaStarInverse(U, sigmaError, test);
        debug("sigmaStarInverse", sigmaStarInverse);

        RealMatrix H1matrix = thetaDiff.transpose().multiply(T1).multiply(thetaDiff).multiply(sigmaStarInverse);
        debug("H1matrix = thetaDiff transpose * T1 * thetaDiff * sigmaStarInverse", H1matrix);

        H1 = H1matrix.getTrace();
        debug("H1 = " + H1);

        // Matrix which represents the non-centrality parameter as a linear combination of chi-squared r.v.'s.
        S = FT1.transpose().multiply(thetaDiff).multiply(sigmaStarInverse).multiply(thetaDiff.transpose())
                .multiply(FT1).scalarMultiply(1 / H1);
        debug("S = FT1 transpose * thetaDiff * sigmaStar inverse * thetaDiff transpose * FT1 * (1/H1)", S);

        // We use the S matrix to generate the F-critical, numerical df's, and denominator df's
        // for a central F distribution.  The resulting F distribution is used as an approximation
        // for the distribution of the non-centrality parameter.
        // See formulas 18-21 and A8,A10 from Glueck & Muller (2003) for details.
        EigenDecomposition sEigenDecomp = new EigenDecomposition(S);
        sEigenValues = sEigenDecomp.getRealEigenvalues();
        // calculate H0
        if (sEigenValues.length > 0)
            H0 = H1 * (1 - sEigenValues[0]);
        if (H0 <= 0)
            H0 = 0;

        // count the # of positive eigen values
        for (double value : sEigenValues) {
            if (value > 0)
                sStar++;
        }
        // TODO: throw error if sStar is <= 0
        // TODO: NO: throw error if sStar != sEigenValues.length instead???
        double stddevG = Math.sqrt(sigmaG.getEntry(0, 0));
        RealMatrix svec = sEigenDecomp.getVT();
        mzSq = svec.multiply(FT1.transpose()).multiply(CGaussian).scalarMultiply(1 / stddevG);
        for (int i = 0; i < mzSq.getRowDimension(); i++) {
            for (int j = 0; j < mzSq.getColumnDimension(); j++) {
                double entry = mzSq.getEntry(i, j);
                mzSq.setEntry(i, j, entry * entry); // TODO: is there an apache function to do this?
            }
        }

        debug("exiting initialize normally");
    } catch (RuntimeException e) {
        LOGGER.warn("exiting initialize abnormally", e);

        throw new PowerException(e.getMessage(), PowerErrorEnum.INVALID_DISTRIBUTION_NONCENTRALITY_PARAMETER);
    }
}

From source file:edu.dfci.cccb.mev.domain.Heatmap.java

private RealMatrix transpose(final RealMatrix original) {
    return new AbstractRealMatrix() {

        @Override/*from  ww  w  .j  a  va2s . c  o  m*/
        public void setEntry(int row, int column, double value) throws OutOfRangeException {
            original.setEntry(column, row, value);
        }

        @Override
        public int getRowDimension() {
            return original.getColumnDimension();
        }

        @Override
        public double getEntry(int row, int column) throws OutOfRangeException {
            return original.getEntry(column, row);
        }

        @Override
        public int getColumnDimension() {
            return original.getRowDimension();
        }

        @Override
        public RealMatrix createMatrix(int rowDimension, int columnDimension)
                throws NotStrictlyPositiveException {
            return original.createMatrix(rowDimension, columnDimension);
        }

        @Override
        public RealMatrix copy() {
            final RealMatrix result = createMatrix(getRowDimension(), getColumnDimension());
            walkInOptimizedOrder(new RealMatrixPreservingVisitor() {

                @Override
                public void visit(int row, int column, double value) {
                    result.setEntry(row, column, value);
                }

                @Override
                public void start(int rows, int columns, int startRow, int endRow, int startColumn,
                        int endColumn) {
                }

                @Override
                public double end() {
                    return NaN;
                }
            });
            return result;
        }
    };
}

From source file:edu.cudenver.bios.power.GLMMPowerCalculator.java

/**
 * Perform any preliminary calculations / updates on the input matrices
 * @param params input parameters//from  w  ww  .  j a  v a2s. com
 */
private void initialize(GLMMPowerParameters params) {
    debug("entering initialize");

    // TODO: why isn't this done in PowerResourceHelper.studyDesignToPowerParameters?
    // if no power methods are specified, add conditional as a default
    if (params.getPowerMethodList().size() <= 0)
        params.addPowerMethod(PowerMethod.CONDITIONAL_POWER);

    // update sigma error and beta if we have a baseline covariate
    RealMatrix sigmaG = params.getSigmaGaussianRandom();
    debug("sigmaG", sigmaG);

    int numRandom = sigmaG != null ? sigmaG.getRowDimension() : 0;
    if (numRandom == 1) {
        // set the sigma error matrix to [sigmaY - sigmaYG * sigmaG-1 * sigmaGY]
        RealMatrix sigmaY = params.getSigmaOutcome();
        debug("sigmaY", sigmaY);

        RealMatrix sigmaYG = params.getSigmaOutcomeGaussianRandom();
        debug("sigmaYG", sigmaYG);

        RealMatrix sigmaGY = sigmaYG.transpose();
        debug("sigmaYG transpose", sigmaGY);

        RealMatrix sigmaGInverse = new LUDecomposition(sigmaG).getSolver().getInverse();
        debug("sigmaG inverse", sigmaGInverse);

        RealMatrix sigmaError = forceSymmetric(
                sigmaY.subtract(sigmaYG.multiply(sigmaGInverse.multiply(sigmaGY))));
        debug("sigmaError = sigmaY - sigmaYG * sigmaG inverse * sigmaYG transpose", sigmaError);

        if (!MatrixUtils.isPositiveSemidefinite(sigmaError)) {
            throw new IllegalArgumentException(SIGMA_ERROR_NOT_POSITIVE_SEMIDEFINITE_MESSAGE);
        }
        params.setSigmaError(sigmaError);

        // calculate the betaG matrix and fill in the placeholder row for the random predictor
        FixedRandomMatrix beta = params.getBeta();
        beta.updateRandomMatrix(sigmaGInverse.multiply(sigmaGY));
        debug("beta with random part updated to sigmaG inverse * sigmaYG transpose", beta.getCombinedMatrix());
    }

    debug("exiting initialize");
}