# Disclaimer

Some things that I must state:

1. This is just for practice. I want to know the various algorithms in linear algebra in detail.
2. I will not use this in a practical scenario. Except for verifying my assignment solutions. Consequently, there is a slim probability that I’ll be provided with inputs that lead to numerical instability. However, I’d like my code to be able to handle it.

# The code

## Matrix.java

Contains the entire codebase. Probably needs to be split up.

package com.github.subh0m0y.matrix;

import java.util.Arrays;
import java.util.Objects;
import java.util.Random;

import static com.github.subh0m0y.matrix.Standards.EPSILON;

@SuppressWarnings("WeakerAccess")
public class Matrix {
public static final String FORMAT_STRING = "%+.2e";
private final int rows;
private final int cols;
private final double[][] data;

/**
* Creates a new matrix with the given data
*
* @param data         The raw data to encapsulate.
* @param makeDeepCopy Whether to copy every element or use the given reference.
*/
Matrix(final double[][] data, final boolean makeDeepCopy) {
this.rows = data.length;
this.cols = this.rows == 0 ? 0 : data[0].length;

if (makeDeepCopy) {
this.data = new double[rows][cols];
for (int i = 0; i < rows; i++) {
System.arraycopy(data[i], 0, this.data[i], 0, cols);
}
} else {
this.data = data;
}
}

/**
* Creates a new zero matrix of order rows x cols.
*
* @param rows The number of rows required.
* @param cols The number of columns required.
*/
Matrix(final int rows, final int cols) {
if (rows < 0) {
throw new IllegalArgumentException("Invalid number of rows : " + rows);
}
if (cols < 0) {
throw new IllegalArgumentException("Invalid number of columns : " + cols);
}
this.rows = rows;
this.cols = cols;
data = new double[rows][cols];
}

/**
* Creates a new Matrix that is a copy of the given
* matrix and is independent of the original.
*
* @param matrix The matrix to copy.
*/
public Matrix(final Matrix matrix) {
rows = matrix.rows;
cols = matrix.cols;
data = new double[rows][cols];
for (int i = 0; i < rows; i++) {
System.arraycopy(matrix.data[i], 0, data[i], 0, cols);
}
}

public static Matrix fromArray(final double[][] data) {
return new Matrix(data, true);
}

public static Matrix zero(final int rows, final int cols) {
return new Matrix(rows, cols);
}

public static Matrix identity(final int order) {
double[][] data = new double[order][order];
for (int i = 0; i < order; i++) {
data[i][i] = 1;
}
return new Matrix(data, false);
}

public static Matrix fromLinearArray(final int rows, final int cols,
final double... elements) {
if (elements.length != rows * cols) {
throw new IllegalArgumentException("Invalid number of elements: " + elements.length
+ " Expected: " + rows * cols);
}
double[][] data = new double[rows][cols];
for (int i = 0; i < rows; i++) {
System.arraycopy(elements, i * cols, data[i], 0, cols);
}
return new Matrix(data, false);
}

public static Matrix random(final int rows, final int cols) {
Matrix matrix = new Matrix(rows, cols);
Random random = new Random();
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
matrix.data[i][j] = random.nextDouble();
}
}
return matrix;
}

public int getRows() {
return rows;
}

public int getCols() {
return cols;
}

public double get(final int i, final int j) throws IndexOutOfBoundsException {
if (i < 0 || i >= rows) {
throw new IndexOutOfBoundsException("Invalid row index : " + i);
}
if (j < 0 || j >= cols) {
throw new IndexOutOfBoundsException("Invalid column index : " + i);
}
return data[i][j];
}

public double[] getRow(final int row) throws IllegalArgumentException {
throwIfInvalidIndex(row, rows, "row");
return Arrays.copyOf(data[row], cols);
}

public double[] getColumn(final int column) throws IllegalArgumentException {
throwIfInvalidIndex(column, cols, "column");
double[] col = new double[rows];
for (int i = 0; i < rows; i++) {
col[i] = data[i][column];
}
return col;
}

@Override
public boolean equals(Object o) {
if (this == o) return true;
if (!(o instanceof Matrix)) return false;
Matrix matrix = (Matrix) o;
return rows == matrix.rows &&
cols == matrix.cols &&
Arrays.deepEquals(data, matrix.data);
}

@Override
public int hashCode() {
int result = Objects.hash(rows, cols);
result = 31 * result + Arrays.deepHashCode(data);
return result;
}

@Override
public String toString() {
if (rows == 0 || cols == 0) {
return "";
}
StringBuilder builder = new StringBuilder();
for (int i = 0; i < rows; i++) {
if (i == 0) {
builder.append(" /");
} else if (i == rows - 1) {
builder.append(" \");
} else {
builder.append("| ");
}
builder.append(String.format(FORMAT_STRING, data[i][0]));
for (int j = 1; j < cols; j++) {
builder.append(", ");
builder.append(String.format(FORMAT_STRING, data[i][j]));
}
if (i == 0) {
builder.append("\");
} else if (i == rows - 1) {
builder.append("/");
} else {
builder.append(" |");
}
builder.append("n");
}
return builder.toString();
}

public boolean isZero() {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
if (data[i][j] != 0) {
return false;
}
}
}
return true;
}

public boolean isSquare() {
return rows == cols;
}

public boolean isUpperTriangular() {
if (!isSquare()) {
return false;
}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < i; j++) {
if (Math.abs(data[i][j]) > EPSILON) {
return false;
}
}
}
return true;
}

