Menu
▾
▴
[cilib-devel] [PATCH 4/4] Permutations and LUDecomposition for determinant.
|
From: Gary P. <gpa...@gm...> - 2009-08-13 06:46:21
|
Added the option of obtaining a list of permutations for a list and to
create a LU decomposition matrix in order to get the determinant of a
n x n matrix class.
Signed-off-by: Gary Pampara <gpa...@gm...>
---
.../java/net/sourceforge/cilib/math/Maths.java | 104 +++++++++++++++
.../problem/dataset/MatrixDataSetBuilder.java | 2 +-
.../problem/mappingproblem/MappingProblem.java | 2 +-
.../cilib/type/types/container/Matrix.java | 134 +++++++++++++++++---
.../java/net/sourceforge/cilib/math/MathsTest.java | 67 ++++++----
.../cilib/type/types/container/MatrixTest.java | 65 +++++-----
6 files changed, 299 insertions(+), 75 deletions(-)
diff --git a/src/main/java/net/sourceforge/cilib/math/Maths.java b/src/main/java/net/sourceforge/cilib/math/Maths.java
index fce85e5..db2e489 100644
--- a/src/main/java/net/sourceforge/cilib/math/Maths.java
+++ b/src/main/java/net/sourceforge/cilib/math/Maths.java
@@ -21,6 +21,9 @@
*/
package net.sourceforge.cilib.math;
+import java.util.ArrayList;
+import java.util.Iterator;
+import java.util.List;
import net.sourceforge.cilib.math.random.generator.MersenneTwister;
import net.sourceforge.cilib.math.random.generator.Random;
@@ -91,6 +94,107 @@ public final class Maths {
return factorial(n) / factorial(n-r);
}
+ public static <T> Iterator<List<T>> permutation(final List<T> input, final int number) {
+ return new Iterator<List<T>>() {
+ private List<T> internalList = new ArrayList<T>(input); // Keep our own copy
+ private int n = input.size();
+ private int m = number;
+ private int[] index = initialize();
+ private boolean hasMore = true;
+
+ @Override
+ public boolean hasNext() {
+ return this.hasMore;
+ }
+
+ @Override
+ public List<T> next() {
+ if (!this.hasMore){
+ return null;
+ }
+ List<T> list = new ArrayList<T>(this.m);
+ for (int i = 0; i < this.m; i++) {
+ int thisIndexI = this.index[i];
+ T element = internalList.get(thisIndexI);
+ list.add(element);
+ }
+ moveIndex();
+ return list;
+ }
+
+ @Override
+ public void remove() {
+ throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ private int[] initialize() {
+ if (!(this.n >= m && m >= 0))
+ throw new IllegalStateException("Permutation error! n >= m");
+
+ int[] tmp = new int[this.n];
+ for (int i = 0; i < this.n; i++) {
+ tmp[i] = i;
+ }
+
+ reverseAfter(tmp, m - 1);
+ return tmp;
+ }
+
+ /**
+ * Reverse the index elements to the right of the specified index.
+ */
+ private void reverseAfter(int[] indicies, int i) {
+ int start = i + 1;
+ int end = this.n - 1;
+ while (start < end) {
+ int t = indicies[start];
+ indicies[start] = indicies[end];
+ indicies[end] = t;
+ start++;
+ end--;
+ }
+ }
+
+ private void moveIndex(){
+ // find the index of the first element that dips
+ int i = rightmostDip();
+ if (i < 0) {
+ this.hasMore = false;
+ return;
+ }
+
+ // find the least greater element to the right of the dip
+ int leastToRightIndex = i + 1;
+ for (int j = i + 2; j < this.n; j++){
+ if (this.index[j] < this.index[leastToRightIndex] && this.index[j] > this.index[i]) {
+ leastToRightIndex = j;
+ }
+ }
+
+ // switch dip element with least greater element to its right
+ int t = this.index[i];
+ this.index[i] = this.index[leastToRightIndex];
+ this.index[leastToRightIndex] = t;
+
+ if (this.m - 1 > i) {
+ // reverse the elements to the right of the dip
+ reverseAfter(this.index, i);
+ // reverse the elements to the right of m - 1
+ reverseAfter(this.index, this.m - 1);
+ }
+ }
+
+ private int rightmostDip() {
+ for (int i = this.n - 2; i >= 0; i--){
+ if (this.index[i] < this.index[i+1]){
+ return i;
+ }
+ }
+ return -1;
+ }
+ };
+ }
+
/**
* Determine if a "flip" would occur given the provided probability value.
* @param probability The provided probability value. This value must be in [0,1]
diff --git a/src/main/java/net/sourceforge/cilib/problem/dataset/MatrixDataSetBuilder.java b/src/main/java/net/sourceforge/cilib/problem/dataset/MatrixDataSetBuilder.java
index 43d97ea..2712b29 100644
--- a/src/main/java/net/sourceforge/cilib/problem/dataset/MatrixDataSetBuilder.java
+++ b/src/main/java/net/sourceforge/cilib/problem/dataset/MatrixDataSetBuilder.java
@@ -81,7 +81,7 @@ public class MatrixDataSetBuilder extends BinaryDataSetBuilder {
if (m <= 0)
throw new IllegalStateException("Need to have a positive number as the input dimensions");
- matrixBuilder = Matrix.builder().rows(numvectors).columns(m);
+ matrixBuilder = Matrix.builder().dimensions(numvectors, m);
if (tok.nextToken() != StreamTokenizer.TT_NUMBER)
throw new IllegalStateException("Expected an integer number as the third token in the dataset");
diff --git a/src/main/java/net/sourceforge/cilib/problem/mappingproblem/MappingProblem.java b/src/main/java/net/sourceforge/cilib/problem/mappingproblem/MappingProblem.java
index 71a6f73..57d529f 100644
--- a/src/main/java/net/sourceforge/cilib/problem/mappingproblem/MappingProblem.java
+++ b/src/main/java/net/sourceforge/cilib/problem/mappingproblem/MappingProblem.java
@@ -195,7 +195,7 @@ public abstract class MappingProblem extends OptimisationProblemAdapter {
MatrixDataSetBuilder matrixDataSetBuilder = (MatrixDataSetBuilder) dataSetBuilder;
inputs = matrixDataSetBuilder.getMatrix();
- Matrix.Builder matrixBuilder = Matrix.builder().rows(numvectors).columns(numvectors);
+ Matrix.Builder matrixBuilder = Matrix.builder().dimensions(numvectors, numvectors);
for(int i = 0; i < numvectors; i++) {
for(int j = 0; j < i; j++) {
diff --git a/src/main/java/net/sourceforge/cilib/type/types/container/Matrix.java b/src/main/java/net/sourceforge/cilib/type/types/container/Matrix.java
index 807f7ff..12ca4db 100644
--- a/src/main/java/net/sourceforge/cilib/type/types/container/Matrix.java
+++ b/src/main/java/net/sourceforge/cilib/type/types/container/Matrix.java
@@ -198,11 +198,22 @@ public final class Matrix implements Type {
}
/**
- * Obtain the determinant of the current matrix.
+ * Obtain the determinant of the current matrix. The determinant is obtained
+ * by creating a LU decomposition matrix. The determinant is then calculated by:
+ * <p>
+ * <pre>
+ * \begin{equation}
+ * \mbox{det($A$)} = |L| \dot |U|
+ * \end{equation}
+ * </pre>
+ * where |L| = 1 and |U| is the product of the diagonal elements.
* @return The determinant value.
*/
public double determinant() {
- return contents[0][0]*contents[1][1] - contents[1][0]*contents[0][1];
+ Preconditions.checkState(isSquare(), "Cannot obtain determinant of a non-square matrix");
+
+ LUDecomposition decomposition = new LUDecomposition(this);
+ return decomposition.determinant();
}
/**
@@ -260,22 +271,19 @@ public final class Matrix implements Type {
}
/**
- * Define the number of rows that the built up {@code Matrix} will contain.
- * @param rows The required number of rows.
- * @return The current {@code Builder}.
+ * Define the dimensions (rows and coloums) that the bult up {@code Matrix} will
+ * contain.
+ * @param rows The number of rows.
+ * @param columns The number of columns.
+ * @return The current {@code Builder)
+ * @throws IllegalArgumentException if {@code rows} or {@code columns} are less than 1.
*/
- public Builder rows(int rows) {
- this.rowNumber = rows;
- return this;
- }
+ public Builder dimensions(int rows, int columns) {
+ Preconditions.checkArgument(rows >= 1);
+ Preconditions.checkArgument(columns >= 1);
- /**
- * Define the number of coloumns that the built up {@code Matrix} will contain.
- * @param colummns The required number of columns.
- * @return The current {@code Builder}.
- */
- public Builder columns(int colummns) {
- this.colNumber = colummns;
+ this.rowNumber = rows;
+ this.colNumber = columns;
return this;
}
@@ -410,4 +418,98 @@ public final class Matrix implements Type {
return y;
}
}
+
+ /**
+ * This class is an implementation of the Doolittle / Crout algorithms.
+ * <p>
+ * Please refer to <a href=http://en.wikipedia.org/wiki/LU_decomposition#Doolittle_algorithm>
+ * the Wikipedia entry</a> for more information.
+ * <p>
+ * Additional help and code was obtained from the JAMA project (which seems to
+ * be stagnant) and from the commons-math.
+ */
+ private static final class LUDecomposition {
+ private double[][] lu;
+
+ private int m;
+ private int n;
+ private int pivsign;
+ private int[] piv;
+
+ private LUDecomposition(Matrix matrix) {
+ m = matrix.getRows();
+ n = matrix.getColumns();
+ pivsign = 1;
+
+ lu = new double[m][n];
+ for (int i = 0; i < m; i++) {
+ for (int j = 0; j < n; j++) {
+ lu[i][j] = matrix.contents[i][j];
+ }
+ }
+
+ piv = new int[m];
+ for (int i = 0; i < m; i++)
+ piv[i] = i;
+
+ double[] LUrowi;
+ double[] LUcolj = new double[m];
+
+ // Outer-loop
+ for (int j = 0; j < n; j++) {
+ // Make a copy of the j-th column.
+ for (int i = 0; i < m; i++) {
+ LUcolj[i] = lu[i][j];
+ }
+
+ // Apply previous transformations
+ for (int i = 0; i < m; i++) {
+ LUrowi = lu[i];
+
+ int kmax = Math.min(i, j);
+ double s = 0.0;
+ for (int k = 0; k < kmax; k++) {
+ s += LUrowi[k]*LUcolj[k];
+ }
+
+ LUcolj[i] -= s;
+ LUrowi[j] = LUcolj[i];
+ }
+
+ // Find the pivot and exchange if needed
+ int p = j;
+ for (int i = j+1; i < m; i++) {
+ if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
+ p = i;
+ }
+ }
+
+ if (p != j) {
+ for (int k = 0; k < n; k++) {
+ double t = lu[p][k]; lu[p][k] = lu[j][k]; lu[j][k] = t;
+ }
+ int k = piv[p]; piv[p] = piv[j]; piv[j] = k;
+ pivsign = -pivsign;
+ }
+
+ // Compute multipliers.
+ if (j < m & lu[j][j] != 0.0) {
+ for (int i = j+1; i < m; i++) {
+ lu[i][j] /= lu[j][j];
+ }
+ }
+ }
+ }
+
+ public double determinant() {
+ Preconditions.checkState(m == n, "Matrix must be square.");
+
+ double d = Integer.valueOf(pivsign).doubleValue();
+ for (int j = 0; j < n; j++) {
+ d *= lu[j][j];
+ }
+
+ return d;
+ }
+ }
}
diff --git a/src/test/java/net/sourceforge/cilib/math/MathsTest.java b/src/test/java/net/sourceforge/cilib/math/MathsTest.java
index 96f98e1..a774623 100644
--- a/src/test/java/net/sourceforge/cilib/math/MathsTest.java
+++ b/src/test/java/net/sourceforge/cilib/math/MathsTest.java
@@ -21,10 +21,15 @@
*/
package net.sourceforge.cilib.math;
-import static org.junit.Assert.assertEquals;
-import static org.junit.Assert.fail;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Iterator;
+import java.util.List;
+import org.junit.Assert;
import org.junit.Test;
+import static org.hamcrest.CoreMatchers.is;
+
/**
*
@@ -39,38 +44,50 @@ public class MathsTest {
@Test
public void testFactorial() {
- assertEquals(1.0, Maths.factorial(0.0), Double.MIN_NORMAL);
- assertEquals(1.0, Maths.factorial(1.0), Double.MIN_NORMAL);
- assertEquals(6.0, Maths.factorial(3), Double.MIN_NORMAL);
- assertEquals(720.0, Maths.factorial(6), Double.MIN_NORMAL);
- assertEquals(9.33262154439441E157, Maths.factorial(100), Double.MIN_NORMAL);
+ Assert.assertEquals(1.0, Maths.factorial(0.0), Double.MIN_NORMAL);
+ Assert.assertEquals(1.0, Maths.factorial(1.0), Double.MIN_NORMAL);
+ Assert.assertEquals(6.0, Maths.factorial(3), Double.MIN_NORMAL);
+ Assert.assertEquals(720.0, Maths.factorial(6), Double.MIN_NORMAL);
+ Assert.assertEquals(9.33262154439441E157, Maths.factorial(100), Double.MIN_NORMAL);
}
@Test
public void testCombination() {
- assertEquals(792.0, Maths.combination(12, 5), Double.MIN_NORMAL);
+ Assert.assertEquals(792.0, Maths.combination(12, 5), Double.MIN_NORMAL);
+ }
- try {
- Maths.combination(-1, -5);
- fail("Invalid input!");
- }
- catch (Exception e) {}
+ @Test(expected=IllegalArgumentException.class)
+ public void combinationInvalidN() {
+ Maths.combination(-1, 5);
+ }
- try {
- Maths.combination(-1, 5);
- fail("Invalid input!");
- }
- catch (Exception e) {}
+ @Test(expected=IllegalArgumentException.class)
+ public void combinationInvalidR() {
+ Maths.combination(1, -5);
+ }
- try {
- Maths.combination(1, -5);
- fail("Invalid input!");
+ @Test
+ public void combinationSpecialCase() {
+ Assert.assertEquals(1.0, Maths.combination(0, 0), Double.MIN_NORMAL);
+ Assert.assertEquals(1.0, Maths.combination(1, 0), Double.MIN_NORMAL);
+ Assert.assertEquals(1.0, Maths.combination(1, 1), Double.MIN_NORMAL);
+ }
+
+ @Test
+ public void listPermutation() {
+ List<Integer> numbers = Arrays.asList(1, 2);
+ List<List<Integer>> permutationList = new ArrayList<List<Integer>>();
+
+ for (Iterator<List<Integer>> permutations = Maths.permutation(numbers, 2); permutations.hasNext(); ) {
+ permutationList.add(permutations.next());
}
- catch (Exception e) {}
- assertEquals(1.0, Maths.combination(0, 0), Double.MIN_NORMAL);
- assertEquals(1.0, Maths.combination(1, 0), Double.MIN_NORMAL);
- assertEquals(1.0, Maths.combination(1, 1), Double.MIN_NORMAL);
+ Assert.assertThat(permutationList.size(), is(2));
+
+ List<Integer> expected1 = Arrays.asList(1, 2);
+ List<Integer> expected2 = Arrays.asList(2, 1);
+ Assert.assertTrue(permutationList.contains(expected1));
+ Assert.assertTrue(permutationList.contains(expected2));
}
}
diff --git a/src/test/java/net/sourceforge/cilib/type/types/container/MatrixTest.java b/src/test/java/net/sourceforge/cilib/type/types/container/MatrixTest.java
index 0c2c98e..498951d 100644
--- a/src/test/java/net/sourceforge/cilib/type/types/container/MatrixTest.java
+++ b/src/test/java/net/sourceforge/cilib/type/types/container/MatrixTest.java
@@ -35,29 +35,29 @@ public class MatrixTest {
@Test(expected=IllegalArgumentException.class)
public void constructionZeroRow() {
- Matrix.builder().rows(0).build();
+ Matrix.builder().dimensions(0, 1).build();
}
@Test(expected=IllegalArgumentException.class)
public void constructionZeroColumn() {
- Matrix.builder().columns(0).build();
+ Matrix.builder().dimensions(1, 0).build();
}
@Test
public void square() {
- Matrix a = Matrix.builder().rows(2).columns(2).build();
+ Matrix a = Matrix.builder().dimensions(2, 2).build();
Assert.assertTrue(a.isSquare());
}
@Test
public void notSquare() {
- Matrix a = Matrix.builder().rows(3).columns(4).build();
+ Matrix a = Matrix.builder().dimensions(3, 4).build();
Assert.assertFalse(a.isSquare());
}
@Test
public void valueAt() {
- Matrix a = Matrix.builder().rows(1).columns(2)
+ Matrix a = Matrix.builder().dimensions(1, 2)
.addRow(1.0, 2.0)
.build();
@@ -67,13 +67,13 @@ public class MatrixTest {
@Test(expected=IndexOutOfBoundsException.class)
public void invalidValueOf() {
- Matrix a = Matrix.builder().rows(1).columns(1).build();
+ Matrix a = Matrix.builder().dimensions(1, 1).build();
a.valueAt(1, 2);
}
@Test
public void getRow() {
- Matrix a = Matrix.builder().rows(2).columns(2)
+ Matrix a = Matrix.builder().dimensions(2, 2)
.addRow(1.0, 1.0)
.addRow(2.0, 2.0)
.build();
@@ -86,23 +86,23 @@ public class MatrixTest {
@Test
public void rowNumber() {
- Matrix a = Matrix.builder().rows(4).columns(5).build();
+ Matrix a = Matrix.builder().dimensions(4, 5).build();
Assert.assertThat(a.getRows(), is(4));
}
@Test
public void columnNumber() {
- Matrix a = Matrix.builder().rows(5).columns(8).build();
+ Matrix a = Matrix.builder().dimensions(5, 8).build();
Assert.assertThat(a.getColumns(), is(8));
}
@Test
public void addition() {
- Matrix a = Matrix.builder().rows(2).columns(2)
+ Matrix a = Matrix.builder().dimensions(2, 2)
.addRow(1.0, 2.0)
.addRow(3.0, 4.0)
.build();
- Matrix b = Matrix.builder().rows(2).columns(2)
+ Matrix b = Matrix.builder().dimensions(2, 2)
.addRow(1.0, 2.0)
.addRow(3.0, 4.0)
.build();
@@ -115,18 +115,18 @@ public class MatrixTest {
@Test(expected=IllegalArgumentException.class)
public void invalidAddition() {
- Matrix a = Matrix.builder().rows(3).columns(2).build();
- Matrix b = Matrix.builder().rows(1).columns(2).build();
+ Matrix a = Matrix.builder().dimensions(3, 2).build();
+ Matrix b = Matrix.builder().dimensions(1, 2).build();
a.plus(b);
}
@Test
public void subtraction() {
- Matrix a = Matrix.builder().rows(2).columns(2)
+ Matrix a = Matrix.builder().dimensions(2, 2)
.addRow(2.0, 4.0)
.addRow(6.0, 8.0)
.build();
- Matrix b = Matrix.builder().rows(2).columns(2)
+ Matrix b = Matrix.builder().dimensions(2, 2)
.addRow(1.0, 2.0)
.addRow(3.0, 4.0)
.build();
@@ -139,18 +139,18 @@ public class MatrixTest {
@Test(expected=IllegalArgumentException.class)
public void invalidSubtraction() {
- Matrix a = Matrix.builder().rows(2).columns(2).build();
- Matrix b = Matrix.builder().rows(2).columns(3).build();
+ Matrix a = Matrix.builder().dimensions(2, 2).build();
+ Matrix b = Matrix.builder().dimensions(2, 3).build();
a.minus(b);
}
@Test
public void multiplication() {
- Matrix a = Matrix.builder().rows(2).columns(2)
+ Matrix a = Matrix.builder().dimensions(2, 2)
.addRow(2.0, 4.0)
.addRow(6.0, 8.0)
.build();
- Matrix b = Matrix.builder().rows(2).columns(2)
+ Matrix b = Matrix.builder().dimensions(2, 2)
.addRow(1.0, 2.0)
.addRow(3.0, 4.0)
.build();
@@ -163,7 +163,7 @@ public class MatrixTest {
@Test
public void squareTranspose() {
- Matrix a = Matrix.builder().rows(2).columns(2)
+ Matrix a = Matrix.builder().dimensions(2, 2)
.addRow(2.0, 4.0)
.addRow(6.0, 8.0)
.build();
@@ -176,7 +176,7 @@ public class MatrixTest {
@Test
public void transposeRowVector() {
- Matrix a = Matrix.builder().rows(1).columns(2)
+ Matrix a = Matrix.builder().dimensions(1, 2)
.addRow(2.0, 4.0)
.build();
@@ -188,7 +188,7 @@ public class MatrixTest {
@Test
public void transposeColumnVector() {
- Matrix a = Matrix.builder().rows(2).columns(1)
+ Matrix a = Matrix.builder().dimensions(2, 1)
.addRow(2.0)
.addRow(4.0)
.build();
@@ -201,7 +201,7 @@ public class MatrixTest {
@Test
public void identity() {
- Matrix identity = Matrix.builder().rows(4).columns(4).identity().build();
+ Matrix identity = Matrix.builder().dimensions(4, 4).identity().build();
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
@@ -215,7 +215,7 @@ public class MatrixTest {
public void rotation() {
double angle = Math.PI / 4.0;
- Matrix matrix = Matrix.builder().rows(2).columns(2).addRow(1.0, 1.0).addRow(1.0, 1.0).build();
+ Matrix matrix = Matrix.builder().dimensions(2, 2).addRow(1.0, 1.0).addRow(1.0, 1.0).build();
Matrix result = matrix.rotate(angle);
Assert.assertThat(result.valueAt(0, 0), is(1.414213562373095));
@@ -226,18 +226,19 @@ public class MatrixTest {
@Test
public void determinant() {
- Matrix matrix = Matrix.builder().rows(2).columns(2).addRow(1.0, 1.0).addRow(1.0, 1.0).build();
- Assert.assertThat(matrix.determinant(), is(0.0));
+ Matrix matrix = Matrix.builder().dimensions(3, 3).addRow(2.0, 1.0, 0.0).addRow(1.0, 2.0, -1.0).addRow(3.0, 2.0, 1.0).build();
+ double value = matrix.determinant();
+ Assert.assertEquals(4.0, value, 0.00001);
}
@Test(expected=IllegalStateException.class)
public void invalidIdentity() {
- Matrix.builder().rows(2).columns(5).identity().build();
+ Matrix.builder().dimensions(2, 5).identity().build();
}
@Test
public void uniquePositionSetting() {
- Matrix a = Matrix.builder().rows(2).columns(2)
+ Matrix a = Matrix.builder().dimensions(2, 2)
.valueAt(0, 0, 3.0)
.build();
@@ -250,8 +251,8 @@ public class MatrixTest {
@Test
public void equal() {
Matrix.Builder builder = Matrix.builder();
- Matrix a = builder.rows(1).columns(1).valueAt(0, 0, 2.0).build();
- Matrix b = builder.rows(1).columns(1).valueAt(0, 0, 2.0).build();
+ Matrix a = builder.dimensions(1, 1).valueAt(0, 0, 2.0).build();
+ Matrix b = builder.dimensions(1, 1).valueAt(0, 0, 2.0).build();
Assert.assertTrue(a.equals(b));
}
@@ -259,8 +260,8 @@ public class MatrixTest {
@Test
public void hash() {
Matrix.Builder builder = Matrix.builder();
- Matrix a = builder.rows(1).columns(1).valueAt(0, 0, 2.0).build();
- Matrix b = builder.rows(1).columns(1).valueAt(0, 0, 2.0).build();
+ Matrix a = builder.dimensions(1, 1).valueAt(0, 0, 2.0).build();
+ Matrix b = builder.dimensions(1, 1).valueAt(0, 0, 2.0).build();
Assert.assertTrue(a.hashCode() == b.hashCode());
}
--
1.6.4
|