Java Code Examples for org.ejml.data.DenseMatrix64F#get()

public PermutationIndices(DenseMatrix64F matrix) {

        this.matrix = matrix;
        dim = matrix.getNumCols();
        assert (dim == matrix.getNumRows());

        for (int i = 0; i < dim; ++i) {
            double diagonal = matrix.get(i, i);
            if (Double.isInfinite(diagonal)) {
                ++infiniteCount;
            } else if (diagonal == 0.0) {
                ++zeroCount;
            } else {
                ++nonZeroFiniteCount;
            }
        }
    }
 

public static double det(DenseMatrix64F mat) {
    int numCol = mat.getNumCols();
    int numRow = mat.getNumRows();
    if (numCol != numRow) {
        throw new IllegalArgumentException("Must be a square matrix.");
    } else if (numCol <= 6) {
        return numCol >= 2 ? UnrolledDeterminantFromMinor.det(mat) : mat.get(0);
    } else {
        LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64();
        if (alg.inputModified()) {
            mat = mat.copy();
        }

        return !alg.decompose(mat) ? 0.0D : alg.computeDeterminant().real;
    }
}
 

public static double invertAndGetDeterminant(DenseMatrix64F mat, DenseMatrix64F result, boolean log) {

        final int numCol = mat.getNumCols();
        final int numRow = mat.getNumRows();
        if (numCol != numRow) {
            throw new IllegalArgumentException("Must be a square matrix.");
        }

        if (numCol <= 5) {

            if (numCol >= 2) {
                UnrolledInverseFromMinor.inv(mat, result);
            } else {
                result.set(0, 1.0D / mat.get(0));
            }

            double det = numCol >= 2 ?
                    UnrolledDeterminantFromMinor.det(mat) :
                    mat.get(0);
            return log ? Math.log(det) : det;

        } else {

            LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64();
            LinearSolverLu_D64 solver = new LinearSolverLu_D64(alg);
            if (solver.modifiesA()) {
                mat = mat.copy();
            }

            if (!solver.setA(mat)) {
                return Double.NaN;
            }

            solver.invert(result);

            return log ? computeLogDeterminant(alg) : alg.computeDeterminant().real;

        }
    }
 

@Override
void setMatrixRatio(DenseMatrix64F numeratorMatrix, DenseMatrix64F otherMatrix,
                    DenseMatrix64F destination) {

    for (int i = 0; i < destination.numRows; i++) {

        double val = numeratorMatrix.get(i, i) / (numeratorMatrix.get(i, i) + otherMatrix.get(i, i));
        destination.set(i, i, val);
        for (int j = i + 1; j < destination.numRows; j++) {

            double rg = numeratorMatrix.get(i, j);
            double vi = numeratorMatrix.get(i, i) + otherMatrix.get(i, i);
            double vj = numeratorMatrix.get(j, j) + otherMatrix.get(j, j);

            val = rg / sqrt(vi * vj);

            destination.set(i, j, val);
            destination.set(j, i, val);

        }
    }
}
 

public MatrixParameter getEigenVectorsStrengthOfSelection() {
    int dim = getDimTrait();
    DenseMatrix64F V = EigenOps.createMatrixV(eigenDecompositionStrengthOfSelection);

    Parameter[] eigenVectors = new Parameter[getDimTrait()];
    double[] column = new double[dim - 1];
    double norm;
    double sign;
    for (int i = 0; i < dim; i++) {
        norm = 0.0;
        sign = 1.0;
        for (int j = 0; j < dim; j++) {
            norm += V.get(j, i) * V.get(j, i);
        }
        norm = Math.sqrt(norm);
        if (V.get(dim - 1, i) < 0) sign = -1.0;
        for (int j = 0; j < dim - 1; j++) {
            column[j] = V.get(j, i) / norm * sign;
        }
        eigenVectors[i] = new Parameter.Default(column);
    }
    return new MatrixParameter("attenuationMatrix", eigenVectors);
}
 