public boolean isLowerTriangular() {
if (!isSquare()) {
return false;
}
for (int i = 0; i < rows; i++) {
for (int j = i + 1; j < cols; j++) {
if (Math.abs(data[i][j]) > EPSILON) {
return false;
}
}
}
return true;
}

public boolean isDiagonal() {
return isLowerTriangular() && isUpperTriangular();
}

public boolean isIdentity() {
if (!isDiagonal()) {
return false;
}
for (int i = 0; i < rows; i++) {
if (Math.abs(data[i][i] - 1) > EPSILON) {
return false;
}
}
return true;
}

public void transposeInPlace() {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < i; j++) {
double temp = data[i][j];
data[i][j] = data[j][i];
data[j][i] = temp;
}
}
}

public Matrix transpose() {
Matrix transpose = new Matrix(rows, cols);
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
transpose.data[j][i] = data[i][j];
}
}
return transpose;
}

public void zeroFill() {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
if (Math.abs(data[i][j]) < EPSILON) {
data[i][j] = 0;
}
}
}
}

public void scaleInPlace(final double scale) {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
data[i][j] *= scale;
}
}
}

public Matrix scale(final double scale) {
Matrix temp = new Matrix(this);
temp.scaleInPlace(scale);
return temp;
}

private void throwIncompatible(String operation) throws IllegalArgumentException {
throw new IllegalArgumentException("Given matrix is not compatible with the invoking matrix for "
+ operation);
}

}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
}
}
}

Matrix temp = new Matrix(this);
return temp;
}

public void subtractInPlace(final Matrix addend) throws IllegalArgumentException {
throwIncompatible("subtraction");
}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
}
}
}

public Matrix subtract(final Matrix addend) throws IllegalArgumentException {
Matrix temp = new Matrix(this);
return temp;
}

public void multiplyInPlace(final Matrix multiplicand) throws IllegalArgumentException {
if (multiplicand.rows != rows || multiplicand.cols != cols) {
}
double[][] product = new double[rows][cols];
storeProduct(multiplicand, product);
for (int i = 0; i < rows; i++) {
System.arraycopy(product[i], 0, data[i], 0, cols);
}
}

public Matrix multiply(final Matrix multiplicand) throws IllegalArgumentException {
if (multiplicand.rows != cols) {
throwIncompatible("multiplication");
}
double[][] product = new double[rows][cols];
storeProduct(multiplicand, product);
return new Matrix(product, false);
}

private void storeProduct(Matrix multiplicand, double[][] product) {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < multiplicand.cols; j++) {
product[i][j] = 0;
for (int k = 0; k < cols; k++) {
product[i][j] += data[i][k] * multiplicand.data[k][j];
}
}
}
}

public void elementMultiplyInPlace(final Matrix matrix) throws IllegalArgumentException {
if (matrix.rows != rows || matrix.cols != cols) {
throwIncompatible("element-wise multiplication");
}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
data[i][j] *= matrix.data[i][j];
}
}
}

public Matrix elementMultiply(final Matrix matrix) throws IllegalArgumentException {
Matrix temp = new Matrix(this);
temp.elementMultiplyInPlace(matrix);
return temp;
}

public void elementDivideInPlace(final Matrix matrix) throws IllegalArgumentException {
if (matrix.rows != rows || matrix.cols != cols) {
throwIncompatible("element-wise multiplication");
}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
data[i][j] /= matrix.data[i][j];
}
}
}

public Matrix elementDivide(final Matrix matrix) throws IllegalArgumentException {
Matrix temp = new Matrix(this);
temp.elementDivideInPlace(matrix);
return temp;
}

public Matrix exponentiate(int power) throws IllegalArgumentException {
if (!isSquare()) {
throwIncompatible("exponentiation");
}
if (power < 0) {
throw new IllegalArgumentException("Power cannot be negative.");
}
Matrix product = Matrix.identity(rows);
Matrix x = new Matrix(this);
while (power > 0) {
if ((power & 1) == 1) {
product.multiplyInPlace(x);
}
x.multiplyInPlace(x);
power >>= 1;
}
return product;
}

public boolean isOrthogonal() {
Matrix value = multiply(transpose());
value.zeroFill();
return value.isIdentity();
}

public boolean isInvolutary() {
Matrix value = multiply(this);
value.zeroFill();
return value.isIdentity();
}

public Matrix appendRight(final Matrix matrix) {
if (matrix.rows != rows) {
throwIncompatible("appending right");
}
Matrix store = new Matrix(rows, cols + matrix.cols);
for (int i = 0; i < rows; i++) {
System.arraycopy(data[i], 0, store.data[i], 0, cols);
System.arraycopy(matrix.data[i], 0, store.data[i], cols, matrix.cols);
}
return store;
}

public Matrix appendBottom(final Matrix matrix) throws IllegalArgumentException {
if (matrix.cols != cols) {
throwIncompatible("appending right");
}
Matrix store = new Matrix(rows + matrix.rows, cols);
for (int i = 0; i < rows; i++) {
System.arraycopy(data[i], 0, store.data[i], 0, cols);
}
for (int i = 0; i < matrix.rows; i++) {
System.arraycopy(matrix.data[i], 0, store.data[i + rows], 0, cols);
}
return store;
}

private void throwIfInvalidIndex(int value, int limit, String quantity) throws IllegalArgumentException {
if (value <= 0 || value > limit) {
throw new IllegalArgumentException("Invalid " + quantity + " index : " + value);
}
}

