Java Code Examples for org.apache.commons.math3.linear.RealMatrix#setSubMatrix()

@Nonnull
public static RealMatrix combinedMatrices(@Nonnull final RealMatrix[][] grid,
        final int dimensions) {
    Preconditions.checkArgument(grid.length >= 1, "The number of rows must be greater than 1");
    Preconditions.checkArgument(grid[0].length >= 1,
        "The number of cols must be greater than 1");
    Preconditions.checkArgument(dimensions > 0, "Dimension should be more than 0: ",
        dimensions);

    final int rows = grid.length;
    final int cols = grid[0].length;

    final RealMatrix combined = new BlockRealMatrix(rows * dimensions, cols * dimensions);
    for (int row = 0; row < grid.length; row++) {
        for (int col = 0; col < grid[row].length; col++) {
            combined.setSubMatrix(grid[row][col].getData(), row * dimensions, col * dimensions);
        }
    }
    return combined;
}
 

@Nonnull
public static RealMatrix combinedMatrices(@Nonnull final RealMatrix[] grid) {
    Preconditions.checkArgument(grid.length >= 1,
        "The number of rows must be greater than 0: " + grid.length);

    final int rows = grid.length;
    final int rowDims = grid[0].getRowDimension();
    final int colDims = grid[0].getColumnDimension();

    final RealMatrix combined = new BlockRealMatrix(rows * rowDims, colDims);
    for (int row = 0; row < grid.length; row++) {
        RealMatrix cell = grid[row];
        Preconditions.checkArgument(cell.getRowDimension() == rowDims,
            "Mismatch in row dimensions at row ", row);
        Preconditions.checkArgument(cell.getColumnDimension() == colDims,
            "Mismatch in col dimensions at row ", row);
        combined.setSubMatrix(cell.getData(), row * rowDims, 0);
    }
    return combined;
}
 

/**
 * QR decomposition for a tridiagonal matrix T.
 *
 * @see https://gist.github.com/lightcatcher/8118181
 * @see http://www.ericmart.in/blog/optimizing_julia_tridiag_qr
 * @param T target tridiagonal matrix
 * @param R output matrix for R which is the same shape as T
 * @param Qt output matrix for Q.T which is the same shape an T
 */
public static void tridiagonalQR(@Nonnull final RealMatrix T, @Nonnull final RealMatrix R,
        @Nonnull final RealMatrix Qt) {
    int n = T.getRowDimension();
    Preconditions.checkArgument(n == R.getRowDimension() && n == R.getColumnDimension(),
        "T and R must be the same shape");
    Preconditions.checkArgument(n == Qt.getRowDimension() && n == Qt.getColumnDimension(),
        "T and Qt must be the same shape");

    // initial R = T
    R.setSubMatrix(T.getData(), 0, 0);

    // initial Qt = identity
    Qt.setSubMatrix(eye(n), 0, 0);

    for (int i = 0; i < n - 1; i++) {
        // Householder projection for a vector x
        // https://en.wikipedia.org/wiki/Householder_transformation
        RealVector x = T.getSubMatrix(i, i + 1, i, i).getColumnVector(0);
        x = unitL2norm(x);

        RealMatrix subR = R.getSubMatrix(i, i + 1, 0, n - 1);
        R.setSubMatrix(
            subR.subtract(x.outerProduct(subR.preMultiply(x)).scalarMultiply(2)).getData(), i,
            0);

        RealMatrix subQt = Qt.getSubMatrix(i, i + 1, 0, n - 1);
        Qt.setSubMatrix(
            subQt.subtract(x.outerProduct(subQt.preMultiply(x)).scalarMultiply(2)).getData(), i,
            0);
    }
}
 

/**
 * Find eigenvalues and eigenvectors of given tridiagonal matrix T.
 *
 * @see http://web.csulb.edu/~tgao/math423/s94.pdf
 * @see http://stats.stackexchange.com/questions/20643/finding-matrix-eigenvectors-using-qr-
 *      decomposition
 * @param T target tridiagonal matrix
 * @param nIter number of iterations for the QR method
 * @param eigvals eigenvalues are stored here
 * @param eigvecs eigenvectors are stored here
 */