private double[][] getXtXInverse(MatrixFormulation matrixFormulation) {
    DenseMatrix64F XtXInvMatrix = matrixFormulation.XtXInv.copy();
    int dim = XtXInvMatrix.getNumCols();
    double[][] XtXInvArray = new double[dim][dim];
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            XtXInvArray[i][j] = XtXInvMatrix.get(i, j);
        }
    }
    return XtXInvArray;
}
 

public static void blockUnwrap(final DenseMatrix64F block, final double[] destination,
                               final int destinationOffset,
                               final int offsetRow, final int offsetCol,
                               final int nCols) {
    for (int i = 0; i < block.getNumRows(); i++) { // Rows
        for (int j = 0; j < block.getNumCols(); j++) {
            destination[destinationOffset + (i + offsetRow) * nCols + j + offsetCol] = block.get(i, j);
        }
    }
}
 

private static void rotateNd(double[] x, int dim) {

        // Get first `dim` locations
        DenseMatrix64F matrix = new DenseMatrix64F(dim, dim);
        for (int row = 0; row < dim; ++row) {
            for (int col = 0; col < dim; ++col) {
                matrix.set(row, col, x[col * dim + row]);
            }
        }

        // Do a QR decomposition
        QRDecomposition<DenseMatrix64F> qr = DecompositionFactory.qr(dim, dim);
        qr.decompose(matrix);
        DenseMatrix64F qm = qr.getQ(null, true);
        DenseMatrix64F rm = qr.getR(null, true);

        // Reflection invariance
        if (rm.get(0,0) < 0) {
            CommonOps.scale(-1, rm);
            CommonOps.scale(-1, qm);
        }

        // Compute Q^{-1}
        DenseMatrix64F qInv = new DenseMatrix64F(dim, dim);
        CommonOps.transpose(qm, qInv);

        // Apply to each location
        for (int location = 0; location < x.length / dim; ++location) {
            WrappedVector locationVector = new WrappedVector.Raw(x, location * dim, dim);
            MissingOps.matrixVectorMultiple(qInv, locationVector, locationVector, dim);
        }
    }
 

@Override
public double[] chainRule(BranchSufficientStatistics statistics, NodeRef node,
                          DenseMatrix64F gradQInv, DenseMatrix64F gradN) {

    final double rate = branchRateModel.getBranchRate(tree, node);
    final double differential = branchRateModel.getBranchRateDifferential(tree, node);
    final double scaling = differential / rate;

    // Q_i w.r.t. rate
    DenseMatrix64F gradMatQInv = matrixJacobianQInv;
    CommonOps.scale(scaling, statistics.getBranch().getRawVariance(), gradMatQInv);

    double[] gradient = new double[1];
    for (int i = 0; i < gradMatQInv.getNumElements(); i++) {
        gradient[0] += gradMatQInv.get(i) * gradQInv.get(i);
    }

    // n_i w.r.t. rate
    // TODO: Fix delegate to (possibly) un-link drift from arbitrary rate
    DenseMatrix64F gradMatN = matrixJacobianN;
    CommonOps.scale(scaling, statistics.getBranch().getRawDisplacement(), gradMatN);
    for (int i = 0; i < gradMatN.numRows; i++) {
        gradient[0] += gradMatN.get(i) * gradN.get(i);
    }

    return gradient;

}
 

private void blockUnwrap(final DenseMatrix64F YY, final DenseMatrix64F XX, final DenseMatrix64F XY, final DenseMatrix64F YX, final double[] destination, final int offset) {
    for (int i = 0; i < dimProcess; i++) { // Rows
        for (int j = 0; j < dimProcess; j++) {
            destination[offset + i * dimTrait + j] = YY.get(i, j);
            destination[offset + (i + dimProcess) * dimTrait + j + dimProcess] = XX.get(i, j);
        }
        for (int j = 0; j < dimProcess; j++) {
            destination[offset + i * dimTrait + j + dimProcess] = YX.get(i, j);
            destination[offset + (i + dimProcess) * dimTrait + j] = XY.get(i, j);
        }
    }
}
 