public Matrix[] splitAtColumn(final int column) throws IllegalArgumentException {
throwIfInvalidIndex(column, cols, "column");
int columnResidue = cols - column;
Matrix left = new Matrix(rows, column);
Matrix right = new Matrix(rows, columnResidue);
for (int i = 0; i < rows; i++) {
System.arraycopy(data[i], 0, left.data[i], 0, column);
System.arraycopy(data[i], column, right.data[i], 0, columnResidue);
}
return new Matrix[]{left, right};
}

public Matrix[] splitAtRow(final int row) throws IllegalArgumentException {
throwIfInvalidIndex(row, rows, "row");
int rowResidue = rows - row;
Matrix top = new Matrix(row, cols);
Matrix bottom = new Matrix(rowResidue, cols);
for (int i = 0; i < row; i++) {
System.arraycopy(data[i], 0, top.data[i], 0, cols);
}
for (int i = 0; i < rowResidue; i++) {
System.arraycopy(data[i + row], 0, bottom.data[i], 0, cols);
}
return new Matrix[]{top, bottom};
}

public void swapRowsInPlace(final int row1, final int row2) throws IllegalArgumentException {
throwIfInvalidIndex(row1, rows, "row");
throwIfInvalidIndex(row2, rows, "row");
double[] row = new double[cols];
System.arraycopy(data[row1], 0, row, 0, cols);
System.arraycopy(data[row2], 0, data[row1], 0, cols);
System.arraycopy(row, 0, data[row2], 0, cols);
}

public Matrix swapRows(final int row1, final int row2) throws IllegalArgumentException {
Matrix temp = new Matrix(this);
temp.swapRowsInPlace(row1, row2);
return temp;
}

public void swapColumnsInPlace(int col1, int col2) throws IllegalArgumentException {
throwIfInvalidIndex(col1, rows, "column");
throwIfInvalidIndex(col2, rows, "column");
double[] column = new double[rows];
for (int i = 0; i < rows; i++) {
column[i] = data[i][col1];
}
for (int i = 0; i < rows; i++) {
data[i][col1] = data[i][col2];
}
for (int i = 0; i < rows; i++) {
data[i][col2] = column[i];
}
}

public Matrix swapColumns(final int col1, final int col2) throws IllegalArgumentException {
Matrix temp = new Matrix(this);
temp.swapColumnsInPlace(col1, col2);
return temp;
}

private void checkLen(double[] doubles, int len) {
if (doubles.length != len) {
throw new IllegalArgumentException(
"Invalid number of elements. Expected : " + len + " Found : " + doubles.length
);
}
}

private void addScaledRow(double scale, double[] row, int index) {
for (int i = 0; i < cols; i++) {
data[index][i] += row[i] * scale;
}
}

private void scaleRow(double scale, int index) {
for (int i = 0; i < cols; i++) {
data[index][i] *= scale;
}
}

public void convertToReducedRowEchelon() {
int limit = Math.min(rows, cols);
// Cascade from the top left corner downwards and rightwards
for (int i = 0; i < limit; i++) {
double factor = data[i][i];
for (int j = i + 1; j < rows; j++) {
// Make all the rows below the ith row have zeros in the ith column
}
// Normalize to make sure that the leading number in every row is 1
scaleRow(1 / data[i][i], i);
}
}

public void convertEchelonToNormal() {
int limit = Math.min(rows, cols);
// Cascade from the bottom right part of original matrix upwards and leftwards
for (int i = limit - 1; i >= 0; i--) {
double factor = data[i][i];
for (int j = i - 1; j >= 0; j--) {
// Make all the rows above the ith row have zeros in the ith column
}
}
}

private Matrix rowReducedForm = null;

public int getRank() {
if (rowReducedForm == null) {
rowReducedForm = new Matrix(this);
rowReducedForm.convertToReducedRowEchelon();
}
int zeroRows = 0, limit = Math.min(rows, cols);
outer:
for (double[] row : rowReducedForm.data) {
for (int i = 0; i < limit; i++) {
if (Math.abs(row[i]) > EPSILON) {
continue outer;
}
}
zeroRows++;
}
return limit - zeroRows;
}

private Matrix inverse = null;

public Matrix getInverse() throws ArithmeticException {
if (inverse != null) {
return inverse;
}
if (!isSquare()) {
throw new ArithmeticException("Cannot find inverse of a non-square matrix.");
}
Matrix augmented = appendRight(Matrix.identity(rows));
int rank = augmented.getRank();
if (rank < rows) {
throw new ArithmeticException("Cannot find inverse of singular matrix.");
}
augmented.convertEchelonToNormal();
Matrix[] matrices = augmented.splitAtColumn(cols);
rowReducedForm = matrices[0];
inverse = matrices[1];
return inverse;
}
}


## Standards.java

Contains accessible constants used throughout the system. Will be expanded. Or removed completely.

package com.github.subh0m0y.matrix;

public class Standards {
public static double EPSILON = 1.0e-14;
}


# Test cases

There should probably be more cases with more corner cases. I used the TestNG framework here.

## BasicMatrixTest.java

package com.github.subh0m0y.matrix;

import org.testng.annotations.BeforeMethod;
import org.testng.annotations.Test;

import java.util.Random;

import static org.testng.Assert.*;