public static void tridiagonalEigen(@Nonnull final RealMatrix T, @Nonnull final int nIter,
        @Nonnull final double[] eigvals, @Nonnull final RealMatrix eigvecs) {
    Preconditions.checkArgument(Arrays.deepEquals(T.getData(), T.transpose().getData()),
        "Target matrix T must be a symmetric (tridiagonal) matrix");
    Preconditions.checkArgument(eigvecs.getRowDimension() == eigvecs.getColumnDimension(),
        "eigvecs must be a square matrix");
    Preconditions.checkArgument(T.getRowDimension() == eigvecs.getRowDimension(),
        "T and eigvecs must be the same shape");
    Preconditions.checkArgument(eigvals.length == eigvecs.getRowDimension(),
        "Number of eigenvalues and eigenvectors must be same");

    int nEig = eigvals.length;

    // initialize eigvecs as an identity matrix
    eigvecs.setSubMatrix(eye(nEig), 0, 0);

    RealMatrix T_ = T.copy();

    for (int i = 0; i < nIter; i++) {
        // QR decomposition for the tridiagonal matrix T
        RealMatrix R = new Array2DRowRealMatrix(nEig, nEig);
        RealMatrix Qt = new Array2DRowRealMatrix(nEig, nEig);
        tridiagonalQR(T_, R, Qt);

        RealMatrix Q = Qt.transpose();
        T_ = R.multiply(Q);
        eigvecs.setSubMatrix(eigvecs.multiply(Q).getData(), 0, 0);
    }

    // diagonal elements correspond to the eigenvalues
    for (int i = 0; i < nEig; i++) {
        eigvals[i] = T_.getEntry(i, i);
    }
}
 

@Override
public Set<Model> flatten(MatrixStack matrixStack) {
	Set<Model> models = new HashSet<>();

	matrixStack.pushMatrix();
	matrixStack.transform(matrix.getMatrix());
	//Create a new model with transformation applied.
	MeshModel transformedModel = clone();
	// correct formula for Normal Matrix is transpose(inverse(mat3(model_mat))
	// we have to augemnt that to 4x4
	RealMatrix normalMatrix3x3 = new LUDecomposition(matrixStack.getMatrix().getSubMatrix(0, 2, 0, 2), 1e-5).getSolver().getInverse().transpose();
	RealMatrix normalMatrix = MatrixUtils.createRealMatrix(4, 4);
	normalMatrix.setSubMatrix(normalMatrix3x3.getData(), 0, 0);
	normalMatrix.setEntry(3, 3, 1);

	transformedModel.faces.stream().forEach(f -> {
			f.normal = TransformUtil.transform(f.normal, normalMatrix);
			f.vertices.forEach(v -> {
				v.vec = matrixStack.apply(v.vec);
				v.normal = v.normal.map(n -> TransformUtil.transform(n, normalMatrix));
			});
		}
	);

	models.add(transformedModel);
	//Flatten child models
	models.addAll(children.stream().flatMap(m -> m.flatten(matrixStack).stream()).collect(Collectors.toSet()));
	matrixStack.popMatrix();
	return models;
}
 

public static RealMatrix augment(RealMatrix matrix, int rows, int columns) {
    if (matrix.getRowDimension() > rows)
        throw new DimensionMismatchException(of("rows: {0} !< {1}"), matrix.getRowDimension(), rows);
    if (matrix.getColumnDimension() > columns)
        throw new DimensionMismatchException(of("columns: {0} !< {1}"), matrix.getColumnDimension(), columns);

    RealMatrix augmented = MatrixUtils.createRealMatrix(rows, columns);
    augmented.setSubMatrix(matrix.getData(), 0, 0);
    return augmented;
}
 

/**
 * Rotates the current matrix
 *
 * @param rotation The rotation to aply
 * @return The rorated matrix
 */
public MatrixStack rotate(Rotation rotation) {
	RealMatrix rotMat = MatrixUtils.createRealMatrix(4, 4);
	rotMat.setSubMatrix(rotation.getMatrix(), 0, 0);
	rotMat.setEntry(3, 3, 1);
	current = current.preMultiply(rotMat);
	return this;
}
 

@Override
public Set<Model> flatten(MatrixStack matrixStack) {
	Set<Model> models = new HashSet<>();

	matrixStack.pushMatrix();
	matrixStack.transform(matrix.getMatrix());
	//Create a new model with transformation applied.
	MeshModel transformedModel = clone();
	// correct formula for Normal Matrix is transpose(inverse(mat3(model_mat))
	// we have to augemnt that to 4x4
	RealMatrix normalMatrix3x3 = new LUDecomposition(matrixStack.getMatrix().getSubMatrix(0, 2, 0, 2), 1e-5).getSolver().getInverse().transpose();
	RealMatrix normalMatrix = MatrixUtils.createRealMatrix(4, 4);
	normalMatrix.setSubMatrix(normalMatrix3x3.getData(), 0, 0);
	normalMatrix.setEntry(3, 3, 1);

	transformedModel.faces.stream().forEach(f -> {
			f.normal = TransformUtil.transform(f.normal, normalMatrix);
			f.vertices.forEach(v -> v.vec = matrixStack.apply(v.vec));
		}
	);

	models.add(transformedModel);
	//Flatten child models
	matrixStack.pushMatrix();
	matrixStack.translate(0.5, 0.5, 0.5);
	models.addAll(children.stream().flatMap(m -> m.flatten(matrixStack).stream()).collect(Collectors.toSet()));
	matrixStack.popMatrix().popMatrix();
	return models;
}
 