private double sampleNode(NodeRef node, WrappedNormalSufficientStatistics statistics) {

        final int thisNumber = node.getNumber();
        final int[] obsInds = maskDelegate.getObservedIndices(node);
        final int obsDim = obsInds.length;

        if (obsDim == dim) {
            return 0;
        }

        ReadableVector mean = statistics.getMean();
        ReadableMatrix thisP = statistics.getPrecision();
        double precisionScalar = statistics.getPrecisionScalar();
        for (int i = 0; i < mean.getDim(); ++i) {
            fcdMean[i] = mean.get(i);
        }

        for (int i = 0; i < mean.getDim(); ++i) {
            for (int j = 0; j < mean.getDim(); ++j) {
                fcdPrecision[i][j] = thisP.get(i, j) * precisionScalar;
            }
        }

        if (repeatedMeasuresModel != null) {
            //TODO: preallocate memory for all these matrices/vectors
            DenseMatrix64F Q = new DenseMatrix64F(fcdPrecision); //storing original fcd precision
            double[] tipPartial = repeatedMeasuresModel.getTipPartial(thisNumber, false);
            int offset = dim;
            DenseMatrix64F P = MissingOps.wrap(tipPartial, offset, dim, dim);

            DenseMatrix64F preOrderMean = new DenseMatrix64F(dim, 1);
            for (int i = 0; i < dim; i++) {
                preOrderMean.set(i, 0, fcdMean[i]);
            }
            DenseMatrix64F dataMean = new DenseMatrix64F(dim, 1);
            for (int i = 0; i < dim; i++) {
                dataMean.set(i, 0, tipPartial[i]);
            }

            D1Matrix64F bufferMean = new DenseMatrix64F(dim, 1);

            mult(Q, preOrderMean, bufferMean); //bufferMean = Q * preOrderMean
            multAdd(P, dataMean, bufferMean); //bufferMean = Q * preOderMean + P * dataMean


            CommonOps.addEquals(P, Q); //P = P + Q
            DenseMatrix64F V = new DenseMatrix64F(dim, dim);
            CommonOps.invert(P, V); //V = inv(P + Q)
            mult(V, bufferMean, dataMean); // dataMean = inv(P + Q) * (Q * preOderMean + P * dataMean)

            for (int i = 0; i < dim; i++) {
                fcdMean[i] = dataMean.get(i);
                for (int j = 0; j < dim; j++) {
                    fcdPrecision[i][j] = P.get(i, j);
                }

            }
        }

        MultivariateNormalDistribution fullDistribution = new MultivariateNormalDistribution(fcdMean, fcdPrecision); //TODO: should this not be declared until 'else' statement?
        MultivariateNormalDistribution drawDistribution;

        if (mask != null && obsDim > 0) {
            addMaskOnContiuousTraitsPrecisionSpace(thisNumber);
            drawDistribution = new MultivariateNormalDistribution(maskedMean, maskedPrecision);
        } else {
            drawDistribution = fullDistribution;
        }

        double[] oldValue = getNodeTrait(node);

        int attempt = 0;
        boolean validTip = false;

        while (!validTip & attempt < maxAttempts) {

            setNodeTrait(node, drawDistribution.nextMultivariateNormal());

            if (latentLiability.validTraitForTip(thisNumber)) {
                validTip = true;
            }
            attempt++;
        }

        if (attempt == maxAttempts) {
            return Double.NEGATIVE_INFINITY;
        }

        double[] newValue = getNodeTrait(node);

        if (latentLiability instanceof DummyLatentTruncationProvider) {
            return Double.POSITIVE_INFINITY;
        } else {
            return fullDistribution.logPdf(oldValue) - fullDistribution.logPdf(newValue);
        }
    }
 