public class BasicMatrixTest {
private static final int ROWS = 100;
private static final int COLS = 100;
private static Random random;

@BeforeMethod
public void setUp() {
random = new Random();
}

private static void changeElement(double[][] data, Random random) {
data[random.nextInt(ROWS)][random.nextInt(COLS)] *= random.nextGaussian() * 2;
data[random.nextInt(ROWS)][random.nextInt(COLS)] += random.nextGaussian() + 1;
}

@Test
public void testEqualsAndHashCode() {
double[][] data = new double[ROWS][COLS];
Utilities.populate(data, random);
Matrix matrix1 = new Matrix(data, false);
Matrix matrix2 = new Matrix(matrix1);
assertTrue(matrix1.equals(matrix2));
assertEquals(matrix1.hashCode(), matrix2.hashCode());
changeElement(data, random);
assertFalse(matrix1.equals(matrix2));
assertNotEquals(matrix1.hashCode(), matrix2.hashCode());
}

@Test
public void testConstructors() {
double[][] data = new double[ROWS][COLS];
Matrix matrixZero1 = new Matrix(ROWS, COLS);
Matrix matrixZero2 = new Matrix(data, true);

assertTrue(matrixZero1.isZero());
assertTrue(matrixZero2.isZero());
assertEquals(matrixZero1, matrixZero2);

Utilities.populate(data, random);

Matrix matrix1 = new Matrix(data, true);
Matrix matrix2 = new Matrix(data, false);
assertEquals(matrix1, matrix2);

changeElement(data, random);
assertNotEquals(matrix1, matrix2);

Matrix matrix3 = new Matrix(matrix1);
assertEquals(matrix1, matrix3);
}

@Test
public void testFromArray() {
double[][] data = new double[ROWS][COLS];
Utilities.populate(data, random);
double[] linear = new double[ROWS * COLS];
for (int i = 0; i < ROWS; i++) {
System.arraycopy(data[i], 0, linear, i * COLS, COLS);
}
Matrix matrix1 = new Matrix(data, false);
Matrix matrix2 = Matrix.fromArray(data);
Matrix matrix3 = Matrix.fromLinearArray(ROWS, COLS, linear);
assertEquals(matrix1, matrix2);
assertEquals(matrix2, matrix3);
changeElement(data, random);
assertNotEquals(matrix1, matrix2);
assertNotEquals(matrix1, matrix3);
}

@Test
public void testGet() {
double[][] data = new double[ROWS][COLS];
Utilities.populate(data, random);

Matrix matrix1 = new Matrix(data, false);
Matrix matrix2 = new Matrix(matrix1);

assertEquals(matrix1.getRows(), ROWS);
assertEquals(matrix1.getCols(), COLS);
assertEquals(matrix2.getRows(), ROWS);
assertEquals(matrix2.getCols(), COLS);

assertEquals(matrix1.getRows(), matrix2.getRows());
assertEquals(matrix1.getCols(), matrix2.getCols());

for (int i = 0; i < ROWS; i++) {
for (int j = 0; j < COLS; j++) {
assertEquals(matrix1.get(i, j), matrix2.get(i, j));
}
}
assertThrows(() -> matrix1.get(-1, 0));
assertThrows(() -> matrix2.get(0, COLS));
}

@Test
public void testIsZero() {
double[][] data = new double[ROWS][COLS];
Matrix matrix1 = new Matrix(data, false);
assertTrue(matrix1.isZero());
changeElement(data, random);
assertFalse(matrix1.isZero());
Matrix matrix2 = Matrix.zero(ROWS, COLS);
assertTrue(matrix2.isZero());
}

@Test
public void testIsIdentity() {
double[][] data = new double[ROWS][ROWS];
Matrix matrix1 = new Matrix(data, false);
assertFalse(matrix1.isIdentity());
for (int i = 0; i < ROWS; i++) {
data[i][i] = 1;
}
assertTrue(matrix1.isSquare());
assertTrue(matrix1.isIdentity());
changeElement(data, random);
assertFalse(matrix1.isIdentity());

Matrix matrix2 = Matrix.identity(ROWS);
assertTrue(matrix2.isIdentity());
}

@Test
public void testZeroFill() {
double[][] data = new double[ROWS][COLS];
for (int i = 0; i < ROWS; i++) {
for (int j = 0; j < COLS; j++) {
data[i][j] = random.nextDouble() * Standards.EPSILON;
}
}
Matrix matrix = new Matrix(data, false);
assertFalse(matrix.isZero());
matrix.zeroFill();
assertTrue(matrix.isZero());
}

@Test
public void testTranspose() {
double[][] data = new double[ROWS][COLS];
Utilities.populate(data, random);
Matrix matrix1 = new Matrix(data, false);
Matrix matrix2 = matrix1.transpose();
matrix1.transposeInPlace();
assertEquals(matrix1, matrix2);
}

@Test
public void testAppendingAndSplitting() {
Matrix matrix1 = Matrix.random(ROWS, COLS);
Matrix matrix2 = Matrix.random(ROWS, COLS);
Matrix matrixHorizontal = matrix1.appendRight(matrix2);
Matrix matrixVertical = matrix1.appendBottom(matrix2);
Matrix[] horizontal = matrixHorizontal.splitAtColumn(COLS);
Matrix[] vertical = matrixVertical.splitAtRow(ROWS);
assertEquals(matrix1, horizontal[0]);
assertEquals(matrix2, horizontal[1]);
assertEquals(matrix1, vertical[0]);
assertEquals(matrix2, vertical[1]);
}

@Test
public void testSwapAndGetRow() {
Matrix matrix1 = Matrix.random(ROWS, COLS);
Matrix matrix2 = new Matrix(matrix1);
int row1 = random.nextInt(ROWS);
int row2 = random.nextInt(ROWS);
matrix1.swapRowsInPlace(row1, row2);
Matrix matrix3 = matrix2.swapRows(row1, row2);
assertEquals(matrix1, matrix3);
assertEquals(matrix1.getRow(row1), matrix2.getRow(row2));
assertEquals(matrix1.getRow(row2), matrix2.getRow(row1));
}

@Test
public void testSwapAndGetColumn() {
Matrix matrix1 = Matrix.random(ROWS, COLS);
Matrix matrix2 = new Matrix(matrix1);
int col1 = random.nextInt(ROWS);
int col2 = random.nextInt(ROWS);
matrix1.swapColumnsInPlace(col1, col2);
Matrix matrix3 = matrix2.swapColumns(col1, col2);
assertEquals(matrix1, matrix3);
assertEquals(matrix1.getColumn(col1), matrix2.getColumn(col2));
assertEquals(matrix1.getColumn(col2), matrix2.getColumn(col1));
}
}