@Override
public Set<Model> flatten(MatrixStack matrixStack) {
	Set<Model> models = new HashSet<>();

	matrixStack.pushMatrix();
	matrixStack.transform(matrix.getMatrix());
	//Create a new model with transformation applied.
	MeshModel transformedModel = clone();
	// correct formula for Normal Matrix is transpose(inverse(mat3(model_mat))
	// we have to augemnt that to 4x4
	RealMatrix normalMatrix3x3 = new LUDecomposition(matrixStack.getMatrix().getSubMatrix(0, 2, 0, 2), 1e-5).getSolver().getInverse().transpose();
	RealMatrix normalMatrix = MatrixUtils.createRealMatrix(4, 4);
	normalMatrix.setSubMatrix(normalMatrix3x3.getData(), 0, 0);
	normalMatrix.setEntry(3, 3, 1);

	transformedModel.faces.stream().forEach(f -> {
			f.normal = TransformUtil.transform(f.normal, normalMatrix);
			f.vertices.forEach(v -> {
				v.vec = matrixStack.apply(v.vec);
				v.normal = v.normal.map(n -> TransformUtil.transform(n, normalMatrix));
			});
		}
	);

	models.add(transformedModel);
	//Flatten child models
	matrixStack.pushMatrix();
	matrixStack.translate(0.5, 0.5, 0.5);
	models.addAll(children.stream().flatMap(m -> m.flatten(matrixStack).stream()).collect(Collectors.toSet()));
	matrixStack.popMatrix().popMatrix();
	return models;
}
 

/**
 * Lanczos tridiagonalization for a symmetric matrix C to make s * s tridiagonal matrix T.
 *
 * @see http://www.cas.mcmaster.ca/~qiao/publications/spie05.pdf
 * @param C target symmetric matrix
 * @param a initial vector
 * @param T result is stored here
 */
public static void lanczosTridiagonalization(@Nonnull final RealMatrix C,
        @Nonnull final double[] a, @Nonnull final RealMatrix T) {
    Preconditions.checkArgument(Arrays.deepEquals(C.getData(), C.transpose().getData()),
        "Target matrix C must be a symmetric matrix");
    Preconditions.checkArgument(C.getColumnDimension() == a.length,
        "Column size of A and length of a should be same");
    Preconditions.checkArgument(T.getRowDimension() == T.getColumnDimension(),
        "T must be a square matrix");

    int s = T.getRowDimension();

    // initialize T with zeros
    T.setSubMatrix(new double[s][s], 0, 0);

    RealVector a0 = new ArrayRealVector(a.length);
    RealVector r = new ArrayRealVector(a);

    double beta0 = 1.d;

    for (int i = 0; i < s; i++) {
        RealVector a1 = r.mapDivide(beta0);
        RealVector Ca1 = C.operate(a1);

        double alpha1 = a1.dotProduct(Ca1);

        r = Ca1.add(a1.mapMultiply(-1.d * alpha1)).add(a0.mapMultiply(-1.d * beta0));

        double beta1 = r.getNorm();

        T.setEntry(i, i, alpha1);
        if (i - 1 >= 0) {
            T.setEntry(i, i - 1, beta0);
        }
        if (i + 1 < s) {
            T.setEntry(i, i + 1, beta1);
        }

        a0 = a1.copy();
        beta0 = beta1;
    }
}
 

private void computeADFStatistics() {
	double[] y = diff(ts);
	RealMatrix designMatrix = null;
	int k = lag+1;
	int n = ts.length - 1;
	
	RealMatrix z = MatrixUtils.createRealMatrix(laggedMatrix(y, k)); //has rows length(ts) - 1 - k + 1
	RealVector zcol1 = z.getColumnVector(0); //has length length(ts) - 1 - k + 1
	double[] xt1 = subsetArray(ts, k-1, n-1);  //ts[k:(length(ts) - 1)], has length length(ts) - 1 - k + 1
	double[] trend = sequence(k,n); //trend k:n, has length length(ts) - 1 - k + 1
	if (k > 1) {
		RealMatrix yt1 = z.getSubMatrix(0, ts.length - 1 - k, 1, k-1); //same as z but skips first column
		//build design matrix as cbind(xt1, 1, trend, yt1)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3 + k - 1);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
		designMatrix.setSubMatrix(yt1.getData(), 0, 3);
		
	} else {
		//build design matrix as cbind(xt1, 1, tt)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
	}
	/*OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
	regression.setNoIntercept(true);
	regression.newSampleData(zcol1.toArray(), designMatrix.getData());
	double[] beta = regression.estimateRegressionParameters();
	double[] sd = regression.estimateRegressionParametersStandardErrors();
	*/
	RidgeRegression regression = new RidgeRegression(designMatrix.getData(), zcol1.toArray());
	regression.updateCoefficients(.0001);
	double[] beta = regression.getCoefficients();
	double[] sd = regression.getStandarderrors();
	
	double t = beta[0] / sd[0];
	if (t <= PVALUE_THRESHOLD) {
		this.needsDiff = true;
	} else {
		this.needsDiff = false;
	}	
}