@Override
public double[][] getJointVariance(final double priorSampleSize,
                                   final double[][] treeVariance, final double[][] treeSharedLengths,
                                   final double[][] traitVariance) {

    double[] eigVals = this.getEigenValuesStrengthOfSelection();
    DenseMatrix64F V = wrap(this.getEigenVectorsStrengthOfSelection(), 0, dim, dim);
    DenseMatrix64F Vinv = new DenseMatrix64F(dim, dim);
    CommonOps.invert(V, Vinv);

    DenseMatrix64F transTraitVariance = new DenseMatrix64F(traitVariance);

    DenseMatrix64F tmp = new DenseMatrix64F(dim, dim);
    CommonOps.mult(Vinv, transTraitVariance, tmp);
    CommonOps.multTransB(tmp, Vinv, transTraitVariance);

    // Computation of matrix
    int ntaxa = tree.getExternalNodeCount();
    double ti;
    double tj;
    double tij;
    double ep;
    double eq;
    double var;
    DenseMatrix64F varTemp = new DenseMatrix64F(dim, dim);
    double[][] jointVariance = new double[dim * ntaxa][dim * ntaxa];
    for (int i = 0; i < ntaxa; ++i) {
        for (int j = 0; j < ntaxa; ++j) {
            ti = treeSharedLengths[i][i];
            tj = treeSharedLengths[j][j];
            tij = treeSharedLengths[i][j];
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    ep = eigVals[p];
                    eq = eigVals[q];
                    var = tij / ep / eq;
                    var += (1 - Math.exp(ep * tij)) * Math.exp(-ep * ti) / ep / ep / eq;
                    var += (1 - Math.exp(eq * tij)) * Math.exp(-eq * tj) / ep / eq / eq;
                    var -= (1 - Math.exp((ep + eq) * tij)) * Math.exp(-ep * ti) * Math.exp(-eq * tj) / ep / eq / (ep + eq);
                    var += (1 - Math.exp(-ep * ti)) * (1 - Math.exp(-eq * tj)) / ep / eq / priorSampleSize;
                    var += 1 / priorSampleSize;
                    varTemp.set(p, q, var * transTraitVariance.get(p, q));
                }
            }
            CommonOps.mult(V, varTemp, tmp);
            CommonOps.multTransB(tmp, V, varTemp);
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    jointVariance[i * dim + p][j * dim + q] = varTemp.get(p, q);
                }
            }
        }
    }
    return jointVariance;
}
 

private double[][] getJointVarianceFull(final double priorSampleSize,
                                        final double[][] treeVariance, final double[][] treeSharedLengths,
                                        final double[][] traitVariance) {

    double[] eigVals = this.getEigenValuesStrengthOfSelection();
    DenseMatrix64F V = wrap(this.getEigenVectorsStrengthOfSelection(), 0, dim, dim);
    DenseMatrix64F Vinv = new DenseMatrix64F(dim, dim);
    CommonOps.invert(V, Vinv);

    DenseMatrix64F transTraitVariance = new DenseMatrix64F(traitVariance);

    DenseMatrix64F tmp = new DenseMatrix64F(dim, dim);
    CommonOps.mult(Vinv, transTraitVariance, tmp);
    CommonOps.multTransB(tmp, Vinv, transTraitVariance);

    // inverse of eigenvalues
    double[][] invEigVals = new double[dim][dim];
    for (int p = 0; p < dim; ++p) {
        for (int q = 0; q < dim; ++q) {
            invEigVals[p][q] = 1 / (eigVals[p] + eigVals[q]);
        }
    }

    // Computation of matrix
    int ntaxa = tree.getExternalNodeCount();
    double ti;
    double tj;
    double tij;
    double ep;
    double eq;
    DenseMatrix64F varTemp = new DenseMatrix64F(dim, dim);
    double[][] jointVariance = new double[dim * ntaxa][dim * ntaxa];
    for (int i = 0; i < ntaxa; ++i) {
        for (int j = 0; j < ntaxa; ++j) {
            ti = treeSharedLengths[i][i];
            tj = treeSharedLengths[j][j];
            tij = treeSharedLengths[i][j];
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    ep = eigVals[p];
                    eq = eigVals[q];
                    varTemp.set(p, q, Math.exp(-ep * ti) * Math.exp(-eq * tj) * (invEigVals[p][q] * (Math.exp((ep + eq) * tij) - 1) + 1 / priorSampleSize) * transTraitVariance.get(p, q));
                }
            }
            CommonOps.mult(V, varTemp, tmp);
            CommonOps.multTransB(tmp, V, varTemp);
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    jointVariance[i * dim + p][j * dim + q] = varTemp.get(p, q);
                }
            }
        }
    }
    return jointVariance;
}
 