## MatrixMathTest.java

package com.github.subh0m0y.matrix;

import org.testng.annotations.BeforeMethod;
import org.testng.annotations.Test;

import java.util.Random;

import static org.testng.Assert.*;

public class MatrixMathTest {

private static final int ROWS = 100;
private static final int COLS = 100;
private static final int POWER_BOUND = 1000;
private static final int COUNT = 10;

private static Random random;

@BeforeMethod
public void setUp() {
random = new Random();
}

@Test
double[][] data1 = new double[ROWS][COLS];
double[][] data2 = new double[ROWS][COLS];
Utilities.populate(data1, random);
Utilities.populate(data2, random);
Matrix matrix1 = new Matrix(data1, false);
Matrix matrix2 = new Matrix(data1, false);
assertEquals(matrix1, matrix3);
}

@Test
public void testSubtraction() {
double[][] data1 = new double[ROWS][COLS];
double[][] data2 = new double[ROWS][COLS];
Utilities.populate(data1, random);
Utilities.populate(data2, random);
Matrix matrix1 = new Matrix(data1, false);
Matrix matrix2 = new Matrix(data1, false);
Matrix matrix3 = matrix1.subtract(matrix2);
matrix1.subtractInPlace(matrix2);
assertEquals(matrix1, matrix3);
}

@Test
public void testElementMultiply() {
double[][] data1 = new double[ROWS][COLS];
double[][] data2 = new double[ROWS][COLS];
Utilities.populate(data1, random);
Utilities.populate(data2, random);
Matrix matrix1 = new Matrix(data1, false);
Matrix matrix2 = new Matrix(data1, false);
Matrix matrix3 = matrix1.elementMultiply(matrix2);
matrix1.elementMultiplyInPlace(matrix2);
assertEquals(matrix1, matrix3);
}

@Test
public void testElementDivide() {
double[][] data1 = new double[ROWS][COLS];
double[][] data2 = new double[ROWS][COLS];
Utilities.populate(data1, random);
Utilities.populate(data2, random);
Matrix matrix1 = new Matrix(data1, false);
Matrix matrix2 = new Matrix(data1, false);
Matrix matrix3 = matrix1.elementDivide(matrix2);
matrix1.elementDivideInPlace(matrix2);
assertEquals(matrix1, matrix3);
}

@Test
public void testMultiplication() {
double[][] data1 = new double[ROWS][ROWS];
double[][] data2 = new double[ROWS][ROWS];
Utilities.populate(data1, random);
Utilities.populate(data2, random);
Matrix matrix1 = new Matrix(data1, false);
Matrix matrix2 = new Matrix(data1, false);
Matrix matrix3 = matrix1.multiply(matrix2);
matrix1.multiplyInPlace(matrix2);
assertEquals(matrix1, matrix3);
}

@Test
public void testScale() {
double[][] data1 = new double[ROWS][ROWS];
double scale = random.nextDouble();
Utilities.populate(data1, random);
Matrix matrix1 = new Matrix(data1, false);
Matrix matrix2 = matrix1.scale(scale);
matrix1.scaleInPlace(scale);
assertEquals(matrix1, matrix2);
}

@Test
public void testExponentiation() {
double[][] data = new double[ROWS][ROWS];
Utilities.populate(data, random);
int power = random.nextInt(POWER_BOUND);

Matrix matrix1 = new Matrix(data, false);
Matrix matrix2 = new Matrix(data, true);
Matrix product = Matrix.identity(ROWS);
for (int i = 0; i < power; i++) {
product.multiplyInPlace(matrix1);
}
matrix2 = matrix2.exponentiate(power);
assertEquals(product, matrix2);
}

@Test
public void testOrthogonal() {
for (int i = 0; i < COUNT; i++) {
Matrix matrix = Matrix.identity(ROWS);
assertTrue(matrix.isOrthogonal());
double theta = random.nextDouble() * Math.PI * 2;
double cos = Math.cos(theta);
double sin = Math.sin(theta);
matrix = Matrix.fromLinearArray(2, 2, cos, -sin, sin, cos);
assertTrue(matrix.isOrthogonal());
matrix = Matrix.fromLinearArray(2, 2, cos, sin, sin, -cos);
assertTrue(matrix.isOrthogonal());
}
}
}


# Requests

Besides all general aspects, I’d like some guidance on what to do next? Pointers on how to implement a method or algorithm in a more effective way will also be appreciated. Things to do to avoid numerical stability is also on my mind. Suggestions for more test cases are also welcome.

2. Building software – Apache Maven – Link
3. Testing framework – TestNG – Link

Note that 2 and 3 are well integrated into 1.

