jidt/java/unittests/infodynamics/utils/MatrixUtilsTest.java

571 lines
22 KiB
Java
Executable File

/*
* Java Information Dynamics Toolkit (JIDT)
* Copyright (C) 2012, Joseph T. Lizier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package infodynamics.utils;
import java.util.List;
import junit.framework.TestCase;
/**
* Test functionality of MatrixUtils methods.
*
* @author Joseph Lizier joseph.lizier_at_gmail.com
*
*/
public class MatrixUtilsTest extends TestCase {
private static double OCTAVE_RESOLUTION = 0.00001;
public void testIdentityMatrix() {
// Test identity matrix generation, including for size 0
for (int n = 0; n < 10; n++) {
double[][] I = MatrixUtils.identityMatrix(n);
assertNotNull(I);
assertEquals(n, I.length);
for (int i = 0; i < n; i++) {
assertEquals(n, I[i].length);
for (int j = 0; j < n; j++) {
assertEquals(i == j ? 1 : 0, I[i][j], 0.000000001);
}
}
}
}
public void testCovariance() throws Exception {
// Load a test data file which contains a time series and
// covariance matrix as computed by matlab:
// We'll just take the first two columns from this data set
OctaveFileReader ofr = new OctaveFileReader(
"demos/data/Network-GaussianLinear-N100-T100-p0.04-b0.50-c0.50-dir-disc-repeat1.txt");
double[][] data = ofr.getDouble2DMatrix("timeseries");
double[][] expectedCovariance = ofr.getDouble2DMatrix("empiricalCovariance");
double[][] computedCovariance = MatrixUtils.covarianceMatrix(data);
checkMatrix(expectedCovariance, computedCovariance, 0.00001);
// And test that it's still correct if we supply the data in 2 separate
// parts:
double[][] part1 = MatrixUtils.selectColumns(data,
MatrixUtils.range(0, 49));
double[][] part2 = MatrixUtils.selectColumns(data,
MatrixUtils.range(50, 99));
double[][] split2ComputedCovariance = MatrixUtils.covarianceMatrix(
part1, part2);
checkMatrix(expectedCovariance, split2ComputedCovariance, 0.00001);
// And test that it's still correct if we supply the data in 3 separate
// parts:
double[][] part2a = MatrixUtils.selectColumns(data,
MatrixUtils.range(50, 74));
double[][] part2b = MatrixUtils.selectColumns(data,
MatrixUtils.range(75, 99));
double[][] split3ComputedCovariance = MatrixUtils.covarianceMatrix(
part1, part2a, part2b);
checkMatrix(expectedCovariance, split3ComputedCovariance, 0.00001);
}
/**
* Test our Cholesky decomposition implementation
*
* @throws Exception
*/
public void testCholesky() throws Exception {
// Check some ordinary Cholesky decompositions:
double[][] A = {{6, 2, 3}, {2, 5, 1}, {3, 1, 4}};
// Expected result from Octave:
double[][] expectedL = {{2.44949, 0, 0}, {0.81650, 2.08167, 0},
{1.22474, 0, 1.58114}};
double[][] L = MatrixUtils.CholeskyDecomposition(A);
checkMatrix(expectedL, L, OCTAVE_RESOLUTION);
double[][] A2 = {{6, 2, 3, 1}, {2, 5, 1, 0.5}, {3, 1, 4, 2}, {1, 0.5, 2, 3}};
// Expected result from Octave:
double[][] expectedL2 = {{2.44949, 0, 0, 0}, {0.81650, 2.08167, 0, 0},
{1.22474, 0, 1.58114, 0}, {0.40825, 0.08006, 0.94868, 1.38814}};
double[][] L2 = MatrixUtils.CholeskyDecomposition(A2);
checkMatrix(expectedL2, L2, OCTAVE_RESOLUTION);
// Now check that it picks up asymmetric A:
double[][] asymmetricA = {{6, 2, 3}, {2, 5, 1}, {3, 1.0001, 4}};
boolean flaggedException = false;
try {
MatrixUtils.CholeskyDecomposition(asymmetricA);
} catch (Exception e) {
flaggedException = true;
}
assertTrue(flaggedException);
}
public void testCholeskyNotPositiveDefinite() throws Exception {
// Now check that it picks up if A is not positive definite:
// Variable 2 is a scaled version of variable 1 here
double[][] notpositiveDefiniteA = {{1, 2}, {2, 4}};
boolean flaggedException = false;
int problematicRow = -1;
try {
MatrixUtils.CholeskyDecomposition(notpositiveDefiniteA);
} catch (NonPositiveDefiniteMatrixException npdme) {
flaggedException = true;
problematicRow = npdme.problematicRow;
}
assertTrue(flaggedException);
assertEquals(1, problematicRow);
// Test that we pick up a variable which duplicates another
double[][] copiedVariableA = {{1, 1.9, 1}, {1.9, 4, 1.9}, {1, 1.9, 1}};
problematicRow = -1;
flaggedException = false;
try {
MatrixUtils.CholeskyDecomposition(copiedVariableA);
} catch (NonPositiveDefiniteMatrixException npdme) {
flaggedException = true;
problematicRow = npdme.problematicRow;
}
assertTrue(flaggedException);
assertEquals(2, problematicRow);
double[][] copiedVariableRow1A = {{1, 1, 1.9}, {1, 1, 1.9}, {1.9, 1.9, 4}};
problematicRow = -1;
flaggedException = false;
try {
MatrixUtils.CholeskyDecomposition(copiedVariableRow1A);
} catch (NonPositiveDefiniteMatrixException npdme) {
flaggedException = true;
problematicRow = npdme.problematicRow;
}
assertTrue(flaggedException);
assertEquals(1, problematicRow);
// Now do this with some real data:
RandomGenerator rg = new RandomGenerator();
int N = 1000;
double[][] data = new double[N][3];
MatrixUtils.copyIntoColumn(data, 0, rg.generateNormalData(N, 0, 1));
MatrixUtils.copyIntoColumn(data, 1, rg.generateNormalData(N, 0, 1));
MatrixUtils.copyIntoColumn(data, 2, rg.generateNormalData(N, 0, 1));
double[][] covarianceMatrix = MatrixUtils.covarianceMatrix(data);
// This one should be fine and not throw an exception:
MatrixUtils.CholeskyDecomposition(covarianceMatrix);
// But if we substitute in a linearly redundant variable:
MatrixUtils.copyIntoColumn(data, 2,
MatrixUtils.add(MatrixUtils.selectColumn(data, 0),
MatrixUtils.selectColumn(data, 1)));
covarianceMatrix = MatrixUtils.covarianceMatrix(data);
problematicRow = -1;
flaggedException = false;
double[][] L;
try {
// This should throw an exception:
L = MatrixUtils.CholeskyDecomposition(covarianceMatrix);
// Debug prints in case it does not (because tolerance in the CholeskyDecomposition was too large):
System.out.println(L[2][2]);
System.out.println("Det " + MatrixUtils.determinantViaCholeskyResult(L) + " > 0, unexpectedly, but it happens - precision is approx: " +
java.lang.Math.ulp(0));
MatrixUtils.printMatrix(System.out, covarianceMatrix);
System.out.printf("%.6f, %.6f\n", covarianceMatrix[0][2] - covarianceMatrix[0][1] - covarianceMatrix[0][0],
covarianceMatrix[1][2] - covarianceMatrix[1][1] - covarianceMatrix[1][0]);
} catch (NonPositiveDefiniteMatrixException npdme) {
flaggedException = true;
problematicRow = npdme.problematicRow;
L = null;
}
assertTrue(flaggedException);
assertEquals(2, problematicRow);
// Now if we swap the second and third variables in the data, it should still complain about the last one
// since that's the one it saw last:
double[] temp = MatrixUtils.selectColumn(data, 2);
MatrixUtils.copyIntoColumn(data, 2, MatrixUtils.selectColumn(data, 1));
MatrixUtils.copyIntoColumn(data, 1, temp);
covarianceMatrix = MatrixUtils.covarianceMatrix(data);
problematicRow = -1;
flaggedException = false;
try {
// This should throw an exception:
L = MatrixUtils.CholeskyDecomposition(covarianceMatrix);
} catch (NonPositiveDefiniteMatrixException npdme) {
flaggedException = true;
problematicRow = npdme.problematicRow;
L = null;
}
assertTrue(flaggedException);
assertEquals(2, problematicRow);
}
/**
* Test the inversion of symmetric positive definite matrices
*
* @throws Exception
*/
public void testInverseOfSymmPosDefMatrices() throws Exception {
// Check some ordinary matrices:
double[][] A = {{6, 2, 3}, {2, 5, 1}, {3, 1, 4}};
// Expected result from Octave:
double[][] expectedInv = {{0.29231, -0.07692, -0.2},
{-0.07692, 0.23077, 0}, {-0.2, 0, 0.4}};
double[][] inv = MatrixUtils.invertSymmPosDefMatrix(A);
checkMatrix(expectedInv, inv, OCTAVE_RESOLUTION);
double[][] A2 = {{6, 2, 3, 1}, {2, 5, 1, 0.5}, {3, 1, 4, 2}, {1, 0.5, 2, 3}};
// Expected result from Octave:
double[][] expectedInv2 = {{0.303393, -0.079840, -0.245509, 0.075848},
{-0.079840, 0.231537, 0.011976, -0.019960},
{-0.245509, 0.011976, 0.586826, -0.311377},
{0.075848, -0.019960, -0.311377, 0.518962}};
double[][] inv2 = MatrixUtils.invertSymmPosDefMatrix(A2);
checkMatrix(expectedInv2, inv2, OCTAVE_RESOLUTION);
// Now check that it picks up asymmetric A:
double[][] asymmetricA = {{6, 2, 3}, {2, 5, 1}, {3, 1.0001, 4}};
boolean flaggedException = false;
try {
MatrixUtils.invertSymmPosDefMatrix(asymmetricA);
} catch (Exception e) {
flaggedException = true;
}
assertTrue(flaggedException);
// Now check that it picks up if A is not positive definite:
double[][] notpositiveDefiniteA = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
flaggedException = false;
try {
MatrixUtils.invertSymmPosDefMatrix(notpositiveDefiniteA);
} catch (Exception e) {
flaggedException = true;
}
assertTrue(flaggedException);
}
/**
* Test the solving of matrix equations via Cholesky decomposition.
* Solving A*X = B
*
* @throws Exception
*/
public void testSolvingMatrixEquationsOfSymmPosDefMatrices() throws Exception {
// Check some ordinary matrices:
double[][] A = {{6, 2, 3}, {2, 5, 1}, {3, 1, 4}};
double[][] B = {{5}, {4}, {3}};
// Expected result from Octave:
double[][] expectedX = {{0.55385}, {0.53846}, {0.20000}};
double[][] X = MatrixUtils.solveViaCholeskyResult(
MatrixUtils.CholeskyDecomposition(A), B);
checkMatrix(expectedX, X, OCTAVE_RESOLUTION);
// Check more complicated example
double[][] A2 = {{6, 2, 3, 1}, {2, 5, 1, 0.5}, {3, 1, 4, 2}, {1, 0.5, 2, 3}};
double[][] B2 = {{10, 5, 4, 12}, {4, 6, -1, 4.3}, {20, 1, 0, -5}, {6, 3, 2, 1}};
double[][] expectedX2 = {{-1.740519, 1.019960, 1.445110, 4.600798},
{0.247505, 0.942116, -0.590818, -0.042315},
{7.461078, -1.502994, -1.616766, -6.140120},
{-2.435130, 1.504990, 1.361277, 2.900200}};
double[][] X2 = MatrixUtils.solveViaCholeskyResult(
MatrixUtils.CholeskyDecomposition(A2), B2);
checkMatrix(expectedX2, X2, OCTAVE_RESOLUTION);
// TODO Check error conditions
}
public void testDeterminant() throws Exception {
// test some error conditions:
double[][] AnonSquare = {{6, 2}, {2, 5, 1}, {3, 1, 4}};
boolean flaggedException = false;
try {
MatrixUtils.determinant(AnonSquare);
} catch (Exception e) {
flaggedException = true;
}
assertTrue(flaggedException);
// Test some simple examples:
double[][] A1 = {{3.445454}};
assertEquals(3.445454, MatrixUtils.determinant(A1), OCTAVE_RESOLUTION);
double[][] A2 = {{6, 2}, {2, 5}};
assertEquals(26, MatrixUtils.determinant(A2), OCTAVE_RESOLUTION);
// Check against value computed by Octave:
double[][] A = {{6, 2, 3}, {2, 5, 1}, {3, 1, 4}};
assertEquals(65, MatrixUtils.determinant(A), OCTAVE_RESOLUTION);
// Check zero determinant case
double[][] AzeroDet = {{6, 2, 3}, {2, 5, 1}, {10, -14, 5}};
assertEquals(0, MatrixUtils.determinant(AzeroDet), OCTAVE_RESOLUTION);
}
public void testDeterminantSymmPosDef() throws Exception {
// Test some simple examples:
double[][] A1 = {{3.445454}};
assertEquals(3.445454,
MatrixUtils.determinantSymmPosDefMatrix(A1), OCTAVE_RESOLUTION);
double[][] A2 = {{6, 2}, {2, 5}};
assertEquals(26, MatrixUtils.determinantSymmPosDefMatrix(A2), OCTAVE_RESOLUTION);
// Check against value computed by Octave:
double[][] A = {{6, 2, 3}, {2, 5, 1}, {3, 1, 4}};
assertEquals(65, MatrixUtils.determinantSymmPosDefMatrix(A), OCTAVE_RESOLUTION);
// Now check that it picks up asymmetric A:
double[][] asymmetricA = {{6, 2, 3}, {2, 5, 1}, {3, 1.0001, 4}};
boolean flaggedException = false;
try {
MatrixUtils.determinantSymmPosDefMatrix(asymmetricA);
} catch (Exception e) {
flaggedException = true;
}
assertTrue(flaggedException);
// Now check that it picks up if A is not positive definite:
double[][] notpositiveDefiniteA = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
flaggedException = false;
try {
MatrixUtils.determinantSymmPosDefMatrix(notpositiveDefiniteA);
} catch (Exception e) {
flaggedException = true;
}
assertTrue(flaggedException);
}
public void testSortIndices() {
double[] array1 = {0.1, 0.2, 0.3, 0.4, 0.5};
checkArray(new int[] {0, 1, 2, 3, 4}, MatrixUtils.sortIndices(array1));
double[] array2 = {0.5, 0.4, 0.3, 0.2, 0.1};
checkArray(new int[] {4, 3, 2, 1, 0}, MatrixUtils.sortIndices(array2));
double[] array3 = {0.3, 0.1, 0.5, 0.4, 0.2};
checkArray(new int[] {1, 4, 0, 3, 2}, MatrixUtils.sortIndices(array3));
}
public void testDelayEmbeddings() throws Exception {
double[] array1 = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
// Do a standard delay embedding with tau 1
checkMatrix(new double[][] { {4, 3, 2, 1, 0}, {5, 4, 3, 2, 1},
{6, 5, 4, 3, 2}, {7, 6, 5, 4, 3},
{8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} },
MatrixUtils.makeDelayEmbeddingVector(array1, 5, 4, 6),
0.00001);
// Now specify tau explicitly
checkMatrix(new double[][] { {4, 3, 2, 1, 0}, {5, 4, 3, 2, 1},
{6, 5, 4, 3, 2}, {7, 6, 5, 4, 3},
{8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} },
MatrixUtils.makeDelayEmbeddingVector(array1, 5, 1, 4, 6),
0.00001);
// Do same standard delay embedding but starting at an offset
checkMatrix(new double[][] { {8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} },
MatrixUtils.makeDelayEmbeddingVector(array1, 5, 8, 2),
0.00001);
// Now specify tau explicitly
checkMatrix(new double[][] { {8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} },
MatrixUtils.makeDelayEmbeddingVector(array1, 5, 1, 8, 2),
0.00001);
// Try with tau 2
checkMatrix(new double[][] { {8, 6, 4, 2, 0}, {9, 7, 5, 3, 1} },
MatrixUtils.makeDelayEmbeddingVector(array1, 5, 2, 8, 2),
0.00001);
// Try with tau 3
checkMatrix(new double[][] { {6, 3, 0}, {7, 4, 1}, {8, 5, 2}, {9, 6, 3} },
MatrixUtils.makeDelayEmbeddingVector(array1, 3, 3, 6, 4),
0.00001);
}
/**
* Check that all entries in the given matrix match those of the expected
* matrix
*
* @param expected
* @param actual
* @param resolution
*/
public static void checkMatrix(double[][] expected, double[][] actual, double resolution) {
for (int r = 0; r < expected.length; r++) {
for (int c = 0; c < expected[r].length; c++) {
assertEquals(expected[r][c], actual[r][c], resolution);
}
}
}
/**
* Check that all entries in the given array match those of the expected
* array
*
* @param expected
* @param actual
*/
public static void checkArray(int[] expected, int[] actual) {
for (int r = 0; r < expected.length; r++) {
assertEquals(expected[r], actual[r]);
}
}
public void test2DArrayCopy() {
double[][] temp = {{1,2,3,4}, {4,5,6,7}, {7,8,9,10}, {10,11,12,13}};
double[][] newMatrix = new double[10][10];
MatrixUtils.arrayCopy(temp, 1, 1, newMatrix, 3, 3, 3, 3);
for (int r = 0; r < 3; r++) {
for (int c = 0; c < 3; c++) {
assertEquals(temp[r+1][c+1], newMatrix[r+3][c+3]);
}
}
}
public void testAddAndRemoveFromMean() {
RandomGenerator rg = new RandomGenerator();
int N = 1000;
int cols = 3;
int rowsToRemove = 10;
double[][] data = rg.generateNormalData(N, cols, 0, 1);
for (int c = 0; c < cols; c++) {
double mean = MatrixUtils.mean(data, c);
// Remove mean from row c as well
double meanAfterRemoval = MatrixUtils.removeFromColumnMean(
data, c, c, rowsToRemove, mean, N);
// And now try manually:
double sum = mean * (double) N;
for (int r = c; r < c + rowsToRemove; r++) {
sum -= data[r][c];
}
double meanAfterRemovalManual = sum / (double) (N - rowsToRemove);
assertEquals(meanAfterRemovalManual, meanAfterRemoval, 1e-5);
// Then add back in:
double meanAfterAddingBackIn = MatrixUtils.addToColumnMean(
data, c, c, rowsToRemove, meanAfterRemoval, N - rowsToRemove);
assertEquals(mean, meanAfterAddingBackIn, 1e-5);
// Take one row out then try swapping it with another:
double meanAfterOneRemoval = MatrixUtils.removeFromColumnMean(
data, c, c, 1, mean, N);
double meanAfterSwappedRemoval = MatrixUtils.swapIntoColumnMean(
data, c, c, c+1, meanAfterOneRemoval, N-1);
double meanAfterSwapManual = (mean * (double) N - data[c+1][c]) / ((double) (N - 1));
assertEquals(meanAfterSwapManual, meanAfterSwappedRemoval, 1e-5);
}
}
public void testAddAndRemoveFromCovariance() {
RandomGenerator rg = new RandomGenerator();
// Run a short test:
double[][] shortData = new double[10][1];
shortData[0][0] = 1; // Only item that contributes to the variance (otherwise it will be zero)
double[][] shortCovariance = MatrixUtils.covarianceMatrix(shortData);
double initialCovariance = shortCovariance[0][0];
System.out.printf("Initial covariance = %.5f\n", initialCovariance);
MatrixUtils.removeFromCovarianceMatrix(shortCovariance,
shortData, null, null,
new int[] {0}, new int[] {}, new int[] {},
new double[] {0.1}, new double[] {}, new double[] {},
0, 1, 10);
// and now check what happens if we remove this row manually:
double[][] cutShortData = MatrixUtils.selectRows(shortData, 1, 9);
double[][] covShortAfterManualCut = MatrixUtils.covarianceMatrix(cutShortData);
double covMatRemoval = MatrixUtils.computeCovarianceMatrixRemoval(shortData, shortData,
0, 0, 0.1, 0.1, 10, 0, 1);
System.out.printf("computeCovarianceMatrixRemoval returns %.5f\n", covMatRemoval);
assertEquals(covShortAfterManualCut[0][0], shortCovariance[0][0], 1e-6);
assertEquals(0, shortCovariance[0][0], 1e-6);
// or indeed if and now check what happens if we add this row back in:
MatrixUtils.addToCovarianceMatrix(shortCovariance,
shortData, null, null,
new int[] {0}, new int[] {}, new int[] {},
new double[] {0}, new double[] {}, new double[] {},
0, 1, 9);
assertEquals(initialCovariance, shortCovariance[0][0], 1e-6);
int N = 100;
int cols = 3;
int rowsToRemove = 10;
double[][] data = rg.generateNormalData(N, cols, 0, 1);
double[][] covariances = MatrixUtils.covarianceMatrix(data);
double[][] covariancesOriginal = MatrixUtils.arrayCopy(covariances);
MatrixUtils.removeFromCovarianceMatrix(covariances,
MatrixUtils.selectColumns(data, new int[] {0}), MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {2}),
new int[] {0}, new int[] {0}, new int[] {0},
new double[] {MatrixUtils.mean(data, 0)}, new double[] {MatrixUtils.mean(data, 1)}, new double[] {MatrixUtils.mean(data, 2)},
0, rowsToRemove, N);
// and now check what happens if we remove these rows manually:
double[][] cutData = MatrixUtils.selectRows(data, rowsToRemove, N-rowsToRemove);
double[][] covAfterManualCut = MatrixUtils.covarianceMatrix(cutData);
// Need fairly large tolerance here, because the two are being computed in different
/// ways and so have slightly different numerical errors
checkMatrix(covAfterManualCut, covariances, 1e-6);
// And check that all is resolved once we insert them again
double runningMean0 = MatrixUtils.mean(data, 0, rowsToRemove, N-rowsToRemove);
double runningMean1 = MatrixUtils.mean(data, 1, rowsToRemove, N-rowsToRemove);
double runningMean2 = MatrixUtils.mean(data, 2, rowsToRemove, N-rowsToRemove);
MatrixUtils.addToCovarianceMatrix(covariances,
MatrixUtils.selectColumns(data, new int[] {0}), MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {2}),
new int[] {0}, new int[] {0}, new int[] {0},
new double[] {runningMean0},
new double[] {runningMean1},
new double[] {runningMean2},
0, rowsToRemove, N - rowsToRemove);
// Maybe need larger tolerance here, because the two are being computed in different
/// ways and so have slightly different numerical errors
checkMatrix(covariancesOriginal, covariances, 1e-6);
// Now need to test that we can swap rows in and out ok with new code
List<Integer> rows = MatrixUtils.createArrayList(MatrixUtils.range(rowsToRemove, N-1));
for (int r = rowsToRemove; r < rowsToRemove + 10; r++) {
// Swap row r - 1 in, and row r out:
rows.remove(0); // row r will be the first index for this test
rows.add(r-1);
// First from covariance, since it uses the old means:
MatrixUtils.swapIntoCovarianceMatrix(
covAfterManualCut,
MatrixUtils.selectColumns(data, new int[] {0}), MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {2}),
new int[] {0}, new int[] {0}, new int[] {0},
new double[] {runningMean0}, new double[] {runningMean1}, new double[] {runningMean2},
r - 1, r, N-rowsToRemove);
// Then from the means:
runningMean0 = MatrixUtils.swapIntoColumnMean(data, 0, r-1, r, runningMean0, N-rowsToRemove);
runningMean1 = MatrixUtils.swapIntoColumnMean(data, 1, r-1, r, runningMean1, N-rowsToRemove);
runningMean2 = MatrixUtils.swapIntoColumnMean(data, 2, r-1, r, runningMean2, N-rowsToRemove);
// Now check them:
double[][] extractedData = MatrixUtils.selectRows(data, rows);
double[] means = MatrixUtils.means(extractedData);
assertEquals(means[0], runningMean0, 1e-6);
assertEquals(means[1], runningMean1, 1e-6);
assertEquals(means[2], runningMean2, 1e-6);
double[][] covarianceMatrix = MatrixUtils.covarianceMatrix(extractedData, means);
checkMatrix(covarianceMatrix, covAfterManualCut, 1e-6);
}
}
}