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

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) {
    int dim = destination.numRows;

    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {

            double n = Math.abs(numeratorMatrix.get(i, j));
            double d = Math.abs(otherMatrix.get(i, j));

            if (n == 0 && d == 0) {
                destination.set(i, j, 0);
            } else {
                destination.set(i, j, n / (n + d));
            }

        }
    }

}
 

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

        }
    }
}
 

private static Complex64F[] findRoots(double... coefficients) {
    int N = coefficients.length - 1;

    // Construct the companion matrix. This is a square N x N matrix.
    final DenseMatrix64F c = new DenseMatrix64F(N, N);

    double a = coefficients[N];
    for (int i = 0; i < N; i++) {
        c.set(i, N - 1, -coefficients[i] / a);
    }
    for (int i = 1; i < N; i++) {
        c.set(i, i - 1, 1);
    }

    // Use generalized eigenvalue decomposition to find the roots.
    EigenDecomposition<DenseMatrix64F> evd = DecompositionFactory.eig(N, false);

    evd.decompose(c);

    final Complex64F[] roots = new Complex64F[N];

    for (int i = 0; i < N; i++) {
        roots[i] = evd.getEigenvalue(i);
    }

    return roots;
}
 

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

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

@Override
public double[] getTipPartial(int taxonIndex, boolean fullyObserved) {

    assert (numTraits == 1);
    assert (samplingPrecision.rows() == dimTrait && samplingPrecision.columns() == dimTrait);

    recomputeVariance();

    if (fullyObserved) {
        return new double[dimTrait + 1];
    }

    double[] partial = super.getTipPartial(taxonIndex, fullyObserved);
    DenseMatrix64F V = MissingOps.wrap(partial, dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);

    //TODO: remove diagonalOnly part
    if (diagonalOnly) {
        for (int index = 0; index < dimTrait; index++) {
            V.set(index, index, V.get(index, index) + 1 / samplingPrecision.component(index, index));
        }
    } else {
        for (int i = 0; i < dimTrait; i++) {
            for (int j = 0; j < dimTrait; j++) {
                V.set(i, j, V.get(i, j) + samplingVariance.component(i, j));
            }
        }
    }


    DenseMatrix64F P = new DenseMatrix64F(dimTrait, dimTrait);
    MissingOps.safeInvert2(V, P, false); //TODO this isn't necessary when this is fully observed

    MissingOps.unwrap(P, partial, dimTrait);
    MissingOps.unwrap(V, partial, dimTrait + dimTrait * dimTrait);

    if (DEBUG) {
        System.err.println("taxon " + taxonIndex);
        System.err.println("\tprecision: " + P);
        System.err.println("\tmean: " + new WrappedVector.Raw(partial, 0, dimTrait));
    }

    return partial;
}
 

private static void idMinusA(DenseMatrix64F A) {
    CommonOps.scale(-1.0, A);
    for (int i = 0; i < A.numCols; i++) {
        A.set(i, i, 1.0 + A.get(i, i));
    }
}
 

/**
 * 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);

}
 

/**
 * Loads the PCA matrix, means and eigenvalues matrix (if {@link #doWhitening} is true) from the given
 * file. The file is supposed to be generated by the {@link #savePCAToFile(String)} method.
 * 
 * @param PCAFileName
 *            the PCA file
 * @throws Exception
 */
public void loadPCAFromFile(String PCAFileName) throws Exception {
	// parse the PCA projection file and put the PCA components in a 2-d double array
	BufferedReader in = new BufferedReader(new FileReader(PCAFileName));
	String line = "";

	// parse the first line which contains the training sample means
	line = in.readLine();
	String[] meanString = line.trim().split(" ");
	if (meanString.length != sampleSize) {
		throw new Exception("Means line is wrong!");
	}
	means = new DenseMatrix64F(sampleSize, 1);
	for (int j = 0; j < sampleSize; j++) {
		means.set(j, Double.parseDouble(meanString[j]));
	}

	// parse the first line which contains the eigenvalues and initialize the diagonal eigenvalue matrix W
	line = in.readLine();
	if (doWhitening) {
		String[] eigenvaluesString = line.trim().split(" ");
		if (eigenvaluesString.length < numComponents) {
			throw new Exception("Eigenvalues line is wrong!");
		}
		W = new DenseMatrix64F(numComponents, numComponents);
		for (int i = 0; i < numComponents; i++) {
			// transform the eigenValues
			double eigenvalue = Double.parseDouble(eigenvaluesString[i]);
			eigenvalue = Math.pow(eigenvalue, -0.5);
			W.set(i, i, eigenvalue);
		}
	}

	V_t = new DenseMatrix64F(numComponents, sampleSize);
	for (int i = 0; i < numComponents; i++) {
		if (i % 100 == 0) {
			System.out.println(i + " PCA components loaded.");
		}
		try {
			line = in.readLine();
		} catch (IOException e) {
			throw new Exception(
					"Check whether the given PCA matrix contains the correct number of components!");
		}
		String[] componentString = line.trim().split(" ");
		for (int j = 0; j < sampleSize; j++) {
			double componentElement = Double.parseDouble(componentString[j]);
			V_t.set(i, j, componentElement);
		}
	}

	// if whitening is true then whiten the PCA matrix V_t by multiplying it with W
	if (doWhitening) {
		System.out.print("Whitening the PCA matrix..");
		DenseMatrix64F V_t_w = new DenseMatrix64F(numComponents, sampleSize);
		CommonOps.mult(W, V_t, V_t_w);
		V_t = V_t_w;
	}

	in.close();

	isPcaInitialized = true;
}