Get this bounty!!!

## Behavior:

With the following texture

I get a weird behavior when I try to rotate my particles in the direction they are moving. Essentially, the quad seems to rotate, but just barely so that when they bounce on the ground or when they spawn we can clearly see that they are not aligned upward when the particle is going upward.

## Shader code and process explanation:

particle_billboard.geom:

//This geometry shader takes points as a
//primitive input and return quads built with
//triangle strips

#version 430

layout (points) in;
layout (triangle_strip, max_vertices = 4) out;

layout (std140, binding = 0) uniform CameraInfo {
mat4  ProjectionView; // Projection * View
vec3  eye;    // camera position
};

uniform float particle_size;

in  vec4 ex_Color[];
in  vec3 ex_Direction[]; // Only one direction per particle so index is always direction[0]
out vec2 TexCoord;
out vec4 FragColor;

float angle_between_vector(vec3 v1, vec3 v2)
{
return acos(
dot(v1, v2)/ length(v1)*length(v2)
);
}

vec3 project_vector_on_plane(vec3 v, vec3 normalized_plane_normal)
{
return v - dot(v, normalized_plane_normal) * normalized_plane_normal;
}

{
vec4 qr;
qr.x = axis.x * sin(half_angle);
qr.y = axis.y * sin(half_angle);
qr.z = axis.z * sin(half_angle);
qr.w = cos(half_angle);
return qr;
}

vec3 rotate_vertex_position(vec3 position, vec3 axis, float rad_angle)
{
vec3 v = position.xyz;
return v + 2.0 * cross(q.xyz, cross(q.xyz, v) + q.w * v);
}

void main(void) {
// Built-in variable "gl_in" of the geometry shader
// gives us the index in the primitive. Points only
// have one index, so gl_in[0] it is.

vec3 wPos = gl_in[0].gl_Position.xyz;
vec3 wRotatedPos = wPos;
vec3 originalPos = wPos;
vec3 wPosToCamera = normalize(eye - wPos);
vec3 wUp = vec3(0.0, 1.0, 0.0);
vec3 wRight = cross(wPosToCamera, wUp);

FragColor = ex_Color[0]; // Points only have one vertex

// 0.1 - Project ex_Direction on

// 1- Compute angle between right and direction
float wDirectionAngle = angle_between_vector(project_vector_on_plane(ex_Direction[0], wPosToCamera), wRight);

// Lower left corner
wPos -= (wRight * particle_size/2);
wPos.y -= particle_size/2;
wRotatedPos = rotate_vertex_position(wPos - originalPos, wPosToCamera, wDirectionAngle) + originalPos;
gl_Position = ProjectionView * vec4(wRotatedPos, 1.0);
TexCoord = vec2(0.0, 0.0);
EmitVertex();

// Upper left corner
wPos.y += particle_size;
wRotatedPos = rotate_vertex_position(wPos - originalPos, wPosToCamera, wDirectionAngle) + originalPos;
gl_Position = ProjectionView * vec4(wRotatedPos, 1.0);
TexCoord = vec2(0.0, 1.0);
EmitVertex();

// lower right corner
wPos.y -= particle_size;
wPos += wRight * particle_size;
wRotatedPos = rotate_vertex_position(wPos - originalPos, wPosToCamera, wDirectionAngle) + originalPos;
gl_Position = ProjectionView * vec4(wRotatedPos, 1.0);
TexCoord = vec2(1.0, 0.0);
EmitVertex();

// Upper right corner
wPos.y += particle_size;
wRotatedPos = rotate_vertex_position(wPos - originalPos, wPosToCamera, wDirectionAngle) + originalPos;
gl_Position = ProjectionView * vec4(wRotatedPos, 1.0);
TexCoord = vec2(1.0, 1.0);
EmitVertex();

EndPrimitive();
}


That is pretty much it for the code, I don’t really do much else in the vertex and frag shaders… But I can provide them if necessary.

First of all, I draw the particles with quads that I build in the geometry shader from a single point. The quads are built with a right and up vector a little bit like this tutorial does, mainly for making the quads always face the camera.

Now, I visualize the texture as being aligned with the right vector, so I try to compute the angle (wDirectionAngle) between this vector and the direction vector of the particle projected on the camera plane. The camera plane being the plane with the normal in the same direction as the particle_position-to-camera-vectorwPosToCamera“. Finally, I use the wPosToCamera as a quaternion vector to rotate my generated quads coordinate around the center of the quad. The rotation around the quaternion and the parameters I use are also heavily inspired on this tutorial.

I’ve debugged the direction angle and it appears to be correct, although I’m not 100% sure. What could possibly cause this to happen, is there a flaw in the logical process? I couldn’t really get any inspiration of similar posts, since they are all about 2D.

Get this bounty!!!

## #StackBounty: #c# #mathematics #projectile-physics #trajectory #aiming Projectile Aim Prediction with Acceleration

### Bounty: 50

I’m trying to solve the classic shoot moving object problem but with acceleration attached to that changes it from a quadratic to quartic formula but my math skills are not this good sadly as i prefer to speak in code and not formulas.

i found this https://wiki.beyondunreal.com/Legacy:Projectile_Aiming. i ported it and it’s almost what i’m looking for but missing control over the acceleration as its only made for gravity or no gravity and even this i got working only with tricks and modifying it without knowing what i’m doing is getting me only that far.

I made myself a Prototype Interface for the minimum of what i’m trying to get out of it

