halfdan/eva2
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

removed classpath

Browse Source
This commit is contained in:
Andreas Dräger 2010-08-09 11:03:36 +00:00
parent 0ec063bcff
commit e12ab22e31
4 changed files with 660 additions and 548 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.classpath
View File

@ -1,6 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<classpath>
<classpathentry kind="src" path="src"/>
<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER"/>
<classpathentry kind="output" path="build"/>
</classpath>

374
src/eva2/tools/math/Mathematics.java
View File

@ -23,8 +23,10 @@ public class Mathematics {
* @return * @return
*/ */
public static double[][] adjoint(double[][] a) { public static double[][] adjoint(double[][] a) {
if (a == null) return null; if (a == null)
if (a.length != a[0].length) return null; return null;
if (a.length != a[0].length)
return null;
double[][] b = new double[a.length][a.length]; double[][] b = new double[a.length][a.length];
for (int i = 0; i < a.length; i++) for (int i = 0; i < a.length; i++)
for (int j = 0; j < a.length; j++) for (int j = 0; j < a.length; j++)
@ -52,9 +54,12 @@ public class Mathematics {
* is not square). * is not square).
*/ */
public static double determinant(double[][] matrix) { public static double determinant(double[][] matrix) {
if (matrix == null) return 0; if (matrix == null)
if (matrix.length != matrix[0].length) return 0; return 0;
if (matrix.length == 1) return matrix[0][0]; if (matrix.length != matrix[0].length)
return 0;
if (matrix.length == 1)
return matrix[0][0];
if (matrix.length == 2) if (matrix.length == 2)
return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]; return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
if (matrix.length == 3) if (matrix.length == 3)
@ -66,7 +71,8 @@ public class Mathematics {
double det = 0; double det = 0;
for (int k = 0; k < matrix.length; k++) { for (int k = 0; k < matrix.length; k++) {
if (matrix[0][k] != 0) det += matrix[0][k] * adjoint(matrix, 0, k); if (matrix[0][k] != 0)
det += matrix[0][k] * adjoint(matrix, 0, k);
} }
return det; return det;
} }
@ -83,13 +89,14 @@ public class Mathematics {
* what kind of distance funktion * what kind of distance funktion
* @return the distance of x and y * @return the distance of x and y
* @throws Exception * @throws Exception
* if x and y have different dimensions an exception is * if x and y have different dimensions an exception is thrown.
* thrown.
*/ */
public static double dist(double[] x, double[] y, int root) { public static double dist(double[] x, double[] y, int root) {
if (x.length != y.length) if (x.length != y.length)
throw new RuntimeException("The vectors x and y must have the same dimension"); throw new RuntimeException(
if (root == 0) throw new RuntimeException("There is no 0-root!"); "The vectors x and y must have the same dimension");
if (root == 0)
throw new RuntimeException("There is no 0-root!");
double d = 0; double d = 0;
for (int i = 0; i < x.length; i++) for (int i = 0; i < x.length; i++)
d += Math.pow(Math.abs(x[i] - y[i]), root); d += Math.pow(Math.abs(x[i] - y[i]), root);
@ -107,12 +114,12 @@ public class Mathematics {
* what kind of distance funktion * what kind of distance funktion
* @return the distance of x and y * @return the distance of x and y
* @throws Exception * @throws Exception
* if x and y have different dimensions an exception is * if x and y have different dimensions an exception is thrown.
* thrown.
*/ */
public static double euclidianDist(double[] x, double[] y) { public static double euclidianDist(double[] x, double[] y) {
if (x.length != y.length) if (x.length != y.length)
throw new RuntimeException("The vectors x and y must have the same dimension"); throw new RuntimeException(
"The vectors x and y must have the same dimension");
double d = 0; double d = 0;
for (int i = 0; i < x.length; i++) for (int i = 0; i < x.length; i++)
d += Math.pow(Math.abs(x[i] - y[i]), 2); d += Math.pow(Math.abs(x[i] - y[i]), 2);
@ -120,8 +127,8 @@ public class Mathematics {
} }
/** /**
* Expand a vector to a higher dimension (len) by filling it up * Expand a vector to a higher dimension (len) by filling it up with a
* with a constant value. * constant value.
* *
* @param x * @param x
* @param len * @param len
@ -130,12 +137,14 @@ public class Mathematics {
*/ */
public static double[] expandVector(double[] x, int len, double v) { public static double[] expandVector(double[] x, int len, double v) {
if (len <= x.length) {// not really an error, just perform identity if (len <= x.length) {// not really an error, just perform identity
// System.err.println("Error, invalid length in expandVector, expecting l>" + x.length); // System.err.println("Error, invalid length in expandVector, expecting l>"
// + x.length);
return x; return x;
} else { } else {
double[] expanded = new double[len]; double[] expanded = new double[len];
System.arraycopy(x, 0, expanded, 0, x.length); System.arraycopy(x, 0, expanded, 0, x.length);
for (int i=x.length; i<expanded.length; i++) expanded[i] = v; for (int i = x.length; i < expanded.length; i++)
expanded[i] = v;
return expanded; return expanded;
} }
} }
@ -190,15 +199,18 @@ public class Mathematics {
*/ */
public static double getAvgRange(double[][] range) { public static double getAvgRange(double[][] range) {
double sum = 0.; double sum = 0.;
for (int i=0; i<range.length; i++) sum+=(range[i][1]-range[i][0]); for (int i = 0; i < range.length; i++)
sum += (range[i][1] - range[i][0]);
return sum / range.length; return sum / range.length;
} }
/** /**
* Calculates the norm of the given vector relative to the problem range. * Calculates the norm of the given vector relative to the problem range.
* *
* @param vector a double vector within the range * @param vector
* @param range the range in each dimension * a double vector within the range
* @param range
* the range in each dimension
* @return measure of the length relative to the problem range * @return measure of the length relative to the problem range
*/ */
public static double getRelativeLength(double[] vector, double[][] range) { public static double getRelativeLength(double[] vector, double[][] range) {
@ -221,7 +233,8 @@ public class Mathematics {
* @param j * @param j
* @param w * @param w
*/ */
public static void getRotationEntriesSingleAxis(Matrix tmp, int i, int j, double w) { public static void getRotationEntriesSingleAxis(Matrix tmp, int i, int j,
double w) {
double cosw = Math.cos(w); double cosw = Math.cos(w);
double sinw = Math.sin(w); double sinw = Math.sin(w);
tmp.set(i, i, cosw); tmp.set(i, i, cosw);
@ -250,21 +263,28 @@ public class Mathematics {
} }
/** /**
* Return a matrix A which performs the rotation of vec to (1,0,0,...0) if forward is true, else * Return a matrix A which performs the rotation of vec to (1,0,0,...0) if
* return a matrix B which performs the reverted rotation, where B=A' (transposition). * forward is true, else return a matrix B which performs the reverted
* rotation, where B=A' (transposition).
* *
* @param vec * @param vec
* @return * @return
*/ */
public static Matrix getRotationMatrix(Matrix vec) { public static Matrix getRotationMatrix(Matrix vec) {
Matrix A = Matrix.identity(vec.getRowDimension(), vec.getRowDimension()); Matrix A = Matrix
Matrix tmp = Matrix.identity(vec.getRowDimension(), vec.getRowDimension()); .identity(vec.getRowDimension(), vec.getRowDimension());
Matrix tmp = Matrix.identity(vec.getRowDimension(), vec
.getRowDimension());
Matrix z = (Matrix) vec.clone(); Matrix z = (Matrix) vec.clone();
z.multi(1. / z.norm2()); // normalize z.multi(1. / z.norm2()); // normalize
for (int i = 1; i < vec.getRowDimension(); i++) { for (int i = 1; i < vec.getRowDimension(); i++) {
double w = Math.atan2(z.get(i,0), z.get(0,0));// calc angle between the projection of x and x0 in x0-xi-plane double w = Math.atan2(z.get(i, 0), z.get(0, 0));// calc angle
// between the
// projection of x
// and x0 in
// x0-xi-plane
// System.out.println("deg: "+(w/Math.PI)*180); // System.out.println("deg: "+(w/Math.PI)*180);
// make partial rotation matrix // make partial rotation matrix
@ -281,8 +301,11 @@ public class Mathematics {
/** /**
* This method return a vector from a to b * This method return a vector from a to b
* @param a first vector *
* @param b second vectors * @param a
* first vector
* @param b
* second vectors
* @return the vector from a to b * @return the vector from a to b
*/ */
public static double[] getVectorFromTo(double[] a, double[] b) { public static double[] getVectorFromTo(double[] a, double[] b) {
@ -299,11 +322,13 @@ public class Mathematics {
* @param f1 * @param f1
* @return * @return
*/ */
public static double hyperbolicInterpolation(double x, double x0, double x1, public static double hyperbolicInterpolation(double x, double x0,
double f0, double f1) { double x1, double f0, double f1) {
if (x1 == 0) return lerp(f0, f1, (x - x0) / (-x0)); if (x1 == 0)
return lerp(f0, f1, (x - x0) / (-x0));
double l = lerp(x0 / x1, 1, x); double l = lerp(x0 / x1, 1, x);
if (l == 0) return linearInterpolation(x, x0, x1, f0, f1); if (l == 0)
return linearInterpolation(x, x0, x1, f0, f1);
return lerp(f0, f1, x / l); return lerp(f0, f1, x / l);
} }
@ -314,7 +339,8 @@ public class Mathematics {
* @param makeRange * @param makeRange
* @param destRange * @param destRange
*/ */
public static void intersectRange(double[][] r1, double[][] r2, double[][] destRange) { public static void intersectRange(double[][] r1, double[][] r2,
double[][] destRange) {
for (int i = 0; i < r1.length && i < r2.length; i++) { for (int i = 0; i < r1.length && i < r2.length; i++) {
destRange[i][0] = Math.max(r1[i][0], r2[i][0]); destRange[i][0] = Math.max(r1[i][0], r2[i][0]);
destRange[i][1] = Math.min(r1[i][1], r2[i][1]); destRange[i][1] = Math.min(r1[i][1], r2[i][1]);
@ -329,11 +355,14 @@ public class Mathematics {
* @return * @return
*/ */
public static double[][] inverse(double[][] a) { public static double[][] inverse(double[][] a) {
if (a == null) return null; if (a == null)
if (a.length != a[0].length) return null; return null;
if (a.length != a[0].length)
return null;
double det = determinant(a); double det = determinant(a);
if (det == 0) return null; if (det == 0)
return null;
double[][] b = adjoint(a); double[][] b = adjoint(a);
for (int i = 0; i < a.length; i++) for (int i = 0; i < a.length; i++)
for (int j = 0; j < a.length; j++) for (int j = 0; j < a.length; j++)
@ -342,14 +371,16 @@ public class Mathematics {
} }
/** /**
* Check whether the given value lies within the interval in every dimension. * Check whether the given value lies within the interval in every
* dimension.
* *
* @param x * @param x
* @param range * @param range
* @return true if the vector lies within the range, else false * @return true if the vector lies within the range, else false
*/ */
public static boolean isInRange(double v, double lower, double upper) { public static boolean isInRange(double v, double lower, double upper) {
if (v<lower || (v>upper)) return false; if (v < lower || (v > upper))
return false;
return true; return true;
} }
@ -362,25 +393,30 @@ public class Mathematics {
*/ */
public static boolean isInRange(double[] x, double[][] range) { public static boolean isInRange(double[] x, double[][] range) {
for (int i = 0; i < x.length; i++) { for (int i = 0; i < x.length; i++) {
if (x[i]<range[i][0] || (x[i]>range[i][1])) return false; if (x[i] < range[i][0] || (x[i] > range[i][1]))
return false;
} }
return true; return true;
} }
/** /**
* Returns false if a vector contains NaN, its squared sum is NaN * Returns false if a vector contains NaN, its squared sum is NaN or the
* or the absolute sum is smaller than 10^-18. * absolute sum is smaller than 10^-18.
*
* @param d * @param d
* @return * @return
*/ */
public static boolean isValidVec(double[] d) { public static boolean isValidVec(double[] d) {
double sum = 0; double sum = 0;
for (int i = 0; i < d.length; i++) { for (int i = 0; i < d.length; i++) {
if (Double.isNaN(d[i])) return false; if (Double.isNaN(d[i]))
return false;
sum += Math.pow(d[i], 2); sum += Math.pow(d[i], 2);
} }
if (Double.isNaN(sum)) return false; if (Double.isNaN(sum))
if (Math.abs(sum) < 0.000000000000000001) return false; return false;
if (Math.abs(sum) < 0.000000000000000001)
return false;
return true; return true;
} }
@ -412,20 +448,23 @@ public class Mathematics {
*/ */
public static double linearInterpolation(double x, double x0, double x1, public static double linearInterpolation(double x, double x0, double x1,
double f0, double f1) { double f0, double f1) {
if (x1 == x0) return f0; if (x1 == x0)
return f0;
return lerp(f0, f1, (x - x0) / (x1 - x0)); return lerp(f0, f1, (x - x0) / (x1 - x0));
} }
public static double max(double[] vals) { public static double max(double[] vals) {
double maxVal = vals[0]; double maxVal = vals[0];
for (int i=1; i<vals.length; i++) maxVal = Math.max(maxVal, vals[i]); for (int i = 1; i < vals.length; i++)
maxVal = Math.max(maxVal, vals[i]);
return maxVal; return maxVal;
} }
/** /**
* Computes the mean for an array of doubles. * Computes the mean for an array of doubles.
* *
* @param vector the array * @param vector
* the array
* @return the mean * @return the mean
*/ */
public static double mean(double[] vector) { public static double mean(double[] vector) {
@ -435,8 +474,11 @@ public class Mathematics {
return sum(vector) / (double) vector.length; return sum(vector) / (double) vector.length;
} }
/** This method returns a mean vector from a whole array of vectors. /**
* @param d d[i] the vectors, d[i][j] the jth coordinate of the ith vector * This method returns a mean vector from a whole array of vectors.
*
* @param d
* d[i] the vectors, d[i][j] the jth coordinate of the ith vector
* @return The mean vector. * @return The mean vector.
*/ */
public static double[] meanVect(double[][] d) { public static double[] meanVect(double[][] d) {
@ -455,40 +497,57 @@ public class Mathematics {
/** /**
* Computes the median of a given double vector by sorting x. * Computes the median of a given double vector by sorting x.
* *
* @param x a vector of doubles * @param x
* @param cloneX flag whether x should be cloned before sorting. * a vector of doubles
* @param cloneX
* flag whether x should be cloned before sorting.
* @return the median * @return the median
*/ */
public static double median(double[] x, boolean cloneX) { public static double median(double[] x, boolean cloneX) {
double[] in; double[] in;
if (cloneX) in = (double[]) x.clone(); if (cloneX)
else in = x; in = (double[]) x.clone();
else
in = x;
if (in.length==1) return in[0]; if (in.length == 1)
else if (in.length==2) return (in[0]+in[1])/2.; return in[0];
else if (in.length == 2)
return (in[0] + in[1]) / 2.;
else { else {
Arrays.sort(in); Arrays.sort(in);
if (in.length % 2 != 0) return in[(in.length-1) / 2]; if (in.length % 2 != 0)
else return (in[in.length/2] + in[(in.length/2)+1]) / 2.; return in[(in.length - 1) / 2];
else
return (in[in.length / 2] + in[(in.length / 2) + 1]) / 2.;
} }
} }
/** /**
* Computes the median of a given list of double vectors by sorting it. * Computes the median of a given list of double vectors by sorting it. If
* If the size is even, no direct median is defined - in that case it may be * the size is even, no direct median is defined - in that case it may be
* interpolated by the two closest elements or one of them may be selected (always * interpolated by the two closest elements or one of them may be selected
* the smaller one depending on the comparator. * (always the smaller one depending on the comparator.
* *
* @see #DoubleArrayComparator * @see #DoubleArrayComparator
* @param dblArrList a list of double vectors * @param dblArrList
* @param interpolate flag whether, for even size, the median is interpolated * a list of double vectors
* @param interpolate
* flag whether, for even size, the median is interpolated
* @return the median * @return the median
*/ */
public static double[] median(List<double[]> dblArrList, boolean interpolate) { public static double[] median(List<double[]> dblArrList, boolean interpolate) {
java.util.Collections.sort(dblArrList, new DoubleArrayComparator()); // by default, the comparator uses pareto dominance java.util.Collections.sort(dblArrList, new DoubleArrayComparator()); // by
// default,
// the
// comparator
// uses
// pareto
// dominance
int len = dblArrList.size(); int len = dblArrList.size();
if (len % 2 != 0) return dblArrList.get((len-1) / 2); if (len % 2 != 0)
return dblArrList.get((len - 1) / 2);
else { else {
double[] med = dblArrList.get(len / 2).clone(); double[] med = dblArrList.get(len / 2).clone();
if (interpolate) { if (interpolate) {
@ -501,14 +560,16 @@ public class Mathematics {
public static double min(double[] vals) { public static double min(double[] vals) {
double minVal = vals[0]; double minVal = vals[0];
for (int i=1; i<vals.length; i++) minVal = Math.min(minVal, vals[i]); for (int i = 1; i < vals.length; i++)
minVal = Math.min(minVal, vals[i]);
return minVal; return minVal;
} }
/** /**
* Computes the 2-norm of an array of doubles. * Computes the 2-norm of an array of doubles.
* *
* @param doubles the array of double * @param doubles
* the array of double
* @return the 2-norm of the elements * @return the 2-norm of the elements
*/ */
public static double norm(double[] d) { public static double norm(double[] d) {
@ -520,10 +581,13 @@ public class Mathematics {
} }
/** /**
* Normalizes the doubles in the array by their sum, * Normalizes the doubles in the array by their sum, so that they add up to
* so that they add up to one. * one.
* @param doubles the array of double *
* @exception IllegalArgumentException if sum is Zero or NaN * @param doubles
* the array of double
* @exception IllegalArgumentException
* if sum is Zero or NaN
*/ */
public static double[] normalizeSum(double[] v) { public static double[] normalizeSum(double[] v) {
double[] res = new double[v.length]; double[] res = new double[v.length];
@ -532,11 +596,13 @@ public class Mathematics {
} }
/** /**
* Normalizes the doubles in the array by their sum, * Normalizes the doubles in the array by their sum, so that they add up to
* so that they add up to one. * one.
* *
* @param doubles the array of double * @param doubles
* @exception IllegalArgumentException if sum is Zero or NaN * the array of double
* @exception IllegalArgumentException
* if sum is Zero or NaN
*/ */
public static void normalizeSum(double[] v, double[] res) { public static void normalizeSum(double[] v, double[] res) {
svMult(1. / sum(v), v, res); svMult(1. / sum(v), v, res);
@ -564,6 +630,7 @@ public class Mathematics {
/** /**
* Return the product over a double vector. * Return the product over a double vector.
*
* @param vals * @param vals
* @return * @return
*/ */
@ -576,7 +643,8 @@ public class Mathematics {
} }
// /** // /**
// * Normalizes the doubles in the array using the given value so that they sum up to 1. // * Normalizes the doubles in the array using the given value so that they
// sum up to 1.
// * // *
// * @param doubles the array of double // * @param doubles the array of double
// * @param sum the value by which the doubles are to be normalized // * @param sum the value by which the doubles are to be normalized
@ -588,17 +656,18 @@ public class Mathematics {
// } // }
// if (sum == 0) { // if (sum == 0) {
// // Maybe this should just be a return. // // Maybe this should just be a return.
// throw new IllegalArgumentException("Can't normalize array. Sum is zero."); // throw new
// IllegalArgumentException("Can't normalize array. Sum is zero.");
// } // }
// svMult(1/sum, v, res); // svMult(1/sum, v, res);
// } // }
/** /**
* Project the values in x to the range given. The range must be an vector of 2d-arrays * Project the values in x to the range given. The range must be an vector
* each of which containing lower and upper bound in the i-th dimension. * of 2d-arrays each of which containing lower and upper bound in the i-th
* x must not be longer than the available ranges. * dimension. x must not be longer than the available ranges. Values
* Values exceeding the bounds are set on the bound. * exceeding the bounds are set on the bound. The number of bound violations
* The number of bound violations is returned. * is returned.
* *
* @param x * @param x
* @param range * @param range
@ -606,7 +675,9 @@ public class Mathematics {
*/ */
public static int projectToRange(double[] x, double[][] range) { public static int projectToRange(double[] x, double[][] range) {
int viols = 0; int viols = 0;
if (x.length>range.length) System.err.println("Invalid vector length, x is longer than range! (Mathematics.projectToRange)"); if (x.length > range.length)
System.err
.println("Invalid vector length, x is longer than range! (Mathematics.projectToRange)");
for (int i = 0; i < x.length; i++) { for (int i = 0; i < x.length; i++) {
if (x[i] < range[i][0]) { if (x[i] < range[i][0]) {
viols++; viols++;
@ -632,15 +703,18 @@ public class Mathematics {
return min; return min;
} else if (v > max) { } else if (v > max) {
return max; return max;
} else return v; } else
return v;
} }
/** /**
* Create a random vector, the components will be set to gaussian distributed * Create a random vector, the components will be set to gaussian
* values with mean zero and the given standard deviation. * distributed values with mean zero and the given standard deviation.
* *
* @param dim the desired dimension * @param dim
* @param stdDev the gaussian standard deviation * the desired dimension
* @param stdDev
* the gaussian standard deviation
* @return random vector * @return random vector
*/ */
public static double[] randomVector(int dim, double stdDev) { public static double[] randomVector(int dim, double stdDev) {
@ -665,17 +739,22 @@ public class Mathematics {
for (int i = 0; i < x.length; i++) { for (int i = 0; i < x.length; i++) {
double dimLen = range[i][1] - range[i][0]; double dimLen = range[i][1] - range[i][0];
if (dimLen <= 0.) { if (dimLen <= 0.) {
EVAERROR.errorMsgOnce("Error in reflectBounds: empty range! (possibly multiple errors)"); EVAERROR
.errorMsgOnce("Error in reflectBounds: empty range! (possibly multiple errors)");
} else { } else {
if (x[i] < range[i][0]) { if (x[i] < range[i][0]) {
viols++; viols++;
d = range[i][0] - x[i]; d = range[i][0] - x[i];
while (d > dimLen) d -= dimLen; // avoid violating the other bound immediately while (d > dimLen)
d -= dimLen; // avoid violating the other bound
// immediately
x[i] = range[i][0] + d; x[i] = range[i][0] + d;
} else if (x[i] > range[i][1]) { } else if (x[i] > range[i][1]) {
viols++; viols++;
d = x[i] - range[i][1]; d = x[i] - range[i][1];
while (d>dimLen) d -= dimLen; // avoid violating the other bound immediately while (d > dimLen)
d -= dimLen; // avoid violating the other bound
// immediately
x[i] = range[i][1] - d; x[i] = range[i][1] - d;
} }
} }
@ -684,8 +763,8 @@ public class Mathematics {
} }
/** /**
* Simple version of reflection of a value moving by a step and bouncing * Simple version of reflection of a value moving by a step and bouncing of
* of min and max values like a pool ball. Precondition is min <= val <= max, * min and max values like a pool ball. Precondition is min <= val <= max,
* post condition is min <= retVal <= max. * post condition is min <= retVal <= max.
* *
* @param val * @param val
@ -694,8 +773,10 @@ public class Mathematics {
* @param max * @param max
* @return * @return
*/ */
public static double reflectValue(double val, double step, double min, double max) { public static double reflectValue(double val, double step, double min,
while (step > (max-min)) step -= (max-min); double max) {
while (step > (max - min))
step -= (max - min);
if ((val + step) > max) if ((val + step) > max)
return (2 * max - val - step); return (2 * max - val - step);
if ((val + step) < min) if ((val + step) < min)
@ -721,12 +802,14 @@ public class Mathematics {
public static double relDist(double[] x, double[] y, double def) public static double relDist(double[] x, double[] y, double def)
throws Exception { throws Exception {
if (x.length != y.length) if (x.length != y.length)
throw new Exception("The vectors x and y must have the same dimension"); throw new Exception(
"The vectors x and y must have the same dimension");
double d = 0; double d = 0;
for (int i = 0; i < x.length; i++) for (int i = 0; i < x.length; i++)
if (y[i] != 0) if (y[i] != 0)
d += Math.pow(((x[i] - y[i]) / y[i]), 2); d += Math.pow(((x[i] - y[i]) / y[i]), 2);
else d += def; else
d += def;
return d; return d;
} }
@ -745,7 +828,8 @@ public class Mathematics {
for (int i = 0; i < src.length; i++) { for (int i = 0; i < src.length; i++) {
dst[src.length - i - 1] = src[i]; dst[src.length - i - 1] = src[i];
} }
} else System.err.println("Mismatching array lengths!"); } else
System.err.println("Mismatching array lengths!");
} }
/** /**
@ -764,8 +848,8 @@ public class Mathematics {
} }
/** /**
* Rotate a given double vector using a rotation matrix. If the matrix * Rotate a given double vector using a rotation matrix. If the matrix is
* is null, x will be returned unchanged. Matrix dimensions must fit. * null, x will be returned unchanged. Matrix dimensions must fit.
* *
* @param x * @param x
* @param rotMatrix * @param rotMatrix
@ -776,7 +860,8 @@ public class Mathematics {
Matrix resVec = rotMatrix.times(new Matrix(x, x.length)); Matrix resVec = rotMatrix.times(new Matrix(x, x.length));
x = resVec.getColumnPackedCopy(); x = resVec.getColumnPackedCopy();
return x; return x;
} else return x; } else
return x;
} }
/** /**
@ -787,11 +872,14 @@ public class Mathematics {
* @param alpha * @param alpha
* @param randomize * @param randomize
*/ */
public static void rotateAllAxes(double[] vect, double alpha, boolean randomize) { public static void rotateAllAxes(double[] vect, double alpha,
boolean randomize) {
for (int i = 0; i < vect.length - 1; i++) { for (int i = 0; i < vect.length - 1; i++) {
for (int j = i + 1; j < vect.length; j++) { for (int j = i + 1; j < vect.length; j++) {
if (randomize) rotate(vect, RNG.randomDouble(-alpha,alpha), i, j); if (randomize)
else rotate(vect, alpha, i, j); rotate(vect, RNG.randomDouble(-alpha, alpha), i, j);
else
rotate(vect, alpha, i, j);
} }
} }
} }
@ -811,8 +899,9 @@ public class Mathematics {
} }
/** /**
* Scale a range by the given factor, meaning that the interval in each dimension is * Scale a range by the given factor, meaning that the interval in each
* extended (fact>1) or reduced (fact<1) by the defined ratio around the center. * dimension is extended (fact>1) or reduced (fact<1) by the defined ratio
* around the center.
* *
* @param rangeScaleFact * @param rangeScaleFact
* @param range * @param range
@ -820,7 +909,10 @@ public class Mathematics {
public static void scaleRange(double rangeScaleFact, double[][] range) { public static void scaleRange(double rangeScaleFact, double[][] range) {
double[] intervalLengths = Mathematics.getAbsRange(range); double[] intervalLengths = Mathematics.getAbsRange(range);
double[] tmpInts = Mathematics.svMult(rangeScaleFact, intervalLengths); double[] tmpInts = Mathematics.svMult(rangeScaleFact, intervalLengths);
Mathematics.vvSub(tmpInts, intervalLengths, tmpInts); // this is what must be added to range interval Mathematics.vvSub(tmpInts, intervalLengths, tmpInts); // this is what
// must be added
// to range
// interval
for (int i = 0; i < range.length; i++) { for (int i = 0; i < range.length; i++) {
range[i][0] -= tmpInts[i] / 2; range[i][0] -= tmpInts[i] / 2;
range[i][1] += tmpInts[i] / 2; range[i][1] += tmpInts[i] / 2;
@ -840,8 +932,8 @@ public class Mathematics {
} }
/** /**
* Shift bounds by a constant value in every dimension. The dists * Shift bounds by a constant value in every dimension. The dists must be of
* must be of dimensions as the range. * dimensions as the range.
* *
* @param range * @param range
* @return * @return
@ -861,10 +953,12 @@ public class Mathematics {
// * @param x1 // * @param x1
// * @param f0 // * @param f0
// * @param f1 // * @param f1
// * @return If an error with the spline occurs, a linear interpolation will be // * @return If an error with the spline occurs, a linear interpolation will
// be
// * returned. // * returned.
// */ // */
// /* public static double splineInterpolation(double x, double x0, double x1, // /* public static double splineInterpolation(double x, double x0, double
// x1,
// double f0, double f1) { // double f0, double f1) {
// try { // try {
// double[] t = { x0, x1 }, f = { f0, f1 }; // double[] t = { x0, x1 }, f = { f0, f1 };
@ -885,15 +979,15 @@ public class Mathematics {
* @param x1 * @param x1
* @param f0 * @param f0
* @param f1 * @param f1
* @return If an error with the spline occurs, a linear interpolation will be * @return If an error with the spline occurs, a linear interpolation will
* returned. * be returned.
*/ */
public static double splineInterpolation(double x, double x0, double x1, public static double splineInterpolation(double x, double x0, double x1,
double f0, double f1) { double f0, double f1) {
try { try {
double[] t = { x0, x1 }, f = { f0, f1 }; double[] t = { x0, x1 }, f = { f0, f1 };
SplineInterpolation spline = new SplineInterpolation(new BasicDataSet(t, SplineInterpolation spline = new SplineInterpolation(
f, 1)); new BasicDataSet(t, f, 1));
return spline.getY(x); return spline.getY(x);
} catch (InterpolationException e) { } catch (InterpolationException e) {
e.printStackTrace(); e.printStackTrace();
@ -915,9 +1009,11 @@ public class Mathematics {
int i, j, m = 0, n = 0; int i, j, m = 0, n = 0;
for (i = 0; i < a.length; i++) { for (i = 0; i < a.length; i++) {
if (i == k) continue; if (i == k)
continue;
for (j = 0; j < a[0].length; j++) { for (j = 0; j < a[0].length; j++) {
if (j == l) continue; if (j == l)
continue;
b[m][n++] = a[i][j]; b[m][n++] = a[i][j];
} }
m++; m++;
@ -930,7 +1026,8 @@ public class Mathematics {
/** /**
* Computes the sum of the elements of an array of doubles. * Computes the sum of the elements of an array of doubles.
* *
* @param doubles the array of double * @param doubles
* the array of double
* @return the sum of the elements * @return the sum of the elements
*/ */
public static double sum(double[] doubles) { public static double sum(double[] doubles) {
@ -944,7 +1041,8 @@ public class Mathematics {
/** /**
* Computes the sum of the elements of an array of integers. * Computes the sum of the elements of an array of integers.
* *
* @param ints the array of integers * @param ints
* the array of integers
* @return the sum of the elements * @return the sum of the elements
*/ */
public static int sum(int[] ints) { public static int sum(int[] ints) {
@ -1012,7 +1110,8 @@ public class Mathematics {
} }
/** /**
* Multiplies (scales) every element of the array v with s returning a new vector. * Multiplies (scales) every element of the array v with s returning a new
* vector.
* *
* @param s * @param s
* a scalar * a scalar
@ -1028,12 +1127,13 @@ public class Mathematics {
return res; return res;
} }
/** /**
* Multiplies (scales) every element of the array v with s in place. * Multiplies (scales) every element of the array v with s in place.
* *
* @param s a scalar * @param s
* @param v an array to be multiplied with s. * a scalar
* @param v
* an array to be multiplied with s.
* @return a scaled array. * @return a scaled array.
*/ */
public static void svMult(double s, double[] v, double[] res) { public static void svMult(double s, double[] v, double[] res) {
@ -1050,7 +1150,8 @@ public class Mathematics {
* @param w * @param w
* @return * @return
*/ */
public static void svvAddScaled(double s, double[] v, double[] w, double[] res) { public static void svvAddScaled(double s, double[] v, double[] w,
double[] res) {
for (int i = 0; i < v.length; i++) { for (int i = 0; i < v.length; i++) {
res[i] = s * v[i] + w[i]; res[i] = s * v[i] + w[i];
} }
@ -1058,6 +1159,7 @@ public class Mathematics {
/** /**
* Add vectors returning a new vector c = a + b; * Add vectors returning a new vector c = a + b;
*
* @param a * @param a
* @param b * @param b
* @return a new vector c = a + b * @return a new vector c = a + b
@ -1097,17 +1199,19 @@ public class Mathematics {
} }
/** /**
* Add vectors in place setting with an offset within the target * Add vectors in place setting with an offset within the target vector,
* vector, meaning that res[resOffs+i]=v1[v1Offs+i]+v2[v2Offs+i] for i in length. * meaning that res[resOffs+i]=v1[v1Offs+i]+v2[v2Offs+i] for i in length.
* *
* @param v1 * @param v1
* @param v2 * @param v2
* @return vector addition * @return vector addition
*/ */
public static void vvAddOffs(double[] v1, int v1Offs, double[] v2, int v2Offs, double[] res, int resOffs, int len) { public static void vvAddOffs(double[] v1, int v1Offs, double[] v2,
for (int i = 0; i < len; i++) int v2Offs, double[] res, int resOffs, int len) {
for (int i = 0; i < len; i++) {
res[resOffs + i] = v1[v1Offs + i] + v2[v2Offs + i]; res[resOffs + i] = v1[v1Offs + i] + v2[v2Offs + i];
} }
}
/** /**
* Scalar product of two vectors returning sum_i (a_i * b_i). * Scalar product of two vectors returning sum_i (a_i * b_i).
@ -1124,7 +1228,6 @@ public class Mathematics {
return result; return result;
} }
/** /**
* Component wise multiplication of vectors: res[i]=u[i]*v[i] * Component wise multiplication of vectors: res[i]=u[i]*v[i]
* *
@ -1166,6 +1269,7 @@ public class Mathematics {
/** /**
* Return a vector of given length containing zeroes. * Return a vector of given length containing zeroes.
*
* @param n * @param n
* @return * @return
*/ */
@ -1187,4 +1291,16 @@ public class Mathematics {
Arrays.fill(ret, d); Arrays.fill(ret, d);
return ret; return ret;
} }
/**
* Scales a vector with the given scalar.
*
* @param scale
* @param vec
*/
public static void scale(double scale, double[] vec) {
for (double d : vec) {
d *= scale;
}
}
} }

2
src/eva2/tools/math/des/DESSolver.java
View File

@ -18,7 +18,7 @@ import java.io.Serializable;
* @version Copyright (c) ZBiT, University of T&uuml;bingen, Germany Compiler: * @version Copyright (c) ZBiT, University of T&uuml;bingen, Germany Compiler:
* JDK 1.6.0 * JDK 1.6.0
* @date Sep 10, 2007 * @date Sep 10, 2007
* * @depend - <call> - DESystem
*/ */
public interface DESSolver extends Serializable { public interface DESSolver extends Serializable {

212
src/eva2/tools/math/des/RKSolver.java
View File

@ -12,6 +12,7 @@ import eva2.tools.math.Mathematics;
* @author Andreas Dr&auml;ger * @author Andreas Dr&auml;ger
* @author Marcel Kronfeld * @author Marcel Kronfeld
* @version 1.0 Status: works, but numerical inaccurate * @version 1.0 Status: works, but numerical inaccurate
* @depend - <call> - Mathematics
*/ */
public class RKSolver implements DESSolver, Serializable { public class RKSolver implements DESSolver, Serializable {
/** /**
@ -89,109 +90,6 @@ public class RKSolver implements DESSolver, Serializable {
return unstableFlag; return unstableFlag;
} }
/**
* @param DES
* @param h
* @param x
* @param Ytemp
* @return
*/
public double[] rkTerm(DESystem DES, double h, double x, double[] Ytemp) {
double[][] K = new double[4][];
K[0] = Mathematics.svMult(h, DES.getValue(x, Ytemp));
K[1] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(
Ytemp, Mathematics.svMult(0.5, K[0]))));
K[2] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(
Ytemp, Mathematics.svMult(0.5, K[1]))));
K[3] = Mathematics.svMult(h, DES.getValue(x + h, Mathematics.vvAdd(
Ytemp, K[2])));
double[] change = Mathematics.svDiv(6, Mathematics.vvAdd(K[0],
Mathematics.vvAdd(Mathematics.svMult(2, K[1]), Mathematics
.vvAdd(Mathematics.svMult(2, K[2]), K[3]))));
for (int k = 0; k < change.length; k++) {
if (Double.isNaN(change[k])) {
unstableFlag = true;
change[k] = 0;
// return result;
}
}
return change;
}
/**
* Linearized code for speed-up (no allocations).
*
* @param DES
* @param h
* @param x
* @param Ytemp
* @return
*/
public void rkTerm2(DESystem DES, double h, double x, double[] Ytemp,
double[] res) {
if (kVals == null) { // "static" vectors which are allocated only once
k0tmp = new double[DES.getDESystemDimension()];
k1tmp = new double[DES.getDESystemDimension()];
k2tmp = new double[DES.getDESystemDimension()];
kVals = new double[4][DES.getDESystemDimension()];
}
// double[][] K = new double[4][];
DES.getValue(x, Ytemp, kVals[0]);
Mathematics.svMult(h, kVals[0], kVals[0]);
// K[0] = svMult(h, DES.getValue(x, Ytemp));
Mathematics.svMult(0.5, kVals[0], k0tmp);
Mathematics.vvAdd(Ytemp, k0tmp, k0tmp);
DES.getValue(x + h / 2, k0tmp, kVals[1]);
Mathematics.svMult(h, kVals[1], kVals[1]);
// K[1] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5,
// K[0]))));
Mathematics.svMult(0.5, kVals[1], k1tmp);
Mathematics.vvAdd(Ytemp, k1tmp, k1tmp);
DES.getValue(x + h / 2, k1tmp, kVals[2]);
Mathematics.svMult(h, kVals[2], kVals[2]);
// K[2] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5,
// K[1]))));
Mathematics.vvAdd(Ytemp, kVals[2], k2tmp);
DES.getValue(x + h, k2tmp, k1tmp);
Mathematics.svMult(h, k1tmp, kVals[3]);
// K[3] = svMult(h, DES.getValue(x + h, vvAdd(Ytemp, K[2])));
Mathematics.svMult(2, kVals[2], k0tmp);
Mathematics.vvAdd(k0tmp, kVals[3], k0tmp);
Mathematics.svMult(2, kVals[1], k1tmp);
Mathematics.vvAdd(k1tmp, k0tmp, k2tmp);
Mathematics.vvAdd(kVals[0], k2tmp, k1tmp);
Mathematics.svDiv(6, k1tmp, res);
// double[] change = svDiv(6, vvAdd(K[0], vvAdd(svMult(2, K[1]),
// vvAdd(svMult(2, K[2]), K[3]))));
// for (int i=0; i<res.length; i++) {
// double diff = Math.abs(res[i]-change[i]);
// if (diff > 0.00000001) System.out.println("!!! ");
// }
// double[] change = svdiv(6, vvadd(kVals[0], vvadd(svmult(2, kVals[1]),
// vvadd(svmult(2, kVals[2]), kVals[3]))));
for (int k = 0; k < res.length; k++) {
if (Double.isNaN(res[k])) {
unstableFlag = true;
res[k] = 0;
// return result;
}
}
}
/** /**
* @param stepSize * @param stepSize
*/ */
@ -254,8 +152,9 @@ public class RKSolver implements DESSolver, Serializable {
return solveAtTimePoints(DES, initialValues, timePoints, true); return solveAtTimePoints(DES, initialValues, timePoints, true);
} }
/** /*
* * (non-Javadoc)
* @see eva2.tools.math.des.DESSolver#solveAtTimePointsWithInitialConditions(eva2.tools.math.des.DESystem, double[][], double[])
*/ */
public double[][] solveAtTimePointsWithInitialConditions(DESystem DES, public double[][] solveAtTimePointsWithInitialConditions(DESystem DES,
double[][] initConditions, double[] timePoints) { double[][] initConditions, double[] timePoints) {
@ -319,6 +218,109 @@ public class RKSolver implements DESSolver, Serializable {
return solveByStepSize(DES, initialValues, timeBegin, timeEnd, true); return solveByStepSize(DES, initialValues, timeBegin, timeEnd, true);
} }
/**
* @param DES
* @param h
* @param x
* @param Ytemp
* @return
*/
private double[] rkTerm(DESystem DES, double h, double x, double[] Ytemp) {
double[][] K = new double[4][];
K[0] = Mathematics.svMult(h, DES.getValue(x, Ytemp));
K[1] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(
Ytemp, Mathematics.svMult(0.5, K[0]))));
K[2] = Mathematics.svMult(h, DES.getValue(x + h / 2, Mathematics.vvAdd(
Ytemp, Mathematics.svMult(0.5, K[1]))));
K[3] = Mathematics.svMult(h, DES.getValue(x + h, Mathematics.vvAdd(
Ytemp, K[2])));
double[] change = Mathematics.svDiv(6, Mathematics.vvAdd(K[0],
Mathematics.vvAdd(Mathematics.svMult(2, K[1]), Mathematics
.vvAdd(Mathematics.svMult(2, K[2]), K[3]))));
for (int k = 0; k < change.length; k++) {
if (Double.isNaN(change[k])) {
unstableFlag = true;
change[k] = 0;
// return result;
}
}
return change;
}
/**
* Linearized code for speed-up (no allocations).
*
* @param DES
* @param h
* @param x
* @param Ytemp
* @return
*/
private void rkTerm2(DESystem DES, double h, double x, double[] Ytemp,
double[] res) {
if (kVals == null) { // "static" vectors which are allocated only once
k0tmp = new double[DES.getDESystemDimension()];
k1tmp = new double[DES.getDESystemDimension()];
k2tmp = new double[DES.getDESystemDimension()];
kVals = new double[4][DES.getDESystemDimension()];
}
// double[][] K = new double[4][];
DES.getValue(x, Ytemp, kVals[0]);
Mathematics.svMult(h, kVals[0], kVals[0]);
// K[0] = svMult(h, DES.getValue(x, Ytemp));
Mathematics.svMult(0.5, kVals[0], k0tmp);
Mathematics.vvAdd(Ytemp, k0tmp, k0tmp);
DES.getValue(x + h / 2, k0tmp, kVals[1]);
Mathematics.svMult(h, kVals[1], kVals[1]);
// K[1] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5,
// K[0]))));
Mathematics.svMult(0.5, kVals[1], k1tmp);
Mathematics.vvAdd(Ytemp, k1tmp, k1tmp);
DES.getValue(x + h / 2, k1tmp, kVals[2]);
Mathematics.svMult(h, kVals[2], kVals[2]);
// K[2] = svMult(h, DES.getValue(x + h / 2, vvAdd(Ytemp, svMult(0.5,
// K[1]))));
Mathematics.vvAdd(Ytemp, kVals[2], k2tmp);
DES.getValue(x + h, k2tmp, k1tmp);
Mathematics.svMult(h, k1tmp, kVals[3]);
// K[3] = svMult(h, DES.getValue(x + h, vvAdd(Ytemp, K[2])));
Mathematics.svMult(2, kVals[2], k0tmp);
Mathematics.vvAdd(k0tmp, kVals[3], k0tmp);
Mathematics.svMult(2, kVals[1], k1tmp);
Mathematics.vvAdd(k1tmp, k0tmp, k2tmp);
Mathematics.vvAdd(kVals[0], k2tmp, k1tmp);
Mathematics.svDiv(6, k1tmp, res);
// double[] change = svDiv(6, vvAdd(K[0], vvAdd(svMult(2, K[1]),
// vvAdd(svMult(2, K[2]), K[3]))));
// for (int i=0; i<res.length; i++) {
// double diff = Math.abs(res[i]-change[i]);
// if (diff > 0.00000001) System.out.println("!!! ");
// }
// double[] change = svdiv(6, vvadd(kVals[0], vvadd(svmult(2, kVals[1]),
// vvadd(svmult(2, kVals[2]), kVals[3]))));
for (int k = 0; k < res.length; k++) {
if (Double.isNaN(res[k])) {
unstableFlag = true;
res[k] = 0;
// return result;
}
}
}
/** /**
* When set to <code>TRUE</code>, <code>includeTimes</code> will make the * When set to <code>TRUE</code>, <code>includeTimes</code> will make the
* solver to return a matrix with the first column containing the times. By * solver to return a matrix with the first column containing the times. By