org.ojalgo.matrix.PrimitiveMatrix Java Examples

/**
 * Creates a design matrix used for least squares fitting from a collection of
 * points.
 *
 * @see #solveVector(Collection)
 * @param points points in a 3D space.
 * @return a [points.size()][9] matrix of real values.
 */
private static PrimitiveMatrix createDesignMatrix(
	final Collection<Vector3d> points)
{
	final BasicMatrix.Builder<PrimitiveMatrix> builder = PrimitiveMatrix.FACTORY
		.getBuilder(points.size(), MIN_DATA);
	final Iterator<Vector3d> iterator = points.iterator();
	for (int i = 0; i < points.size(); i++) {
		final Vector3d p = iterator.next();
		builder.set(i, 0, p.x * p.x);
		builder.set(i, 1, p.y * p.y);
		builder.set(i, 2, p.z * p.z);
		builder.set(i, 3, 2 * p.x * p.y);
		builder.set(i, 4, 2 * p.x * p.z);
		builder.set(i, 5, 2 * p.y * p.z);
		builder.set(i, 6, 2 * p.x);
		builder.set(i, 7, 2 * p.y);
		builder.set(i, 8, 2 * p.z);
	}
	return builder.build();
}
 

public static String display(PrimitiveMatrix matrix) {
    long row = matrix.countRows();
    long column = matrix.countColumns();
    StringBuilder sb = new StringBuilder();
    sb.append("size = ").append(matrix.countRows()).append("x")
            .append(matrix.countColumns()).append("\n");
    for (int i=0;i<row;i++){
        for (int j=0;j<column;j++){
            sb.append(matrix.get(i,j));
            if (j!=column-1){
                sb.append(", ");
            }
        }
        sb.append("\n");
    }
    return sb.toString();
}
 

/**
 * Solves the equation for the quadratic surface that best fits the given
 * points.
 * <p>
 * The vector solved is the polynomial Ax<sup>2</sup> + By<sup>2</sup> +
 * Cz<sup>2</sup> + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz, i.e. the general
 * equation of a quadric. The fitting is done with least squares.
 * </p>
 *
 * @param points A collection of points in a 3D space.
 * @return the solution vector of the surface.
 */
private static double[] solveVector(final Collection<Vector3d> points) {
	final int n = points.size();
	// Find (dT * d)^-1
	final PrimitiveMatrix d = createDesignMatrix(points);
	final PrimitiveMatrix dT = d.transpose();
	final PrimitiveMatrix dTDInv = dT.multiply(d).invert();
	// Multiply dT * O, where O = [1, 1, ... 1] (n x 1) matrix
	final PrimitiveMatrix o = PrimitiveMatrix.FACTORY.makeFilled(n, 1,
		new Deterministic(1.0));
	final PrimitiveMatrix dTO = dT.multiply(o);
	// Find solution A = (dT * d)^-1 * (dT * O)
	return dTDInv.multiply(dTO).toRawCopy1D();
}
 

public void eigenValueDecomposition(double[][] data) {
	System.out.println("Performing eigenvector decomposition on " + data.length + " x " + data[data.length - 1].length + " matrix ");
	PrimitiveMatrix matrix = PrimitiveMatrix.FACTORY.rows(data);


	eig = Eigenvalue.make(matrix, true);

	if (!eig.decompose(matrix)) {
		throw new RuntimeException("Decomposition failed");
	}

	// check if eigenvalues are in descending order
	eigenvalues = eig.getEigenvalues().toRawCopy1D();

	for (int i = 0; i < eigenvalues.length - 1; i++) {
		if (eigenvalues[i] < eigenvalues[eigenvalues.length - 1]) {
			ascendingorder = true;
		}
	}


	if (ascendingorder) {
		System.out.println("WARNING: eigenvalues are in ascending order. Will flip them for you.");
		// invert eigenvalues
		double[] tmpeig = new double[eigenvalues.length];
		for (int d = 0; d < eigenvalues.length; d++) {
			tmpeig[d] = eigenvalues[eigenvalues.length - 1 - d];
		}
		eigenvalues = tmpeig;
	} else {
		System.out.println("Eigenvalues are in descending order. ");
	}

	eigenValueMatrix = eig.getV();
}