public class PredictionResult {
public bool IsInRange; // can we even hit the target?
public Vector3[] ShotVelocity; // should be up to 4 possible values
public float[] ShotImpactTime; // how long each shot takes to arrive at the predicted target
public Vector3[] ShotImpactLocation; // somewhat redundant yet still useful if available
}
// only really the delta between start and target matter but its up to the function
// same for startVel and targetVel where startVel is the shooters speed that gets added to the bullet
// the bullet and target can have different accelerations (standing on ground => no gravity or bullet is not affect by gravity)
public PredictionResult ShootAtTarget(Vector3 start, Vector3 target, Vector3 startVel, Vector3 targetVel, Vector3 bulletAccel, Vector3 targetAccel, float bulletSpeed);


any help solving this would be great

Get this bounty!!!

# Introduction

I’m creating a game where the player can obtain 1 to 3 stars for each level based on the score it gets (based on the completion time).
The levels are grouped in “worlds” each of which is unlocked when the user obtains a given number of stars in the previous levels.

For example, world 1 has 5 levels and to unlock world 2 the user needs to gain at least 5 stars (thus at leas one star per level in average).

Here is a basic idea of the words -> n° of levels in that world -> stars to unlock

1 -> 5 -> 0 (of course)
2 -> 5 -> 5 / 15 (33%)
3 -> 5 -> 10 / 30 (33%)
4 -> 7 -> 25 / 45 (55%)
5 -> 7 -> 30 / 66 (45%)
...


To determine the score you should reach, for each level, to get 0 / 1 / 2 / 3 stars, I recorded some play stats from a bunch of beta-testers obtaining a normal distribution of play times for each level.
Given each distribution, I should be able to answer this question for each level:

at what score should I reward n stars in order for x% of the
player to get n stars overall?

# Tuning

So now I can change the thresholds for each star at each level (or group of levels) in order to filter the percentage of players that will obtain a given number of stars at some point.
This way I can set the game difficulty as the difficulty to unlock a given world which is the percentage of people who are good enough to gain enough stars to unlock that world.

For example for early world progress difficulty, I chose a base of 1.22 as exponential base of the percentage reduction (difficulty growth), like this:

world -> difficulty coeff. -> perc. players

1 -> 0 -> 100%
2 -> 1.22 -> 98.78%
3 -> 1.4884 -> 98.5116%
4 -> 1.815848 -> 98.184152%
5 -> 2.21533456 -> 97.78466544%
...


This way, for example, the last world should be reached by about 47% of the players.

Now I want to know how to tune percentages of people gaining 1 / 2 / 3 stars in order to stick to this given percentage progression.
In order to do this, I found the minimal configurations of possible stars obtained in each world level in order to unlock the next one, for example:

world -> n° 1 stars -> n° 2 stars -> n° 3 stars -> stars to unlock next world
1 -> 5 -> 0 -> 0 -> 5
2 -> 10 -> 0 -> 0 -> 10
3 -> 5 -> 10 -> 0 -> 25
4 -> 12 -> 10 -> 0 -> 30
5 -> 18 -> 11 -> 0 -> 40
...


Note that if the user got 10 x 2 stars in the previous world, then it still has those 10 x 2 stars in later worlds, unless it tops them wit 3 stars.

# My Calculations

Now, for example, if I want 98.78% of the players to be able to unlock the second world, given they have to obtain minimum 1 star at each previous level, then p^5 = 0.9878 and p = rad(5, 0.9878) ≈ 0.9975, so 99.75% of the players should be able to get at leas 1 star in each level of the first world.

For the third world, things get a little harder, as the 1 star probabilities for the first 5 levels are locked now.
Players must be able to obtain at least 1 star in each of the 10 levels of the first two worlds with an overall probability of 98.5116%, but the probability to obtain 1 star in the first 5 levels is locked at 99.75%, with an overall probability of obtaining 1 star in each of the first 5 levels of 98.78%.
So I had to solve this equation: p * 0.9878 = 0.985116 so p = 0.985116 / 0.9878 = 0.9973 which is the probability p to get at leas 1 star in all the 5 levels of the second world.
So the probability to obtain 1 star for each single level of the second world is p = rad(5, 0.9973) = 0.9995 which is slightly higher than the previous world.

Fourth world gets even weirder, as the restrictions shifts on the 10 x 2 stars that players need to obtain in the 15 previous levels. To do this, I used binomial distribution to find the probability to extract at least 10 successes over 15 attempts which gave a probability of 58.79% to obtain 2 stars in each of the 15 levels of the first 3 worlds in order to have a 98.184% probability to finish the third world with at leas 10 x 2 stars.

Fifth world differs from fourth only by 7 x 1 star so I just calculated p * 0.98184 = 0.97784 where p is the probability to obtain 1 star in each of the 7 levels of the fifth world, obtaining a probability to obtain 1 star for each single level of the fifth world of 99.50%.

# Problems

Now I’m stuck at world 6 where the user is required to obtain at least 11 x 2 stars. How do I calculate this probability? I can use the binomial distribution, but the probabilities of the events are not the same everywhere as the probability to obtain 2 stars in the first 5 words is locked.

Is there any formula to help me with this?
Is there any simpler / more direct approach I can fallow?
Does any of this make any sense at all?

Get this bounty!!!

## #StackBounty: #books #exercises #mathematics Math book: how to write Exercise and Answers

### Bounty: 50

EDIT
within documentclass[12pt]{book}I want to create chapter-wise exercises and put all the solutions (with or without hints) at the end of the book. The answer should include page No of the exercise as given in the attached jpg file.