private double[][] getJointVarianceDiagonal(final double priorSampleSize,
                                            final double[][] treeVariance, final double[][] treeSharedLengths,
                                            final double[][] traitVariance) {

    // Eigen of strength of selection matrix
    double[] eigVals = this.getEigenValuesStrengthOfSelection();
    int ntaxa = tree.getExternalNodeCount();
    double ti;
    double tj;
    double tij;
    double ep;
    double eq;
    double var;
    DenseMatrix64F varTemp = new DenseMatrix64F(dim, dim);
    double[][] jointVariance = new double[dim * ntaxa][dim * ntaxa];
    for (int i = 0; i < ntaxa; ++i) {
        for (int j = 0; j < ntaxa; ++j) {
            ti = treeSharedLengths[i][i];
            tj = treeSharedLengths[j][j];
            tij = treeSharedLengths[i][j];
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    ep = eigVals[p];
                    eq = eigVals[q];
                    if (ep + eq == 0.0) {
                        var = (tij + 1 / priorSampleSize) * traitVariance[p][q];
                    } else {
                        var = Math.exp(-ep * ti) * Math.exp(-eq * tj) * (Math.expm1((ep + eq) * tij) / (ep + eq) + 1 / priorSampleSize) * traitVariance[p][q];
                    }
                    varTemp.set(p, q, var);
                }
            }
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    jointVariance[i * dim + p][j * dim + q] = varTemp.get(p, q);
                }
            }
        }
    }
    return jointVariance;
}
 

private static boolean allZeroOrInfinite(DenseMatrix64F M) {
    for (int i = 0; i < M.getNumElements(); i++) {
        if (Double.isFinite(M.get(i)) && M.get(i) != 0.0) return false;
    }
    return true;
}
 

/**
 * Computes a basis (the principle components) from the most dominant eigenvectors.
 */
public void computeBasis() {
	if (sampleIndex != numSamples)
		throw new IllegalArgumentException("Not all the data has been added");
	if (numComponents > numSamples)
		throw new IllegalArgumentException(
				"More data needed to compute the desired number of components");

	means = new DenseMatrix64F(sampleSize, 1);
	// compute the mean of all the samples
	for (int i = 0; i < numSamples; i++) {
		for (int j = 0; j < sampleSize; j++) {
			double val = means.get(j);
			means.set(j, val + A.get(i, j));
		}
	}
	for (int j = 0; j < sampleSize; j++) {
		double avg = means.get(j) / numSamples;
		means.set(j, avg);
	}

	// subtract the mean from the original data
	for (int i = 0; i < numSamples; i++) {
		for (int j = 0; j < sampleSize; j++) {
			A.set(i, j, A.get(i, j) - means.get(j));
		}
	}

	// compute SVD and save time by not computing U
	SingularValueDecomposition<DenseMatrix64F> svd = DecompositionFactory.svd(numSamples, sampleSize,
			false, true, compact);
	if (!svd.decompose(A))
		throw new RuntimeException("SVD failed");

	V_t = svd.getV(null, true);
	W = svd.getW(null);

	// singular values are in an arbitrary order initially and need to be sorted in descending order
	SingularOps.descendingOrder(null, false, W, V_t, true);

	// strip off unneeded components and find the basis
	V_t.reshape(numComponents, sampleSize, true);

}