I want to do this in simple and non-tedious way like: For the input of the questions, I just want to add questionfor each question and similar for answer, but all the answer should come at the end of the book.

Exercise style:

Solution stype:

Get this bounty!!!

## #HackerRank: Computing the Correlation

### Problem

You are given the scores of N students in three different subjects – MathematicsPhysics and Chemistry; all of which have been graded on a scale of 0 to 100. Your task is to compute the Pearson product-moment correlation coefficient between the scores of different pairs of subjects (Mathematics and Physics, Physics and Chemistry, Mathematics and Chemistry) based on this data. This data is based on the records of the CBSE K-12 Examination – a national school leaving examination in India, for the year 2013.

Pearson product-moment correlation coefficient

This is a measure of linear correlation described well on this Wikipedia page. The formula, in brief, is given by:

where x and y denote the two vectors between which the correlation is to be measured.

Input Format

The first row contains an integer N.
This is followed by N rows containing three tab-space (‘\t’) separated integers, M P C corresponding to a candidate’s scores in Mathematics, Physics and Chemistry respectively.
Each row corresponds to the scores attained by a unique candidate in these three subjects.

Input Constraints

1 <= N <= 5 x 105
0 <= M, P, C <= 100

Output Format

The output should contain three lines, with correlation coefficients computed
and rounded off correct to exactly 2 decimal places.
The first line should contain the correlation coefficient between Mathematics and Physics scores.
The second line should contain the correlation coefficient between Physics and Chemistry scores.
The third line should contain the correlation coefficient between Chemistry and Mathematics scores.

So, your output should look like this (these values are only for explanatory purposes):

0.12
0.13
0.95


Test Cases

There is one sample test case with scores obtained in Mathematics, Physics and Chemistry by 20 students. The hidden test case contains the scores obtained by all the candidates who appeared for the examination and took all three tests (Mathematics, Physics and Chemistry).
Think: How can you efficiently compute the correlation coefficients within the given time constraints, while handling the scores of nearly 400k students?

Sample Input

20
73  72  76
48  67  76
95  92  95
95  95  96
33  59  79
47  58  74
98  95  97
91  94  97
95  84  90
93  83  90
70  70  78
85  79  91
33  67  76
47  73  90
95  87  95
84  86  95
43  63  75
95  92  100
54  80  87
72  76  90


Sample Output

0.89
0.92
0.81


There is no special library support available for this challenge.

## What is the difference between linear regression on y with x and x with y?

The Pearson correlation coefficient of x and y is the same, whether you compute pearson(x, y) or pearson(y, x). This suggests that doing a linear regression of y given x or x given y should be the same, but that’s the case.

The best way to think about this is to imagine a scatter plot of points with y on the vertical axis and x represented by the horizontal axis. Given this framework, you see a cloud of points, which may be vaguely circular, or may be elongated into an ellipse. What you are trying to do in regression is find what might be called the ‘line of best fit’. However, while this seems straightforward, we need to figure out what we mean by ‘best’, and that means we must define what it would be for a line to be good, or for one line to be better than another, etc. Specifically, we must stipulate a loss function. A loss function gives us a way to say how ‘bad’ something is, and thus, when we minimize that, we make our line as ‘good’ as possible, or find the ‘best’ line.

Traditionally, when we conduct a regression analysis, we find estimates of the slope and intercept so as to minimize the sum of squared errors. These are defined as follows:

In terms of our scatter plot, this means we are minimizing the sum of the vertical distances between the observed data points and the line.

On the other hand, it is perfectly reasonable to regress x onto y, but in that case, we would put x on the vertical axis, and so on. If we kept our plot as is (with x on the horizontal axis), regressing x onto y (again, using a slightly adapted version of the above equation with x and y switched) means that we would be minimizing the sum of the horizontal distances between the observed data points and the line. This sounds very similar, but is not quite the same thing. (The way to recognize this is to do it both ways, and then algebraically convert one set of parameter estimates into the terms of the other. Comparing the first model with the rearranged version of the second model, it becomes easy to see that they are not the same.)

Note that neither way would produce the same line we would intuitively draw if someone handed us a piece of graph paper with points plotted on it. In that case, we would draw a line straight through the center, but minimizing the vertical distance yields a line that is slightly flatter (i.e., with a shallower slope), whereas minimizing the horizontal distance yields a line that is slightly steeper.

A correlation is symmetrical x is as correlated with y as y is with x. The Pearson product-moment correlation can be understood within a regression context, however. The correlation coefficient, r, is the slope of the regression line when both variables have been standardized first. That is, you first subtracted off the mean from each observation, and then divided the differences by the standard deviation. The cloud of data points will now be centered on the origin, and the slope would be the same whether you regressed y onto x, or x onto y.

Now, why does this matter? Using our traditional loss function, we are saying that all of the error is in only one of the variables (viz., y). That is, we are saying that x is measured without error and constitutes the set of values we care about, but that y has sampling error. This is very different from saying the converse. This was important in an interesting historical episode: In the late 70’s and early 80’s in the US, the case was made that there was discrimination against women in the workplace, and this was backed up with regression analyses showing that women with equal backgrounds (e.g., qualifications, experience, etc.) were paid, on average, less than men. Critics (or just people who were extra thorough) reasoned that if this was true, women who were paid equally with men would have to be more highly qualified, but when this was checked, it was found that although the results were ‘significant’ when assessed the one way, they were not ‘significant’ when checked the other way, which threw everyone involved into a tizzy. See here for a famous paper that tried to clear the issue up.