+ For a symmetric, positive definite matrix A, the Cholesky decomposition
+ is an lower triangular matrix L so that A = L*L'.
+
+ If the matrix is not symmetric or positive definite, the constructor
+ returns a partial decomposition and sets an internal flag that may
+ be queried by the isSPD() method.
+ */
+
+public class CholeskyDecomposition implements java.io.Serializable {
+
+/* ------------------------
+ Class variables
+ * ------------------------ */
+
+ /** Array for internal storage of decomposition.
+ @serial internal array storage.
+ */
+ private double[][] L;
+
+ /** Row and column dimension (square matrix).
+ @serial matrix dimension.
+ */
+ private int n;
+
+ /** Symmetric and positive definite flag.
+ @serial is symmetric and positive definite flag.
+ */
+ private boolean isspd;
+
+/* ------------------------
+ Constructor
+ * ------------------------ */
+
+ /** Cholesky algorithm for symmetric and positive definite matrix.
+ @param A Square, symmetric matrix.
+ @return Structure to access L and isspd flag.
+ */
+
+ public CholeskyDecomposition (Matrix Arg) {
+ // Initialize.
+ double[][] A = Arg.getArray();
+ n = Arg.getRowDimension();
+ L = new double[n][n];
+ isspd = (Arg.getColumnDimension() == n);
+ // Main loop.
+ for (int j = 0; j < n; j++) {
+ double[] Lrowj = L[j];
+ double d = 0.0;
+ for (int k = 0; k < j; k++) {
+ double[] Lrowk = L[k];
+ double s = 0.0;
+ for (int i = 0; i < k; i++) {
+ s += Lrowk[i]*Lrowj[i];
+ }
+ Lrowj[k] = s = (A[j][k] - s)/L[k][k];
+ d = d + s*s;
+ isspd = isspd & (A[k][j] == A[j][k]);
+ }
+ d = A[j][j] - d;
+ isspd = isspd & (d > 0.0);
+ L[j][j] = Math.sqrt(Math.max(d,0.0));
+ for (int k = j+1; k < n; k++) {
+ L[j][k] = 0.0;
+ }
+ }
+ }
+
+
+
+// \** Right Triangular Cholesky Decomposition.
+//
+// For a symmetric, positive definite matrix A, the Right Cholesky
+// decomposition is an upper triangular matrix R so that A = R'*R.
+// This constructor computes R with the Fortran inspired column oriented
+// algorithm used in LINPACK and MATLAB. In Java, we suspect a row oriented,
+// lower triangular decomposition is faster. We have temporarily included
+// this constructor here until timing experiments confirm this suspicion.
+// *\
+
+
+ private transient double[][] R;
+
+
+
+ public CholeskyDecomposition (Matrix Arg, int rightflag) {
+ // Initialize.
+ double[][] A = Arg.getArray();
+ n = Arg.getColumnDimension();
+ R = new double[n][n];
+ isspd = (Arg.getColumnDimension() == n);
+ // Main loop.
+ for (int j = 0; j < n; j++) {
+ double d = 0.0;
+ for (int k = 0; k < j; k++) {
+ double s = A[k][j];
+ for (int i = 0; i < k; i++) {
+ s = s - R[i][k]*R[i][j];
+ }
+ R[k][j] = s = s/R[k][k];
+ d = d + s*s;
+ isspd = isspd & (A[k][j] == A[j][k]);
+ }
+ d = A[j][j] - d;
+ isspd = isspd & (d > 0.0);
+ R[j][j] = Math.sqrt(Math.max(d,0.0));
+ for (int k = j+1; k < n; k++) {
+ R[k][j] = 0.0;
+ }
+ }
+ }
+
+
+ public Matrix getR () {
+ return new Matrix(R,n,n);
+ }
+
+
+ /** Is the matrix symmetric and positive definite?
+ @return true if A is symmetric and positive definite.
+ */
+
+ public boolean isSPD () {
+ return isspd;
+ }
+
+ /** Return triangular factor.
+ @return L
+ */
+
+ public Matrix getL () {
+ return new Matrix(L,n,n);
+ }
+
+ /** Solve A*X = B
+ @param B A Matrix with as many rows as A and any number of columns.
+ @return X so that L*L'*X = B
+ @exception IllegalArgumentException Matrix row dimensions must agree.
+ @exception RuntimeException Matrix is not symmetric positive definite.
+ */
+
+ public Matrix solve (Matrix B) {
+ if (B.getRowDimension() != n) {
+ throw new IllegalArgumentException("Matrix row dimensions must agree.");
+ }
+ if (!isspd) {
+ throw new RuntimeException("Matrix is not symmetric positive definite.");
+ }
+
+ // Copy right hand side.
+ double[][] X = B.getArrayCopy();
+ int nx = B.getColumnDimension();
+
+ // Solve L*Y = B;
+ for (int k = 0; k < n; k++) {
+ for (int i = k+1; i < n; i++) {
+ for (int j = 0; j < nx; j++) {
+ X[i][j] -= X[k][j]*L[i][k];
+ }
+ }
+ for (int j = 0; j < nx; j++) {
+ X[k][j] /= L[k][k];
+ }
+ }
+
+ // Solve L'*X = Y;
+ for (int k = n-1; k >= 0; k--) {
+ for (int j = 0; j < nx; j++) {
+ X[k][j] /= L[k][k];
+ }
+ for (int i = 0; i < k; i++) {
+ for (int j = 0; j < nx; j++) {
+ X[i][j] -= X[k][j]*L[k][i];
+ }
+ }
+ }
+ return new Matrix(X,n,nx);
+ }
+}
diff --git a/src/wsi/ra/math/Jama/EigenvalueDecomposition.java b/src/wsi/ra/math/Jama/EigenvalueDecomposition.java
new file mode 100644
index 00000000..e351c4b1
--- /dev/null
+++ b/src/wsi/ra/math/Jama/EigenvalueDecomposition.java
@@ -0,0 +1,956 @@
+package wsi.ra.math.Jama;
+import wsi.ra.math.Jama.util.Maths;
+
+
+/** Eigenvalues and eigenvectors of a real matrix.
+
+ If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
+ diagonal and the eigenvector matrix V is orthogonal.
+ I.e. A = V.times(D.times(V.transpose())) and
+ V.times(V.transpose()) equals the identiCty matrix.
+
+ If A is not symmetric, then the eigenvalue matrix D is block diagonal
+ with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
+ lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
+ columns of V represent the eigenvectors in the sense that A*V = V*D,
+ i.e. A.times(V) equals V.times(D). The matrix V may be badly
+ conditioned, or even singular, so the validity of the equation
+ A = V*D*inverse(V) depends upon V.cond().
+**/
+
+public class EigenvalueDecomposition implements java.io.Serializable {
+
+/* ------------------------
+ Class variables
+ * ------------------------ */
+
+ /** Row and column dimension (square matrix).
+ @serial matrix dimension.
+ */
+ private int n;
+
+ /** Symmetry flag.
+ @serial internal symmetry flag.
+ */
+ private boolean issymmetric;
+
+ /** Arrays for internal storage of eigenvalues.
+ @serial internal storage of eigenvalues.
+ */
+ private double[] d, e;
+
+ /** Array for internal storage of eigenvectors.
+ @serial internal storage of eigenvectors.
+ */
+ private double[][] V;
+
+ /** Array for internal storage of nonsymmetric Hessenberg form.
+ @serial internal storage of nonsymmetric Hessenberg form.
+ */
+ private double[][] H;
+
+ /** Working storage for nonsymmetric algorithm.
+ @serial working storage for nonsymmetric algorithm.
+ */
+ private double[] ort;
+
+/* ------------------------
+ Private Methods
+ * ------------------------ */
+
+ // Symmetric Householder reduction to tridiagonal form.
+
+ private void tred2 () {
+
+ // This is derived from the Algol procedures tred2 by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ for (int j = 0; j < n; j++) {
+ d[j] = V[n-1][j];
+ }
+
+ // Householder reduction to tridiagonal form.
+
+ for (int i = n-1; i > 0; i--) {
+
+ // Scale to avoid under/overflow.
+
+ double scale = 0.0;
+ double h = 0.0;
+ for (int k = 0; k < i; k++) {
+ scale = scale + Math.abs(d[k]);
+ }
+ if (scale == 0.0) {
+ e[i] = d[i-1];
+ for (int j = 0; j < i; j++) {
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ V[j][i] = 0.0;
+ }
+ } else {
+
+ // Generate Householder vector.
+
+ for (int k = 0; k < i; k++) {
+ d[k] /= scale;
+ h += d[k] * d[k];
+ }
+ double f = d[i-1];
+ double g = Math.sqrt(h);
+ if (f > 0) {
+ g = -g;
+ }
+ e[i] = scale * g;
+ h = h - f * g;
+ d[i-1] = f - g;
+ for (int j = 0; j < i; j++) {
+ e[j] = 0.0;
+ }
+
+ // Apply similarity transformation to remaining columns.
+
+ for (int j = 0; j < i; j++) {
+ f = d[j];
+ V[j][i] = f;
+ g = e[j] + V[j][j] * f;
+ for (int k = j+1; k <= i-1; k++) {
+ g += V[k][j] * d[k];
+ e[k] += V[k][j] * f;
+ }
+ e[j] = g;
+ }
+ f = 0.0;
+ for (int j = 0; j < i; j++) {
+ e[j] /= h;
+ f += e[j] * d[j];
+ }
+ double hh = f / (h + h);
+ for (int j = 0; j < i; j++) {
+ e[j] -= hh * d[j];
+ }
+ for (int j = 0; j < i; j++) {
+ f = d[j];
+ g = e[j];
+ for (int k = j; k <= i-1; k++) {
+ V[k][j] -= (f * e[k] + g * d[k]);
+ }
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ }
+ }
+ d[i] = h;
+ }
+
+ // Accumulate transformations.
+
+ for (int i = 0; i < n-1; i++) {
+ V[n-1][i] = V[i][i];
+ V[i][i] = 1.0;
+ double h = d[i+1];
+ if (h != 0.0) {
+ for (int k = 0; k <= i; k++) {
+ d[k] = V[k][i+1] / h;
+ }
+ for (int j = 0; j <= i; j++) {
+ double g = 0.0;
+ for (int k = 0; k <= i; k++) {
+ g += V[k][i+1] * V[k][j];
+ }
+ for (int k = 0; k <= i; k++) {
+ V[k][j] -= g * d[k];
+ }
+ }
+ }
+ for (int k = 0; k <= i; k++) {
+ V[k][i+1] = 0.0;
+ }
+ }
+ for (int j = 0; j < n; j++) {
+ d[j] = V[n-1][j];
+ V[n-1][j] = 0.0;
+ }
+ V[n-1][n-1] = 1.0;
+ e[0] = 0.0;
+ }
+
+ // Symmetric tridiagonal QL algorithm.
+
+ private void tql2 () {
+
+ // This is derived from the Algol procedures tql2, by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ for (int i = 1; i < n; i++) {
+ e[i-1] = e[i];
+ }
+ e[n-1] = 0.0;
+
+ double f = 0.0;
+ double tst1 = 0.0;
+ double eps = Math.pow(2.0,-52.0);
+ for (int l = 0; l < n; l++) {
+
+ // Find small subdiagonal element
+
+ tst1 = Math.max(tst1,Math.abs(d[l]) + Math.abs(e[l]));
+ int m = l;
+ while (m < n) {
+ if (Math.abs(e[m]) <= eps*tst1) {
+ break;
+ }
+ m++;
+ }
+
+ // If m == l, d[l] is an eigenvalue,
+ // otherwise, iterate.
+
+ if (m > l) {
+ int iter = 0;
+ do {
+ iter = iter + 1; // (Could check iteration count here.)
+
+ // Compute implicit shift
+
+ double g = d[l];
+ double p = (d[l+1] - g) / (2.0 * e[l]);
+ double r = Maths.hypot(p,1.0);
+ if (p < 0) {
+ r = -r;
+ }
+ d[l] = e[l] / (p + r);
+ d[l+1] = e[l] * (p + r);
+ double dl1 = d[l+1];
+ double h = g - d[l];
+ for (int i = l+2; i < n; i++) {
+ d[i] -= h;
+ }
+ f = f + h;
+
+ // Implicit QL transformation.
+
+ p = d[m];
+ double c = 1.0;
+ double c2 = c;
+ double c3 = c;
+ double el1 = e[l+1];
+ double s = 0.0;
+ double s2 = 0.0;
+ for (int i = m-1; i >= l; i--) {
+ c3 = c2;
+ c2 = c;
+ s2 = s;
+ g = c * e[i];
+ h = c * p;
+ r = Maths.hypot(p,e[i]);
+ e[i+1] = s * r;
+ s = e[i] / r;
+ c = p / r;
+ p = c * d[i] - s * g;
+ d[i+1] = h + s * (c * g + s * d[i]);
+
+ // Accumulate transformation.
+
+ for (int k = 0; k < n; k++) {
+ h = V[k][i+1];
+ V[k][i+1] = s * V[k][i] + c * h;
+ V[k][i] = c * V[k][i] - s * h;
+ }
+ }
+ p = -s * s2 * c3 * el1 * e[l] / dl1;
+ e[l] = s * p;
+ d[l] = c * p;
+
+ // Check for convergence.
+
+ } while (Math.abs(e[l]) > eps*tst1);
+ }
+ d[l] = d[l] + f;
+ e[l] = 0.0;
+ }
+
+ // Sort eigenvalues and corresponding vectors.
+
+ for (int i = 0; i < n-1; i++) {
+ int k = i;
+ double p = d[i];
+ for (int j = i+1; j < n; j++) {
+ if (d[j] < p) {
+ k = j;
+ p = d[j];
+ }
+ }
+ if (k != i) {
+ d[k] = d[i];
+ d[i] = p;
+ for (int j = 0; j < n; j++) {
+ p = V[j][i];
+ V[j][i] = V[j][k];
+ V[j][k] = p;
+ }
+ }
+ }
+ }
+
+ // Nonsymmetric reduction to Hessenberg form.
+
+ private void orthes () {
+
+ // This is derived from the Algol procedures orthes and ortran,
+ // by Martin and Wilkinson, Handbook for Auto. Comp.,
+ // Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutines in EISPACK.
+
+ int low = 0;
+ int high = n-1;
+
+ for (int m = low+1; m <= high-1; m++) {
+
+ // Scale column.
+
+ double scale = 0.0;
+ for (int i = m; i <= high; i++) {
+ scale = scale + Math.abs(H[i][m-1]);
+ }
+ if (scale != 0.0) {
+
+ // Compute Householder transformation.
+
+ double h = 0.0;
+ for (int i = high; i >= m; i--) {
+ ort[i] = H[i][m-1]/scale;
+ h += ort[i] * ort[i];
+ }
+ double g = Math.sqrt(h);
+ if (ort[m] > 0) {
+ g = -g;
+ }
+ h = h - ort[m] * g;
+ ort[m] = ort[m] - g;
+
+ // Apply Householder similarity transformation
+ // H = (I-u*u'/h)*H*(I-u*u')/h)
+
+ for (int j = m; j < n; j++) {
+ double f = 0.0;
+ for (int i = high; i >= m; i--) {
+ f += ort[i]*H[i][j];
+ }
+ f = f/h;
+ for (int i = m; i <= high; i++) {
+ H[i][j] -= f*ort[i];
+ }
+ }
+
+ for (int i = 0; i <= high; i++) {
+ double f = 0.0;
+ for (int j = high; j >= m; j--) {
+ f += ort[j]*H[i][j];
+ }
+ f = f/h;
+ for (int j = m; j <= high; j++) {
+ H[i][j] -= f*ort[j];
+ }
+ }
+ ort[m] = scale*ort[m];
+ H[m][m-1] = scale*g;
+ }
+ }
+
+ // Accumulate transformations (Algol's ortran).
+
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ V[i][j] = (i == j ? 1.0 : 0.0);
+ }
+ }
+
+ for (int m = high-1; m >= low+1; m--) {
+ if (H[m][m-1] != 0.0) {
+ for (int i = m+1; i <= high; i++) {
+ ort[i] = H[i][m-1];
+ }
+ for (int j = m; j <= high; j++) {
+ double g = 0.0;
+ for (int i = m; i <= high; i++) {
+ g += ort[i] * V[i][j];
+ }
+ // Double division avoids possible underflow
+ g = (g / ort[m]) / H[m][m-1];
+ for (int i = m; i <= high; i++) {
+ V[i][j] += g * ort[i];
+ }
+ }
+ }
+ }
+ }
+
+
+ // Complex scalar division.
+
+ private transient double cdivr, cdivi;
+ private void cdiv(double xr, double xi, double yr, double yi) {
+ double r,d;
+ if (Math.abs(yr) > Math.abs(yi)) {
+ r = yi/yr;
+ d = yr + r*yi;
+ cdivr = (xr + r*xi)/d;
+ cdivi = (xi - r*xr)/d;
+ } else {
+ r = yr/yi;
+ d = yi + r*yr;
+ cdivr = (r*xr + xi)/d;
+ cdivi = (r*xi - xr)/d;
+ }
+ }
+
+
+ // Nonsymmetric reduction from Hessenberg to real Schur form.
+
+ private void hqr2 () {
+
+ // This is derived from the Algol procedure hqr2,
+ // by Martin and Wilkinson, Handbook for Auto. Comp.,
+ // Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ // Initialize
+
+ int nn = this.n;
+ int n = nn-1;
+ int low = 0;
+ int high = nn-1;
+ double eps = Math.pow(2.0,-52.0);
+ double exshift = 0.0;
+ double p=0,q=0,r=0,s=0,z=0,t,w,x,y;
+
+ // Store roots isolated by balanc and compute matrix norm
+
+ double norm = 0.0;
+ for (int i = 0; i < nn; i++) {
+ if (i < low | i > high) {
+ d[i] = H[i][i];
+ e[i] = 0.0;
+ }
+ for (int j = Math.max(i-1,0); j < nn; j++) {
+ norm = norm + Math.abs(H[i][j]);
+ }
+ }
+
+ // Outer loop over eigenvalue index
+
+ int iter = 0;
+ while (n >= low) {
+
+ // Look for single small sub-diagonal element
+
+ int l = n;
+ while (l > low) {
+ s = Math.abs(H[l-1][l-1]) + Math.abs(H[l][l]);
+ if (s == 0.0) {
+ s = norm;
+ }
+ if (Math.abs(H[l][l-1]) < eps * s) {
+ break;
+ }
+ l--;
+ }
+
+ // Check for convergence
+ // One root found
+
+ if (l == n) {
+ H[n][n] = H[n][n] + exshift;
+ d[n] = H[n][n];
+ e[n] = 0.0;
+ n--;
+ iter = 0;
+
+ // Two roots found
+
+ } else if (l == n-1) {
+ w = H[n][n-1] * H[n-1][n];
+ p = (H[n-1][n-1] - H[n][n]) / 2.0;
+ q = p * p + w;
+ z = Math.sqrt(Math.abs(q));
+ H[n][n] = H[n][n] + exshift;
+ H[n-1][n-1] = H[n-1][n-1] + exshift;
+ x = H[n][n];
+
+ // Real pair
+
+ if (q >= 0) {
+ if (p >= 0) {
+ z = p + z;
+ } else {
+ z = p - z;
+ }
+ d[n-1] = x + z;
+ d[n] = d[n-1];
+ if (z != 0.0) {
+ d[n] = x - w / z;
+ }
+ e[n-1] = 0.0;
+ e[n] = 0.0;
+ x = H[n][n-1];
+ s = Math.abs(x) + Math.abs(z);
+ p = x / s;
+ q = z / s;
+ r = Math.sqrt(p * p+q * q);
+ p = p / r;
+ q = q / r;
+
+ // Row modification
+
+ for (int j = n-1; j < nn; j++) {
+ z = H[n-1][j];
+ H[n-1][j] = q * z + p * H[n][j];
+ H[n][j] = q * H[n][j] - p * z;
+ }
+
+ // Column modification
+
+ for (int i = 0; i <= n; i++) {
+ z = H[i][n-1];
+ H[i][n-1] = q * z + p * H[i][n];
+ H[i][n] = q * H[i][n] - p * z;
+ }
+
+ // Accumulate transformations
+
+ for (int i = low; i <= high; i++) {
+ z = V[i][n-1];
+ V[i][n-1] = q * z + p * V[i][n];
+ V[i][n] = q * V[i][n] - p * z;
+ }
+
+ // Complex pair
+
+ } else {
+ d[n-1] = x + p;
+ d[n] = x + p;
+ e[n-1] = z;
+ e[n] = -z;
+ }
+ n = n - 2;
+ iter = 0;
+
+ // No convergence yet
+
+ } else {
+
+ // Form shift
+
+ x = H[n][n];
+ y = 0.0;
+ w = 0.0;
+ if (l < n) {
+ y = H[n-1][n-1];
+ w = H[n][n-1] * H[n-1][n];
+ }
+
+ // Wilkinson's original ad hoc shift
+
+ if (iter == 10) {
+ exshift += x;
+ for (int i = low; i <= n; i++) {
+ H[i][i] -= x;
+ }
+ s = Math.abs(H[n][n-1]) + Math.abs(H[n-1][n-2]);
+ x = y = 0.75 * s;
+ w = -0.4375 * s * s;
+ }
+
+ // MATLAB's new ad hoc shift
+
+ if (iter == 30) {
+ s = (y - x) / 2.0;
+ s = s * s + w;
+ if (s > 0) {
+ s = Math.sqrt(s);
+ if (y < x) {
+ s = -s;
+ }
+ s = x - w / ((y - x) / 2.0 + s);
+ for (int i = low; i <= n; i++) {
+ H[i][i] -= s;
+ }
+ exshift += s;
+ x = y = w = 0.964;
+ }
+ }
+
+ iter = iter + 1; // (Could check iteration count here.)
+
+ // Look for two consecutive small sub-diagonal elements
+
+ int m = n-2;
+ while (m >= l) {
+ z = H[m][m];
+ r = x - z;
+ s = y - z;
+ p = (r * s - w) / H[m+1][m] + H[m][m+1];
+ q = H[m+1][m+1] - z - r - s;
+ r = H[m+2][m+1];
+ s = Math.abs(p) + Math.abs(q) + Math.abs(r);
+ p = p / s;
+ q = q / s;
+ r = r / s;
+ if (m == l) {
+ break;
+ }
+ if (Math.abs(H[m][m-1]) * (Math.abs(q) + Math.abs(r)) <
+ eps * (Math.abs(p) * (Math.abs(H[m-1][m-1]) + Math.abs(z) +
+ Math.abs(H[m+1][m+1])))) {
+ break;
+ }
+ m--;
+ }
+
+ for (int i = m+2; i <= n; i++) {
+ H[i][i-2] = 0.0;
+ if (i > m+2) {
+ H[i][i-3] = 0.0;
+ }
+ }
+
+ // Double QR step involving rows l:n and columns m:n
+
+ for (int k = m; k <= n-1; k++) {
+ boolean notlast = (k != n-1);
+ if (k != m) {
+ p = H[k][k-1];
+ q = H[k+1][k-1];
+ r = (notlast ? H[k+2][k-1] : 0.0);
+ x = Math.abs(p) + Math.abs(q) + Math.abs(r);
+ if (x != 0.0) {
+ p = p / x;
+ q = q / x;
+ r = r / x;
+ }
+ }
+ if (x == 0.0) {
+ break;
+ }
+ s = Math.sqrt(p * p + q * q + r * r);
+ if (p < 0) {
+ s = -s;
+ }
+ if (s != 0) {
+ if (k != m) {
+ H[k][k-1] = -s * x;
+ } else if (l != m) {
+ H[k][k-1] = -H[k][k-1];
+ }
+ p = p + s;
+ x = p / s;
+ y = q / s;
+ z = r / s;
+ q = q / p;
+ r = r / p;
+
+ // Row modification
+
+ for (int j = k; j < nn; j++) {
+ p = H[k][j] + q * H[k+1][j];
+ if (notlast) {
+ p = p + r * H[k+2][j];
+ H[k+2][j] = H[k+2][j] - p * z;
+ }
+ H[k][j] = H[k][j] - p * x;
+ H[k+1][j] = H[k+1][j] - p * y;
+ }
+
+ // Column modification
+
+ for (int i = 0; i <= Math.min(n,k+3); i++) {
+ p = x * H[i][k] + y * H[i][k+1];
+ if (notlast) {
+ p = p + z * H[i][k+2];
+ H[i][k+2] = H[i][k+2] - p * r;
+ }
+ H[i][k] = H[i][k] - p;
+ H[i][k+1] = H[i][k+1] - p * q;
+ }
+
+ // Accumulate transformations
+
+ for (int i = low; i <= high; i++) {
+ p = x * V[i][k] + y * V[i][k+1];
+ if (notlast) {
+ p = p + z * V[i][k+2];
+ V[i][k+2] = V[i][k+2] - p * r;
+ }
+ V[i][k] = V[i][k] - p;
+ V[i][k+1] = V[i][k+1] - p * q;
+ }
+ } // (s != 0)
+ } // k loop
+ } // check convergence
+ } // while (n >= low)
+
+ // Backsubstitute to find vectors of upper triangular form
+
+ if (norm == 0.0) {
+ return;
+ }
+
+ for (n = nn-1; n >= 0; n--) {
+ p = d[n];
+ q = e[n];
+
+ // Real vector
+
+ if (q == 0) {
+ int l = n;
+ H[n][n] = 1.0;
+ for (int i = n-1; i >= 0; i--) {
+ w = H[i][i] - p;
+ r = 0.0;
+ for (int j = l; j <= n; j++) {
+ r = r + H[i][j] * H[j][n];
+ }
+ if (e[i] < 0.0) {
+ z = w;
+ s = r;
+ } else {
+ l = i;
+ if (e[i] == 0.0) {
+ if (w != 0.0) {
+ H[i][n] = -r / w;
+ } else {
+ H[i][n] = -r / (eps * norm);
+ }
+
+ // Solve real equations
+
+ } else {
+ x = H[i][i+1];
+ y = H[i+1][i];
+ q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
+ t = (x * s - z * r) / q;
+ H[i][n] = t;
+ if (Math.abs(x) > Math.abs(z)) {
+ H[i+1][n] = (-r - w * t) / x;
+ } else {
+ H[i+1][n] = (-s - y * t) / z;
+ }
+ }
+
+ // Overflow control
+
+ t = Math.abs(H[i][n]);
+ if ((eps * t) * t > 1) {
+ for (int j = i; j <= n; j++) {
+ H[j][n] = H[j][n] / t;
+ }
+ }
+ }
+ }
+
+ // Complex vector
+
+ } else if (q < 0) {
+ int l = n-1;
+
+ // Last vector component imaginary so matrix is triangular
+
+ if (Math.abs(H[n][n-1]) > Math.abs(H[n-1][n])) {
+ H[n-1][n-1] = q / H[n][n-1];
+ H[n-1][n] = -(H[n][n] - p) / H[n][n-1];
+ } else {
+ cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q);
+ H[n-1][n-1] = cdivr;
+ H[n-1][n] = cdivi;
+ }
+ H[n][n-1] = 0.0;
+ H[n][n] = 1.0;
+ for (int i = n-2; i >= 0; i--) {
+ double ra,sa,vr,vi;
+ ra = 0.0;
+ sa = 0.0;
+ for (int j = l; j <= n; j++) {
+ ra = ra + H[i][j] * H[j][n-1];
+ sa = sa + H[i][j] * H[j][n];
+ }
+ w = H[i][i] - p;
+
+ if (e[i] < 0.0) {
+ z = w;
+ r = ra;
+ s = sa;
+ } else {
+ l = i;
+ if (e[i] == 0) {
+ cdiv(-ra,-sa,w,q);
+ H[i][n-1] = cdivr;
+ H[i][n] = cdivi;
+ } else {
+
+ // Solve complex equations
+
+ x = H[i][i+1];
+ y = H[i+1][i];
+ vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
+ vi = (d[i] - p) * 2.0 * q;
+ if (vr == 0.0 & vi == 0.0) {
+ vr = eps * norm * (Math.abs(w) + Math.abs(q) +
+ Math.abs(x) + Math.abs(y) + Math.abs(z));
+ }
+ cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
+ H[i][n-1] = cdivr;
+ H[i][n] = cdivi;
+ if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
+ H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x;
+ H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x;
+ } else {
+ cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q);
+ H[i+1][n-1] = cdivr;
+ H[i+1][n] = cdivi;
+ }
+ }
+
+ // Overflow control
+
+ t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n]));
+ if ((eps * t) * t > 1) {
+ for (int j = i; j <= n; j++) {
+ H[j][n-1] = H[j][n-1] / t;
+ H[j][n] = H[j][n] / t;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Vectors of isolated roots
+
+ for (int i = 0; i < nn; i++) {
+ if (i < low | i > high) {
+ for (int j = i; j < nn; j++) {
+ V[i][j] = H[i][j];
+ }
+ }
+ }
+
+ // Back transformation to get eigenvectors of original matrix
+
+ for (int j = nn-1; j >= low; j--) {
+ for (int i = low; i <= high; i++) {
+ z = 0.0;
+ for (int k = low; k <= Math.min(j,high); k++) {
+ z = z + V[i][k] * H[k][j];
+ }
+ V[i][j] = z;
+ }
+ }
+ }
+
+
+/* ------------------------
+ Constructor
+ * ------------------------ */
+
+ /** Check for symmetry, then construct the eigenvalue decomposition
+ @param A Square matrix
+ @return Structure to access D and V.
+ */
+
+ public EigenvalueDecomposition (Matrix Arg) {
+ double[][] A = Arg.getArray();
+ n = Arg.getColumnDimension();
+ V = new double[n][n];
+ d = new double[n];
+ e = new double[n];
+
+ issymmetric = true;
+ for (int j = 0; (j < n) & issymmetric; j++) {
+ for (int i = 0; (i < n) & issymmetric; i++) {
+ issymmetric = (A[i][j] == A[j][i]);
+ }
+ }
+
+ if (issymmetric) {
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ V[i][j] = A[i][j];
+ }
+ }
+
+ // Tridiagonalize.
+ tred2();
+
+ // Diagonalize.
+ tql2();
+
+ } else {
+ H = new double[n][n];
+ ort = new double[n];
+
+ for (int j = 0; j < n; j++) {
+ for (int i = 0; i < n; i++) {
+ H[i][j] = A[i][j];
+ }
+ }
+
+ // Reduce to Hessenberg form.
+ orthes();
+
+ // Reduce Hessenberg to real Schur form.
+ hqr2();
+ }
+ }
+
+/* ------------------------
+ Public Methods
+ * ------------------------ */
+
+ /** Return the eigenvector matrix
+ @return V
+ */
+
+ public Matrix getV () {
+ return new Matrix(V,n,n);
+ }
+
+ /** Return the real parts of the eigenvalues
+ @return real(diag(D))
+ */
+
+ public double[] getRealEigenvalues () {
+ return d;
+ }
+
+ /** Return the imaginary parts of the eigenvalues
+ @return imag(diag(D))
+ */
+
+ public double[] getImagEigenvalues () {
+ return e;
+ }
+
+ /** Return the block diagonal eigenvalue matrix
+ @return D
+ */
+
+ public Matrix getD () {
+ Matrix X = new Matrix(n,n);
+ double[][] D = X.getArray();
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ D[i][j] = 0.0;
+ }
+ D[i][i] = d[i];
+ if (e[i] > 0) {
+ D[i][i+1] = e[i];
+ } else if (e[i] < 0) {
+ D[i][i-1] = e[i];
+ }
+ }
+ return X;
+ }
+}
diff --git a/src/wsi/ra/math/Jama/LUDecomposition.java b/src/wsi/ra/math/Jama/LUDecomposition.java
new file mode 100644
index 00000000..0a013056
--- /dev/null
+++ b/src/wsi/ra/math/Jama/LUDecomposition.java
@@ -0,0 +1,318 @@
+package wsi.ra.math.Jama;
+
+
+ /** LU Decomposition.
+
+ For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
+ unit lower triangular matrix L, an n-by-n upper triangular matrix U,
+ and a permutation vector piv of length m so that A(piv,:) = L*U.
+ If m < n, then L is m-by-m and U is m-by-n.
+
+ The LU decompostion with pivoting always exists, even if the matrix is
+ singular, so the constructor will never fail. The primary use of the
+ LU decomposition is in the solution of square systems of simultaneous
+ linear equations. This will fail if isNonsingular() returns false.
+ */
+
+public class LUDecomposition implements java.io.Serializable {
+
+/* ------------------------
+ Class variables
+ * ------------------------ */
+
+ /** Array for internal storage of decomposition.
+ @serial internal array storage.
+ */
+ private double[][] LU;
+
+ /** Row and column dimensions, and pivot sign.
+ @serial column dimension.
+ @serial row dimension.
+ @serial pivot sign.
+ */
+ private int m, n, pivsign;
+
+ /** Internal storage of pivot vector.
+ @serial pivot vector.
+ */
+ private int[] piv;
+
+/* ------------------------
+ Constructor
+ * ------------------------ */
+
+ /** LU Decomposition
+ @param A Rectangular matrix
+ @return Structure to access L, U and piv.
+ */
+
+ public LUDecomposition (Matrix A) {
+
+ // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+// System.out.println("A=");
+// A.print(A.getRowDimension(),A.getColumnDimension());
+
+ LU = A.getArrayCopy();
+ m = A.getRowDimension();
+ n = A.getColumnDimension();
+ piv = new int[m];
+ for (int i = 0; i < m; i++) {
+ piv[i] = i;
+ }
+ pivsign = 1;
+ double[] LUrowi;
+ double[] LUcolj = new double[m];
+
+ // Outer loop.
+
+ for (int j = 0; j < n; j++) {
+
+ // Make a copy of the j-th column to localize references.
+
+ 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];
+
+ // Most of the time is spent in the following dot product.
+
+ int kmax = Math.min(i,j);
+ double s = 0.0;
+ for (int k = 0; k < kmax; k++) {
+ s += LUrowi[k]*LUcolj[k];
+ }
+
+ LUrowi[j] = LUcolj[i] -= s;
+ }
+
+ // Find pivot and exchange if necessary.
+
+ 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];
+ }
+ }
+ }
+ }
+
+/* ------------------------
+ Temporary, experimental code.
+ ------------------------ *\
+
+ \** LU Decomposition, computed by Gaussian elimination.
+
+ This constructor computes L and U with the "daxpy"-based elimination
+ algorithm used in LINPACK and MATLAB. In Java, we suspect the dot-product,
+ Crout algorithm will be faster. We have temporarily included this
+ constructor until timing experiments confirm this suspicion.
+
+ @param A Rectangular matrix
+ @param linpackflag Use Gaussian elimination. Actual value ignored.
+ @return Structure to access L, U and piv.
+ *\
+
+ public LUDecomposition (Matrix A, int linpackflag) {
+ // Initialize.
+ LU = A.getArrayCopy();
+ m = A.getRowDimension();
+ n = A.getColumnDimension();
+ piv = new int[m];
+ for (int i = 0; i < m; i++) {
+ piv[i] = i;
+ }
+ pivsign = 1;
+ // Main loop.
+ for (int k = 0; k < n; k++) {
+ // Find pivot.
+ int p = k;
+ for (int i = k+1; i < m; i++) {
+ if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) {
+ p = i;
+ }
+ }
+ // Exchange if necessary.
+ if (p != k) {
+ for (int j = 0; j < n; j++) {
+ double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t;
+ }
+ int t = piv[p]; piv[p] = piv[k]; piv[k] = t;
+ pivsign = -pivsign;
+ }
+ // Compute multipliers and eliminate k-th column.
+ if (LU[k][k] != 0.0) {
+ for (int i = k+1; i < m; i++) {
+ LU[i][k] /= LU[k][k];
+ for (int j = k+1; j < n; j++) {
+ LU[i][j] -= LU[i][k]*LU[k][j];
+ }
+ }
+ }
+ }
+ }
+
+\* ------------------------
+ End of temporary code.
+ * ------------------------ */
+
+/* ------------------------
+ Public Methods
+ * ------------------------ */
+
+ /** Is the matrix nonsingular?
+ @return true if U, and hence A, is nonsingular.
+ */
+
+ public boolean isNonsingular () {
+ for (int j = 0; j < n; j++) {
+ //System.out.println("LU[j][j]"+LU[j][j]);
+ if (LU[j][j] == 0)
+ return false;
+
+ }
+ return true;
+ }
+
+ /** Return lower triangular factor
+ @return L
+ */
+
+ public Matrix getL () {
+ Matrix X = new Matrix(m,n);
+ double[][] L = X.getArray();
+ for (int i = 0; i < m; i++) {
+ for (int j = 0; j < n; j++) {
+ if (i > j) {
+ L[i][j] = LU[i][j];
+ } else if (i == j) {
+ L[i][j] = 1.0;
+ } else {
+ L[i][j] = 0.0;
+ }
+ }
+ }
+ return X;
+ }
+
+ /** Return upper triangular factor
+ @return U
+ */
+
+ public Matrix getU () {
+ Matrix X = new Matrix(n,n);
+ double[][] U = X.getArray();
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ if (i <= j) {
+ U[i][j] = LU[i][j];
+ } else {
+ U[i][j] = 0.0;
+ }
+ }
+ }
+ return X;
+ }
+
+ /** Return pivot permutation vector
+ @return piv
+ */
+
+ public int[] getPivot () {
+ int[] p = new int[m];
+ for (int i = 0; i < m; i++) {
+ p[i] = piv[i];
+ }
+ return p;
+ }
+
+ /** Return pivot permutation vector as a one-dimensional double array
+ @return (double) piv
+ */
+
+ public double[] getDoublePivot () {
+ double[] vals = new double[m];
+ for (int i = 0; i < m; i++) {
+ vals[i] = (double) piv[i];
+ }
+ return vals;
+ }
+
+ /** Determinant
+ @return det(A)
+ @exception IllegalArgumentException Matrix must be square
+ */
+
+ public double det () {
+ if (m != n) {
+ throw new IllegalArgumentException("Matrix must be square.");
+ }
+ double d = (double) pivsign;
+ for (int j = 0; j < n; j++) {
+ d *= LU[j][j];
+ }
+ return d;
+ }
+
+ /** Solve A*X = B
+ @param B A Matrix with as many rows as A and any number of columns.
+ @return X so that L*U*X = B(piv,:)
+ @exception IllegalArgumentException Matrix row dimensions must agree.
+ @exception RuntimeException Matrix is singular.
+ */
+
+ public Matrix solve (Matrix B) {
+ if (B.getRowDimension() != m) {
+ throw new IllegalArgumentException("Matrix row dimensions must agree.");
+ }
+ if (!this.isNonsingular()) {
+ //System.out.println("B=");
+ //B.print(B.getRowDimension(),B.getColumnDimension());
+ throw new RuntimeException("Matrix is singular.");
+ }
+
+ // Copy right hand side with pivoting
+ int nx = B.getColumnDimension();
+ Matrix Xmat = B.getMatrix(piv,0,nx-1);
+ double[][] X = Xmat.getArray();
+
+ // Solve L*Y = B(piv,:)
+ for (int k = 0; k < n; k++) {
+ for (int i = k+1; i < n; i++) {
+ for (int j = 0; j < nx; j++) {
+ X[i][j] -= X[k][j]*LU[i][k];
+ }
+ }
+ }
+ // Solve U*X = Y;
+ for (int k = n-1; k >= 0; k--) {
+ for (int j = 0; j < nx; j++) {
+ X[k][j] /= LU[k][k];
+ }
+ for (int i = 0; i < k; i++) {
+ for (int j = 0; j < nx; j++) {
+ X[i][j] -= X[k][j]*LU[i][k];
+ }
+ }
+ }
+ return Xmat;
+ }
+}
diff --git a/src/wsi/ra/math/Jama/Matrix.java b/src/wsi/ra/math/Jama/Matrix.java
new file mode 100644
index 00000000..455e2dd6
--- /dev/null
+++ b/src/wsi/ra/math/Jama/Matrix.java
@@ -0,0 +1,1111 @@
+package wsi.ra.math.Jama;
+
+
+import java.io.BufferedReader;
+import java.io.PrintWriter;
+import java.io.StreamTokenizer;
+import java.text.DecimalFormat;
+import java.text.DecimalFormatSymbols;
+import java.text.NumberFormat;
+import java.util.Locale;
+
+import wsi.ra.math.Jama.util.Maths;
+
+
+/**
+ Jama = Java Matrix class.
+
+ The Java Matrix Class provides the fundamental operations of numerical
+ linear algebra. Various constructors create Matrices from two dimensional
+ arrays of double precision floating point numbers. Various "gets" and
+ "sets" provide access to submatrices and matrix elements. Several methods
+ implement basic matrix arithmetic, including matrix addition and
+ multiplication, matrix norms, and element-by-element array operations.
+ Methods for reading and printing matrices are also included. All the
+ operations in this version of the Matrix Class involve real matrices.
+ Complex matrices may be handled in a future version.
+
+ Five fundamental matrix decompositions, which consist of pairs or triples
+ of matrices, permutation vectors, and the like, produce results in five
+ decomposition classes. These decompositions are accessed by the Matrix
+ class to compute solutions of simultaneous linear equations, determinants,
+ inverses and other matrix functions. The five decompositions are:
+
= n) throw new java.io.IOException
+ ("Row " + v.size() + " is too long.");
+ row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
+ } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
+ if (j < n) throw new java.io.IOException
+ ("Row " + v.size() + " is too short.");
+ }
+ int m = v.size(); // Now we've got the number of rows.
+ double[][] A = new double[m][];
+ v.copyInto(A); // copy the rows out of the vector
+ return new Matrix(A);
+ }
+
+
+/* ------------------------
+ Private Methods
+ * ------------------------ */
+
+ /** Check if size(A) == size(B) **/
+
+ private void checkMatrixDimensions (Matrix B) {
+ if (B.m != m || B.n != n) {
+ System.out.println("B.m"+B.m);
+ System.out.println("m"+m);
+ System.out.println("B.n"+B.n);
+ System.out.println("n"+n);
+ throw new IllegalArgumentException("Matrix dimensions must agree.");
+ }
+ }
+
+}
diff --git a/src/wsi/ra/math/Jama/QRDecomposition.java b/src/wsi/ra/math/Jama/QRDecomposition.java
new file mode 100644
index 00000000..76126eba
--- /dev/null
+++ b/src/wsi/ra/math/Jama/QRDecomposition.java
@@ -0,0 +1,218 @@
+package wsi.ra.math.Jama;
+import wsi.ra.math.Jama.util.Maths;
+
+/** QR Decomposition.
+
+ For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
+ orthogonal matrix Q and an n-by-n upper triangular matrix R so that
+ A = Q*R.
+
+ The QR decompostion always exists, even if the matrix does not have
+ full rank, so the constructor will never fail. The primary use of the
+ QR decomposition is in the least squares solution of nonsquare systems
+ of simultaneous linear equations. This will fail if isFullRank()
+ returns false.
+*/
+
+public class QRDecomposition implements java.io.Serializable {
+
+/* ------------------------
+ Class variables
+ * ------------------------ */
+
+ /** Array for internal storage of decomposition.
+ @serial internal array storage.
+ */
+ private double[][] QR;
+
+ /** Row and column dimensions.
+ @serial column dimension.
+ @serial row dimension.
+ */
+ private int m, n;
+
+ /** Array for internal storage of diagonal of R.
+ @serial diagonal of R.
+ */
+ private double[] Rdiag;
+
+/* ------------------------
+ Constructor
+ * ------------------------ */
+
+ /** QR Decomposition, computed by Householder reflections.
+ @param A Rectangular matrix
+ @return Structure to access R and the Householder vectors and compute Q.
+ */
+
+ public QRDecomposition (Matrix A) {
+ // Initialize.
+ QR = A.getArrayCopy();
+ m = A.getRowDimension();
+ n = A.getColumnDimension();
+ Rdiag = new double[n];
+
+ // Main loop.
+ for (int k = 0; k < n; k++) {
+ // Compute 2-norm of k-th column without under/overflow.
+ double nrm = 0;
+ for (int i = k; i < m; i++) {
+ nrm = Maths.hypot(nrm,QR[i][k]);
+ }
+
+ if (nrm != 0.0) {
+ // Form k-th Householder vector.
+ if (QR[k][k] < 0) {
+ nrm = -nrm;
+ }
+ for (int i = k; i < m; i++) {
+ QR[i][k] /= nrm;
+ }
+ QR[k][k] += 1.0;
+
+ // Apply transformation to remaining columns.
+ for (int j = k+1; j < n; j++) {
+ double s = 0.0;
+ for (int i = k; i < m; i++) {
+ s += QR[i][k]*QR[i][j];
+ }
+ s = -s/QR[k][k];
+ for (int i = k; i < m; i++) {
+ QR[i][j] += s*QR[i][k];
+ }
+ }
+ }
+ Rdiag[k] = -nrm;
+ }
+ }
+
+/* ------------------------
+ Public Methods
+ * ------------------------ */
+
+ /** Is the matrix full rank?
+ @return true if R, and hence A, has full rank.
+ */
+
+ public boolean isFullRank () {
+ for (int j = 0; j < n; j++) {
+ if (Rdiag[j] == 0)
+ return false;
+ }
+ return true;
+ }
+
+ /** Return the Householder vectors
+ @return Lower trapezoidal matrix whose columns define the reflections
+ */
+
+ public Matrix getH () {
+ Matrix X = new Matrix(m,n);
+ double[][] H = X.getArray();
+ for (int i = 0; i < m; i++) {
+ for (int j = 0; j < n; j++) {
+ if (i >= j) {
+ H[i][j] = QR[i][j];
+ } else {
+ H[i][j] = 0.0;
+ }
+ }
+ }
+ return X;
+ }
+
+ /** Return the upper triangular factor
+ @return R
+ */
+
+ public Matrix getR () {
+ Matrix X = new Matrix(n,n);
+ double[][] R = X.getArray();
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ if (i < j) {
+ R[i][j] = QR[i][j];
+ } else if (i == j) {
+ R[i][j] = Rdiag[i];
+ } else {
+ R[i][j] = 0.0;
+ }
+ }
+ }
+ return X;
+ }
+
+ /** Generate and return the (economy-sized) orthogonal factor
+ @return Q
+ */
+
+ public Matrix getQ () {
+ Matrix X = new Matrix(m,n);
+ double[][] Q = X.getArray();
+ for (int k = n-1; k >= 0; k--) {
+ for (int i = 0; i < m; i++) {
+ Q[i][k] = 0.0;
+ }
+ Q[k][k] = 1.0;
+ for (int j = k; j < n; j++) {
+ if (QR[k][k] != 0) {
+ double s = 0.0;
+ for (int i = k; i < m; i++) {
+ s += QR[i][k]*Q[i][j];
+ }
+ s = -s/QR[k][k];
+ for (int i = k; i < m; i++) {
+ Q[i][j] += s*QR[i][k];
+ }
+ }
+ }
+ }
+ return X;
+ }
+
+ /** Least squares solution of A*X = B
+ @param B A Matrix with as many rows as A and any number of columns.
+ @return X that minimizes the two norm of Q*R*X-B.
+ @exception IllegalArgumentException Matrix row dimensions must agree.
+ @exception RuntimeException Matrix is rank deficient.
+ */
+
+ public Matrix solve (Matrix B) {
+ if (B.getRowDimension() != m) {
+ throw new IllegalArgumentException("Matrix row dimensions must agree.");
+ }
+ if (!this.isFullRank()) {
+ throw new RuntimeException("Matrix is rank deficient.");
+ }
+
+ // Copy right hand side
+ int nx = B.getColumnDimension();
+ double[][] X = B.getArrayCopy();
+
+ // Compute Y = transpose(Q)*B
+ for (int k = 0; k < n; k++) {
+ for (int j = 0; j < nx; j++) {
+ double s = 0.0;
+ for (int i = k; i < m; i++) {
+ s += QR[i][k]*X[i][j];
+ }
+ s = -s/QR[k][k];
+ for (int i = k; i < m; i++) {
+ X[i][j] += s*QR[i][k];
+ }
+ }
+ }
+ // Solve R*X = Y;
+ for (int k = n-1; k >= 0; k--) {
+ for (int j = 0; j < nx; j++) {
+ X[k][j] /= Rdiag[k];
+ }
+ for (int i = 0; i < k; i++) {
+ for (int j = 0; j < nx; j++) {
+ X[i][j] -= X[k][j]*QR[i][k];
+ }
+ }
+ }
+ return (new Matrix(X,n,nx).getMatrix(0,n-1,0,nx-1));
+ }
+}
diff --git a/src/wsi/ra/math/Jama/SingularValueDecomposition.java b/src/wsi/ra/math/Jama/SingularValueDecomposition.java
new file mode 100644
index 00000000..949c65e2
--- /dev/null
+++ b/src/wsi/ra/math/Jama/SingularValueDecomposition.java
@@ -0,0 +1,540 @@
+package wsi.ra.math.Jama;
+import wsi.ra.math.Jama.util.*;
+
+
+ /** Singular Value Decomposition.
+
+ For an m-by-n matrix A with m >= n, the singular value decomposition is
+ an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
+ an n-by-n orthogonal matrix V so that A = U*S*V'.
+
+ The singular values, sigma[k] = S[k][k], are ordered so that
+ sigma[0] >= sigma[1] >= ... >= sigma[n-1].
+
+ The singular value decompostion always exists, so the constructor will
+ never fail. The matrix condition number and the effective numerical
+ rank can be computed from this decomposition.
+ */
+
+public class SingularValueDecomposition implements java.io.Serializable {
+
+/* ------------------------
+ Class variables
+ * ------------------------ */
+
+ /** Arrays for internal storage of U and V.
+ @serial internal storage of U.
+ @serial internal storage of V.
+ */
+ private double[][] U, V;
+
+ /** Array for internal storage of singular values.
+ @serial internal storage of singular values.
+ */
+ private double[] s;
+
+ /** Row and column dimensions.
+ @serial row dimension.
+ @serial column dimension.
+ */
+ private int m, n;
+
+/* ------------------------
+ Constructor
+ * ------------------------ */
+
+ /** Construct the singular value decomposition
+ @param A Rectangular matrix
+ @return Structure to access U, S and V.
+ */
+
+ public SingularValueDecomposition (Matrix Arg) {
+
+ // Derived from LINPACK code.
+ // Initialize.
+ double[][] A = Arg.getArrayCopy();
+ m = Arg.getRowDimension();
+ n = Arg.getColumnDimension();
+ int nu = Math.min(m,n);
+ s = new double [Math.min(m+1,n)];
+ U = new double [m][nu];
+ V = new double [n][n];
+ double[] e = new double [n];
+ double[] work = new double [m];
+ boolean wantu = true;
+ boolean wantv = true;
+
+ // Reduce A to bidiagonal form, storing the diagonal elements
+ // in s and the super-diagonal elements in e.
+
+ int nct = Math.min(m-1,n);
+ int nrt = Math.max(0,Math.min(n-2,m));
+ for (int k = 0; k < Math.max(nct,nrt); k++) {
+ if (k < nct) {
+
+ // Compute the transformation for the k-th column and
+ // place the k-th diagonal in s[k].
+ // Compute 2-norm of k-th column without under/overflow.
+ s[k] = 0;
+ for (int i = k; i < m; i++) {
+ s[k] = Maths.hypot(s[k],A[i][k]);
+ }
+ if (s[k] != 0.0) {
+ if (A[k][k] < 0.0) {
+ s[k] = -s[k];
+ }
+ for (int i = k; i < m; i++) {
+ A[i][k] /= s[k];
+ }
+ A[k][k] += 1.0;
+ }
+ s[k] = -s[k];
+ }
+ for (int j = k+1; j < n; j++) {
+ if ((k < nct) & (s[k] != 0.0)) {
+
+ // Apply the transformation.
+
+ double t = 0;
+ for (int i = k; i < m; i++) {
+ t += A[i][k]*A[i][j];
+ }
+ t = -t/A[k][k];
+ for (int i = k; i < m; i++) {
+ A[i][j] += t*A[i][k];
+ }
+ }
+
+ // Place the k-th row of A into e for the
+ // subsequent calculation of the row transformation.
+
+ e[j] = A[k][j];
+ }
+ if (wantu & (k < nct)) {
+
+ // Place the transformation in U for subsequent back
+ // multiplication.
+
+ for (int i = k; i < m; i++) {
+ U[i][k] = A[i][k];
+ }
+ }
+ if (k < nrt) {
+
+ // Compute the k-th row transformation and place the
+ // k-th super-diagonal in e[k].
+ // Compute 2-norm without under/overflow.
+ e[k] = 0;
+ for (int i = k+1; i < n; i++) {
+ e[k] = Maths.hypot(e[k],e[i]);
+ }
+ if (e[k] != 0.0) {
+ if (e[k+1] < 0.0) {
+ e[k] = -e[k];
+ }
+ for (int i = k+1; i < n; i++) {
+ e[i] /= e[k];
+ }
+ e[k+1] += 1.0;
+ }
+ e[k] = -e[k];
+ if ((k+1 < m) & (e[k] != 0.0)) {
+
+ // Apply the transformation.
+
+ for (int i = k+1; i < m; i++) {
+ work[i] = 0.0;
+ }
+ for (int j = k+1; j < n; j++) {
+ for (int i = k+1; i < m; i++) {
+ work[i] += e[j]*A[i][j];
+ }
+ }
+ for (int j = k+1; j < n; j++) {
+ double t = -e[j]/e[k+1];
+ for (int i = k+1; i < m; i++) {
+ A[i][j] += t*work[i];
+ }
+ }
+ }
+ if (wantv) {
+
+ // Place the transformation in V for subsequent
+ // back multiplication.
+
+ for (int i = k+1; i < n; i++) {
+ V[i][k] = e[i];
+ }
+ }
+ }
+ }
+
+ // Set up the final bidiagonal matrix or order p.
+
+ int p = Math.min(n,m+1);
+ if (nct < n) {
+ s[nct] = A[nct][nct];
+ }
+ if (m < p) {
+ s[p-1] = 0.0;
+ }
+ if (nrt+1 < p) {
+ e[nrt] = A[nrt][p-1];
+ }
+ e[p-1] = 0.0;
+
+ // If required, generate U.
+
+ if (wantu) {
+ for (int j = nct; j < nu; j++) {
+ for (int i = 0; i < m; i++) {
+ U[i][j] = 0.0;
+ }
+ U[j][j] = 1.0;
+ }
+ for (int k = nct-1; k >= 0; k--) {
+ if (s[k] != 0.0) {
+ for (int j = k+1; j < nu; j++) {
+ double t = 0;
+ for (int i = k; i < m; i++) {
+ t += U[i][k]*U[i][j];
+ }
+ t = -t/U[k][k];
+ for (int i = k; i < m; i++) {
+ U[i][j] += t*U[i][k];
+ }
+ }
+ for (int i = k; i < m; i++ ) {
+ U[i][k] = -U[i][k];
+ }
+ U[k][k] = 1.0 + U[k][k];
+ for (int i = 0; i < k-1; i++) {
+ U[i][k] = 0.0;
+ }
+ } else {
+ for (int i = 0; i < m; i++) {
+ U[i][k] = 0.0;
+ }
+ U[k][k] = 1.0;
+ }
+ }
+ }
+
+ // If required, generate V.
+
+ if (wantv) {
+ for (int k = n-1; k >= 0; k--) {
+ if ((k < nrt) & (e[k] != 0.0)) {
+ for (int j = k+1; j < nu; j++) {
+ double t = 0;
+ for (int i = k+1; i < n; i++) {
+ t += V[i][k]*V[i][j];
+ }
+ t = -t/V[k+1][k];
+ for (int i = k+1; i < n; i++) {
+ V[i][j] += t*V[i][k];
+ }
+ }
+ }
+ for (int i = 0; i < n; i++) {
+ V[i][k] = 0.0;
+ }
+ V[k][k] = 1.0;
+ }
+ }
+
+ // Main iteration loop for the singular values.
+
+ int pp = p-1;
+ int iter = 0;
+ double eps = Math.pow(2.0,-52.0);
+ while (p > 0) {
+ int k,kase;
+
+ // Here is where a test for too many iterations would go.
+
+ // This section of the program inspects for
+ // negligible elements in the s and e arrays. On
+ // completion the variables kase and k are set as follows.
+
+ // kase = 1 if s(p) and e[k-1] are negligible and k
= -1; k--) {
+ if (k == -1) {
+ break;
+ }
+ if (Math.abs(e[k]) <= eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) {
+ e[k] = 0.0;
+ break;
+ }
+ }
+ if (k == p-2) {
+ kase = 4;
+ } else {
+ int ks;
+ for (ks = p-1; ks >= k; ks--) {
+ if (ks == k) {
+ break;
+ }
+ double t = (ks != p ? Math.abs(e[ks]) : 0.) +
+ (ks != k+1 ? Math.abs(e[ks-1]) : 0.);
+ if (Math.abs(s[ks]) <= eps*t) {
+ s[ks] = 0.0;
+ break;
+ }
+ }
+ if (ks == k) {
+ kase = 3;
+ } else if (ks == p-1) {
+ kase = 1;
+ } else {
+ kase = 2;
+ k = ks;
+ }
+ }
+ k++;
+
+ // Perform the task indicated by kase.
+
+ switch (kase) {
+
+ // Deflate negligible s(p).
+
+ case 1: {
+ double f = e[p-2];
+ e[p-2] = 0.0;
+ for (int j = p-2; j >= k; j--) {
+ double t = Maths.hypot(s[j],f);
+ double cs = s[j]/t;
+ double sn = f/t;
+ s[j] = t;
+ if (j != k) {
+ f = -sn*e[j-1];
+ e[j-1] = cs*e[j-1];
+ }
+ if (wantv) {
+ for (int i = 0; i < n; i++) {
+ t = cs*V[i][j] + sn*V[i][p-1];
+ V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
+ V[i][j] = t;
+ }
+ }
+ }
+ }
+ break;
+
+ // Split at negligible s(k).
+
+ case 2: {
+ double f = e[k-1];
+ e[k-1] = 0.0;
+ for (int j = k; j < p; j++) {
+ double t = Maths.hypot(s[j],f);
+ double cs = s[j]/t;
+ double sn = f/t;
+ s[j] = t;
+ f = -sn*e[j];
+ e[j] = cs*e[j];
+ if (wantu) {
+ for (int i = 0; i < m; i++) {
+ t = cs*U[i][j] + sn*U[i][k-1];
+ U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
+ U[i][j] = t;
+ }
+ }
+ }
+ }
+ break;
+
+ // Perform one qr step.
+
+ case 3: {
+
+ // Calculate the shift.
+
+ double scale = Math.max(Math.max(Math.max(Math.max(
+ Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])),
+ Math.abs(s[k])),Math.abs(e[k]));
+ double sp = s[p-1]/scale;
+ double spm1 = s[p-2]/scale;
+ double epm1 = e[p-2]/scale;
+ double sk = s[k]/scale;
+ double ek = e[k]/scale;
+ double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
+ double c = (sp*epm1)*(sp*epm1);
+ double shift = 0.0;
+ if ((b != 0.0) | (c != 0.0)) {
+ shift = Math.sqrt(b*b + c);
+ if (b < 0.0) {
+ shift = -shift;
+ }
+ shift = c/(b + shift);
+ }
+ double f = (sk + sp)*(sk - sp) + shift;
+ double g = sk*ek;
+
+ // Chase zeros.
+
+ for (int j = k; j < p-1; j++) {
+ double t = Maths.hypot(f,g);
+ double cs = f/t;
+ double sn = g/t;
+ if (j != k) {
+ e[j-1] = t;
+ }
+ f = cs*s[j] + sn*e[j];
+ e[j] = cs*e[j] - sn*s[j];
+ g = sn*s[j+1];
+ s[j+1] = cs*s[j+1];
+ if (wantv) {
+ for (int i = 0; i < n; i++) {
+ t = cs*V[i][j] + sn*V[i][j+1];
+ V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
+ V[i][j] = t;
+ }
+ }
+ t = Maths.hypot(f,g);
+ cs = f/t;
+ sn = g/t;
+ s[j] = t;
+ f = cs*e[j] + sn*s[j+1];
+ s[j+1] = -sn*e[j] + cs*s[j+1];
+ g = sn*e[j+1];
+ e[j+1] = cs*e[j+1];
+ if (wantu && (j < m-1)) {
+ for (int i = 0; i < m; i++) {
+ t = cs*U[i][j] + sn*U[i][j+1];
+ U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
+ U[i][j] = t;
+ }
+ }
+ }
+ e[p-2] = f;
+ iter = iter + 1;
+ }
+ break;
+
+ // Convergence.
+
+ case 4: {
+
+ // Make the singular values positive.
+
+ if (s[k] <= 0.0) {
+ s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
+ if (wantv) {
+ for (int i = 0; i <= pp; i++) {
+ V[i][k] = -V[i][k];
+ }
+ }
+ }
+
+ // Order the singular values.
+
+ while (k < pp) {
+ if (s[k] >= s[k+1]) {
+ break;
+ }
+ double t = s[k];
+ s[k] = s[k+1];
+ s[k+1] = t;
+ if (wantv && (k < n-1)) {
+ for (int i = 0; i < n; i++) {
+ t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
+ }
+ }
+ if (wantu && (k < m-1)) {
+ for (int i = 0; i < m; i++) {
+ t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
+ }
+ }
+ k++;
+ }
+ iter = 0;
+ p--;
+ }
+ break;
+ }
+ }
+ }
+
+/* ------------------------
+ Public Methods
+ * ------------------------ */
+
+ /** Return the left singular vectors
+ @return U
+ */
+
+ public Matrix getU () {
+ return new Matrix(U,m,Math.min(m+1,n));
+ }
+
+ /** Return the right singular vectors
+ @return V
+ */
+
+ public Matrix getV () {
+ return new Matrix(V,n,n);
+ }
+
+ /** Return the one-dimensional array of singular values
+ @return diagonal of S.
+ */
+
+ public double[] getSingularValues () {
+ return s;
+ }
+
+ /** Return the diagonal matrix of singular values
+ @return S
+ */
+
+ public Matrix getS () {
+ Matrix X = new Matrix(n,n);
+ double[][] S = X.getArray();
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ S[i][j] = 0.0;
+ }
+ S[i][i] = this.s[i];
+ }
+ return X;
+ }
+
+ /** Two norm
+ @return max(S)
+ */
+
+ public double norm2 () {
+ return s[0];
+ }
+
+ /** Two norm condition number
+ @return max(S)/min(S)
+ */
+
+ public double cond () {
+ return s[0]/s[Math.min(m,n)-1];
+ }
+
+ /** Effective numerical matrix rank
+ @return Number of nonnegligible singular values.
+ */
+
+ public int rank () {
+ double eps = Math.pow(2.0,-52.0);
+ double tol = Math.max(m,n)*s[0]*eps;
+ int r = 0;
+ for (int i = 0; i < s.length; i++) {
+ if (s[i] > tol) {
+ r++;
+ }
+ }
+ return r;
+ }
+}
diff --git a/src/wsi/ra/math/Jama/util/Maths.java b/src/wsi/ra/math/Jama/util/Maths.java
new file mode 100644
index 00000000..bb6620df
--- /dev/null
+++ b/src/wsi/ra/math/Jama/util/Maths.java
@@ -0,0 +1,21 @@
+package wsi.ra.math.Jama.util;
+
+public class Maths {
+
+ /** sqrt(a^2 + b^2) without under/overflow. **/
+
+ public static double hypot(double a, double b) {
+ double r;
+ double aa = Math.abs(a);
+ double bb = Math.abs(b);
+ if (aa > bb) {
+ r = b/a;
+ r = aa*Math.sqrt(1+r*r);
+ } else if (b != 0) {
+ r = a/b;
+ r = bb*Math.sqrt(1+r*r);
+ } else
+ r = 0.0;
+ return r;
+ }
+}
diff --git a/src/wsi/ra/math/RNG.java b/src/wsi/ra/math/RNG.java
new file mode 100644
index 00000000..bba058f9
--- /dev/null
+++ b/src/wsi/ra/math/RNG.java
@@ -0,0 +1,184 @@
+package wsi.ra.math;
+/**
+ * Title: JavaEvA
+ * Description:
+ * Copyright: Copyright (c) 2003
+ * Company: University of Tuebingen, Computer Architecture
+ * @author Holger Ulmer, Felix Streichert, Hannes Planatscher
+ * @version: $Revision: 1.1.1.1 $
+ * $Date: 2003/07/03 14:59:40 $
+ * $Author: ulmerh $
+ */
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+import java.util.Random;
+/*==========================================================================*
+ * CLASS DECLARATION
+ *==========================================================================*/
+/**
+ *
+ */
+public class RNG extends Random {
+ private static Random random;
+ private static long randomSeed;
+ /**
+ *
+ */
+ static {
+ randomSeed=System.currentTimeMillis();
+ random=new Random(randomSeed);
+ }
+ /**
+ *
+ */
+ public static void setseed(long x) {
+ randomSeed=x;
+ if (x==0)
+ randomSeed=System.currentTimeMillis();
+ if (x==999)
+ return;
+ random=new Random(randomSeed);
+ }
+ /**
+ *
+ */
+ public static void setRandomseed() {
+ randomSeed=System.currentTimeMillis();
+ random=new Random(randomSeed);
+ }
+ /**
+ *
+ */
+ public static void setRandom(Random base_random) {
+ random=base_random;
+ }
+ /**
+ *
+ */
+ public static void setRandomSeed(long new_seed){
+ randomSeed=new_seed;
+ random.setSeed(randomSeed);
+ }
+ /**
+ *
+ */
+ public static long getRandomSeed() {
+ return randomSeed;
+ }
+ /**
+ *
+ */
+ public static int randomInt() {
+ return randomInt(0,1);
+ }
+ public static int randomInt(int lo,int hi) {
+ return (Math.abs(random.nextInt())%(hi-lo+1))+lo;
+ }
+ /**
+ *
+ */
+ public static long randomLong() {
+ return randomLong(0,1);
+ }
+ /**
+ *
+ */
+ public static long randomLong(long lo,long hi) {
+ return (Math.abs(random.nextLong())%(hi-lo+1))+lo;
+ }
+ /**
+ *
+ */
+ public static float randomFloat() {
+ return random.nextFloat();
+ }
+ /**
+ *
+ */
+ public static float randomFloat(float lo,float hi) {
+ return (hi-lo)*random.nextFloat()+lo;
+ }
+ /**
+ *
+ */
+ public static double randomDouble() {
+ return random.nextDouble();
+ }
+ /**
+ *
+ */
+ public static double randomDouble(double lo,double hi) {
+ return (hi-lo)*random.nextDouble()+lo;
+ }
+ /**
+ *
+ */
+ public static double[] randomDoubleArray(double[] lo,double[] hi) {
+ double[] xin = new double[lo.length];
+ for (int i=0;i
+ * NOTE: These special functions do not necessarily use the fastest
+ * or most accurate algorithms.
+ *
+ * @version $Revision: 1.1.1.1 $, $Date: 2003/07/03 14:59:40 $
+ * @author Leigh Brookshaw
+ */
+
+
+public final class SpecialFunction extends Object {
+
+ /*
+ ** machine constants
+ */
+ private static final double MACHEP = 1.11022302462515654042E-16;
+ private static final double MAXLOG = 7.09782712893383996732E2;
+ private static final double MINLOG = -7.451332191019412076235E2;
+ private static final double MAXGAM = 171.624376956302725;
+ private static final double SQTPI = 2.50662827463100050242E0;
+ private static final double SQRTH = 7.07106781186547524401E-1;
+ private static final double LOGPI = 1.14472988584940017414;
+
+
+ /*
+ ** Physical Constants in cgs Units
+ */
+
+
+ /**
+ * Boltzman Constant. Units erg/deg(K)
+ */
+ public static final double BOLTZMAN = 1.3807e-16;
+ /**
+ * Elementary Charge. Units statcoulomb
+ */
+ public static final double ECHARGE = 4.8032e-10;
+ /**
+ * Electron Mass. Units g
+ */
+ public static final double EMASS = 9.1095e-28;
+ /**
+ * Proton Mass. Units g
+ */
+ public static final double PMASS = 1.6726e-24;
+ /**
+ * Gravitational Constant. Units dyne-cm^2/g^2
+ */
+ public static final double GRAV = 6.6720e-08;
+ /**
+ * Planck constant. Units erg-sec
+ */
+ public static final double PLANCK = 6.6262e-27;
+ /**
+ * Speed of Light in a Vacuum. Units cm/sec
+ */
+ public static final double LIGHTSPEED = 2.9979e10;
+ /**
+ * Stefan-Boltzman Constant. Units erg/cm^2-sec-deg^4
+ */
+ public static final double STEFANBOLTZ = 5.6703e-5;
+ /**
+ * Avogadro Number. Units 1/mol
+ */
+ public static final double AVOGADRO = 6.0220e23;
+ /**
+ * Gas Constant. Units erg/deg-mol
+ */
+ public static final double GASCONSTANT = 8.3144e07;
+ /**
+ * Gravitational Acceleration at the Earths surface. Units cm/sec^2
+ */
+ public static final double GRAVACC = 980.67;
+
+ /**
+ * Solar Mass. Units g
+ */
+ public static final double SOLARMASS = 1.99e33;
+ /**
+ * Solar Radius. Units cm
+ */
+ public static final double SOLARRADIUS = 6.96e10;
+ /**
+ * Solar Luminosity. Units erg/sec
+ */
+ public static final double SOLARLUM = 3.90e33;
+ /**
+ * Solar Flux. Units erg/cm^2-sec
+ */
+ public static final double SOLARFLUX = 6.41e10;
+ /**
+ * Astronomical Unit (radius of the Earth's orbit). Units cm
+ */
+ public static final double AU = 1.50e13;
+
+
+ /**
+ * Don't let anyone instantiate this class.
+ */
+ private SpecialFunction() {}
+
+ /*
+ ** Function Methods
+ */
+
+ /**
+ * @param x a double value
+ * @return The log10
+ */
+ static public double log10(double x) throws ArithmeticException {
+ if( x <= 0.0 ) throw new ArithmeticException("range exception");
+ return Math.log(x)/2.30258509299404568401;
+ }
+
+
+ /**
+ * @param x a double value
+ * @return the hyperbolic cosine of the argument
+ */
+
+ static public double cosh(double x) throws ArithmeticException {
+ double a;
+ a = x;
+ if( a < 0.0 ) a = Math.abs(x);
+ a = Math.exp(a);
+ return 0.5*(a+1/a);
+ }
+
+ /**
+ * @param x a double value
+ * @return the hyperbolic sine of the argument
+ */
+ static public double sinh(double x) throws ArithmeticException {
+ double a;
+ if(x == 0.0) return x;
+ a = x;
+ if( a < 0.0 ) a = Math.abs(x);
+ a = Math.exp(a);
+ if( x < 0.0 ) return -0.5*(a-1/a);
+ else return 0.5*(a-1/a);
+ }
+
+ /**
+ * @param x a double value
+ * @return the hyperbolic tangent of the argument
+ */
+ static public double tanh(double x) throws ArithmeticException {
+ double a;
+ if( x == 0.0 ) return x;
+ a = x;
+ if( a < 0.0 ) a = Math.abs(x);
+ a = Math.exp(2.0*a);
+ if(x < 0.0 ) return -( 1.0-2.0/(a+1.0) );
+ else return ( 1.0-2.0/(a+1.0) );
+ }
+
+ /**
+ * @param x a double value
+ * @return the hyperbolic arc cosine of the argument
+ */
+
+ static public double acosh(double x) throws ArithmeticException {
+ if( x < 1.0 ) throw new ArithmeticException("range exception");
+ return Math.log( x + Math.sqrt(x*x-1));
+ }
+
+ /**
+ * @param x a double value
+ * @return the hyperbolic arc sine of the argument
+ */
+ static public double asinh(double xx) throws ArithmeticException {
+ double x;
+ int sign;
+ if(xx == 0.0) return xx;
+ if( xx < 0.0 ) {
+ sign = -1;
+ x = -xx;
+ } else {
+ sign = 1;
+ x = xx;
+ }
+ return sign*Math.log( x + Math.sqrt(x*x+1));
+ }
+
+ /**
+ * @param x a double value
+ * @return the hyperbolic arc tangent of the argument
+ */
+ static public double atanh(double x) throws ArithmeticException {
+ if( x > 1.0 || x < -1.0 ) throw
+ new ArithmeticException("range exception");
+ return 0.5 * Math.log( (1.0+x)/(1.0-x) );
+ }
+
+ /**
+ * @param x a double value
+ * @return the Bessel function of order 0 of the argument.
+ */
+
+ static public double j0(double x) throws ArithmeticException {
+ double ax;
+
+ if( (ax=Math.abs(x)) < 8.0 ) {
+ double y=x*x;
+ double ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
+ +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
+ double ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
+ +y*(59272.64853+y*(267.8532712+y*1.0))));
+
+ return ans1/ans2;
+
+ } else {
+ double z=8.0/ax;
+ double y=z*z;
+ double xx=ax-0.785398164;
+ double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
+ +y*(-0.2073370639e-5+y*0.2093887211e-6)));
+ double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
+ +y*(-0.6911147651e-5+y*(0.7621095161e-6
+ -y*0.934935152e-7)));
+
+ return Math.sqrt(0.636619772/ax)*
+ (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
+ }
+ }
+ /**
+ * @param x a double value
+ * @return the Bessel function of order 1 of the argument.
+ */
+
+ static public double j1(double x) throws ArithmeticException {
+
+ double ax;
+ double y;
+ double ans1, ans2;
+
+ if ((ax=Math.abs(x)) < 8.0) {
+ y=x*x;
+ ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+ +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
+ ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+ +y*(99447.43394+y*(376.9991397+y*1.0))));
+ return ans1/ans2;
+ } else {
+ double z=8.0/ax;
+ double xx=ax-2.356194491;
+ y=z*z;
+
+ ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+ +y*(0.2457520174e-5+y*(-0.240337019e-6))));
+ ans2=0.04687499995+y*(-0.2002690873e-3
+ +y*(0.8449199096e-5+y*(-0.88228987e-6
+ +y*0.105787412e-6)));
+ double ans=Math.sqrt(0.636619772/ax)*
+ (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
+ if (x < 0.0) ans = -ans;
+ return ans;
+ }
+ }
+
+ /**
+ * @param n integer order
+ * @param x a double value
+ * @return the Bessel function of order n of the argument.
+ */
+ static public double jn(int n, double x) throws ArithmeticException {
+ int j,m;
+ double ax,bj,bjm,bjp,sum,tox,ans;
+ boolean jsum;
+
+ double ACC = 40.0;
+ double BIGNO = 1.0e+10;
+ double BIGNI = 1.0e-10;
+
+ if(n == 0) return j0(x);
+ if(n == 1) return j1(x);
+
+ ax=Math.abs(x);
+ if(ax == 0.0) return 0.0;
+ else
+ if (ax > (double)n) {
+ tox=2.0/ax;
+ bjm=j0(ax);
+ bj=j1(ax);
+ for (j=1;j0;j--) {
+ bjm=j*tox*bj-bjp;
+ bjp=bj;
+ bj=bjm;
+ if (Math.abs(bj) > BIGNO) {
+ bj *= BIGNI;
+ bjp *= BIGNI;
+ ans *= BIGNI;
+ sum *= BIGNI;
+ }
+ if (jsum) sum += bj;
+ jsum=!jsum;
+ if (j == n) ans=bjp;
+ }
+ sum=2.0*sum-bj;
+ ans /= sum;
+ }
+ return x < 0.0 && n%2 == 1 ? -ans : ans;
+ }
+ /**
+ * @param x a double value
+ * @return the Bessel function of the second kind,
+ * of order 0 of the argument.
+ */
+
+ static public double y0(double x) throws ArithmeticException {
+
+ if (x < 8.0) {
+ double y=x*x;
+
+ double ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
+ +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
+ double ans2=40076544269.0+y*(745249964.8+y*(7189466.438
+ +y*(47447.26470+y*(226.1030244+y*1.0))));
+
+ return (ans1/ans2)+0.636619772*j0(x)*Math.log(x);
+ } else {
+ double z=8.0/x;
+ double y=z*z;
+ double xx=x-0.785398164;
+
+ double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
+ +y*(-0.2073370639e-5+y*0.2093887211e-6)));
+ double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
+ +y*(-0.6911147651e-5+y*(0.7621095161e-6
+ +y*(-0.934945152e-7))));
+ return Math.sqrt(0.636619772/x)*
+ (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
+ }
+ }
+
+ /**
+ * @param x a double value
+ * @return the Bessel function of the second kind,
+ * of order 1 of the argument.
+ */
+ static public double y1(double x) throws ArithmeticException {
+
+ if (x < 8.0) {
+ double y=x*x;
+ double ans1=x*(-0.4900604943e13+y*(0.1275274390e13
+ +y*(-0.5153438139e11+y*(0.7349264551e9
+ +y*(-0.4237922726e7+y*0.8511937935e4)))));
+ double ans2=0.2499580570e14+y*(0.4244419664e12
+ +y*(0.3733650367e10+y*(0.2245904002e8
+ +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
+ return (ans1/ans2)+0.636619772*(j1(x)*Math.log(x)-1.0/x);
+ } else {
+ double z=8.0/x;
+ double y=z*z;
+ double xx=x-2.356194491;
+ double ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+ +y*(0.2457520174e-5+y*(-0.240337019e-6))));
+ double ans2=0.04687499995+y*(-0.2002690873e-3
+ +y*(0.8449199096e-5+y*(-0.88228987e-6
+ +y*0.105787412e-6)));
+ return Math.sqrt(0.636619772/x)*
+ (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
+ }
+ }
+ /**
+ * @param n integer order
+ * @param x a double value
+ * @return the Bessel function of the second kind,
+ * of order n of the argument.
+ */
+ static public double yn(int n, double x) throws ArithmeticException {
+ double by,bym,byp,tox;
+
+ if(n == 0) return y0(x);
+ if(n == 1) return y1(x);
+
+ tox=2.0/x;
+ by=y1(x);
+ bym=y0(x);
+ for (int j=1;j 1) { d *= i--; }
+ if(j < 0) return -d;
+ else return d;
+ }
+
+
+
+
+
+ /**
+ * @param x a double value
+ * @return the Gamma function of the value.
+ *
+ *
+ * Converted to Java from
+ * Cephes Math Library Release 2.2: July, 1992
+ * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ **/
+ static public double gamma(double x) throws ArithmeticException {
+
+ double P[] = {
+ 1.60119522476751861407E-4,
+ 1.19135147006586384913E-3,
+ 1.04213797561761569935E-2,
+ 4.76367800457137231464E-2,
+ 2.07448227648435975150E-1,
+ 4.94214826801497100753E-1,
+ 9.99999999999999996796E-1
+ };
+ double Q[] = {
+ -2.31581873324120129819E-5,
+ 5.39605580493303397842E-4,
+ -4.45641913851797240494E-3,
+ 1.18139785222060435552E-2,
+ 3.58236398605498653373E-2,
+ -2.34591795718243348568E-1,
+ 7.14304917030273074085E-2,
+ 1.00000000000000000320E0
+ };
+ double MAXGAM = 171.624376956302725;
+ double LOGPI = 1.14472988584940017414;
+
+ double p, z;
+ int i;
+
+ double q = Math.abs(x);
+
+ if( q > 33.0 ) {
+ if( x < 0.0 ) {
+ p = Math.floor(q);
+ if( p == q ) throw new ArithmeticException("gamma: overflow");
+ i = (int)p;
+ z = q - p;
+ if( z > 0.5 ) {
+ p += 1.0;
+ z = q - p;
+ }
+ z = q * Math.sin( Math.PI * z );
+ if( z == 0.0 ) throw new ArithmeticException("gamma: overflow");
+ z = Math.abs(z);
+ z = Math.PI/(z * stirf(q) );
+
+ return -z;
+ } else {
+ return stirf(x);
+ }
+ }
+
+ z = 1.0;
+ while( x >= 3.0 ) {
+ x -= 1.0;
+ z *= x;
+ }
+
+ while( x < 0.0 ) {
+ if( x == 0.0 ) {
+ throw new ArithmeticException("gamma: singular");
+ } else
+ if( x > -1.E-9 ) {
+ return( z/((1.0 + 0.5772156649015329 * x) * x) );
+ }
+ z /= x;
+ x += 1.0;
+ }
+
+ while( x < 2.0 ) {
+ if( x == 0.0 ) {
+ throw new ArithmeticException("gamma: singular");
+ } else
+ if( x < 1.e-9 ) {
+ return( z/((1.0 + 0.5772156649015329 * x) * x) );
+ }
+ z /= x;
+ x += 1.0;
+ }
+
+ if( (x == 2.0) || (x == 3.0) ) return z;
+
+ x -= 2.0;
+ p = polevl( x, P, 6 );
+ q = polevl( x, Q, 7 );
+ return z * p / q;
+
+ }
+
+/* Gamma function computed by Stirling's formula.
+ * The polynomial STIR is valid for 33 <= x <= 172.
+
+Cephes Math Library Release 2.2: July, 1992
+Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+ static private double stirf(double x) throws ArithmeticException {
+ double STIR[] = {
+ 7.87311395793093628397E-4,
+ -2.29549961613378126380E-4,
+ -2.68132617805781232825E-3,
+ 3.47222221605458667310E-3,
+ 8.33333333333482257126E-2,
+ };
+ double MAXSTIR = 143.01608;
+
+ double w = 1.0/x;
+ double y = Math.exp(x);
+
+ w = 1.0 + w * polevl( w, STIR, 4 );
+
+ if( x > MAXSTIR ) {
+ /* Avoid overflow in Math.pow() */
+ double v = Math.pow( x, 0.5 * x - 0.25 );
+ y = v * (v / y);
+ } else {
+ y = Math.pow( x, x - 0.5 ) / y;
+ }
+ y = SQTPI * y * w;
+ return y;
+ }
+
+ /**
+ * @param a double value
+ * @param x double value
+ * @return the Complemented Incomplete Gamma function.
+ *
+ *
+ * Converted to Java from
+ * Cephes Math Library Release 2.2: July, 1992
+ * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ **/
+
+ static public double igamc( double a, double x )
+ throws ArithmeticException {
+ double big = 4.503599627370496e15;
+ double biginv = 2.22044604925031308085e-16;
+ double ans, ax, c, yc, r, t, y, z;
+ double pk, pkm1, pkm2, qk, qkm1, qkm2;
+
+ if( x <= 0 || a <= 0 ) return 1.0;
+
+ if( x < 1.0 || x < a ) return 1.0 - igam(a,x);
+
+ ax = a * Math.log(x) - x - lgamma(a);
+ if( ax < -MAXLOG ) return 0.0;
+
+ ax = Math.exp(ax);
+
+ /* continued fraction */
+ y = 1.0 - a;
+ z = x + y + 1.0;
+ c = 0.0;
+ pkm2 = 1.0;
+ qkm2 = x;
+ pkm1 = x + 1.0;
+ qkm1 = z * x;
+ ans = pkm1/qkm1;
+
+ do {
+ c += 1.0;
+ y += 1.0;
+ z += 2.0;
+ yc = y * c;
+ pk = pkm1 * z - pkm2 * yc;
+ qk = qkm1 * z - qkm2 * yc;
+ if( qk != 0 ) {
+ r = pk/qk;
+ t = Math.abs( (ans - r)/r );
+ ans = r;
+ } else
+ t = 1.0;
+
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+ if( Math.abs(pk) > big ) {
+ pkm2 *= biginv;
+ pkm1 *= biginv;
+ qkm2 *= biginv;
+ qkm1 *= biginv;
+ }
+ } while( t > MACHEP );
+
+ return ans * ax;
+ }
+
+
+ /**
+ * @param a double value
+ * @param x double value
+ * @return the Incomplete Gamma function.
+ *
+ *
+ * Converted to Java from
+ * Cephes Math Library Release 2.2: July, 1992
+ * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ **/
+ static public double igam(double a, double x)
+ throws ArithmeticException {
+
+
+ double ans, ax, c, r;
+
+ if( x <= 0 || a <= 0 ) return 0.0;
+
+ if( x > 1.0 && x > a ) return 1.0 - igamc(a,x);
+
+ /* Compute x**a * exp(-x) / gamma(a) */
+ ax = a * Math.log(x) - x - lgamma(a);
+ if( ax < -MAXLOG ) return( 0.0 );
+
+ ax = Math.exp(ax);
+
+ /* power series */
+ r = a;
+ c = 1.0;
+ ans = 1.0;
+
+ do {
+ r += 1.0;
+ c *= x/r;
+ ans += c;
+ }
+ while( c/ans > MACHEP );
+
+ return( ans * ax/a );
+
+ }
+
+ /**
+ * Returns the area under the left hand tail (from 0 to x)
+ * of the Chi square probability density function with
+ * v degrees of freedom.
+ *
+ * @param df degrees of freedom
+ * @param x double value
+ * @return the Chi-Square function.
+ **/
+
+ static public double chisq(double df, double x)
+ throws ArithmeticException {
+
+ if( x < 0.0 || df < 1.0 ) return 0.0;
+
+ return igam( df/2.0, x/2.0 );
+
+ }
+
+ /**
+ * Returns the area under the right hand tail (from x to
+ * infinity) of the Chi square probability density function
+ * with v degrees of freedom:
+ *
+ * @param df degrees of freedom
+ * @param x double value
+ * @return the Chi-Square function.
+ *
+ **/
+
+ static public double chisqc(double df, double x)
+ throws ArithmeticException {
+
+ if( x < 0.0 || df < 1.0 ) return 0.0;
+
+ return igamc( df/2.0, x/2.0 );
+
+ }
+
+ /**
+ * Returns the sum of the first k terms of the Poisson
+ * distribution.
+ * @param k number of terms
+ * @param x double value
+ */
+
+ static public double poisson(int k, double x)
+ throws ArithmeticException {
+
+
+ if( k < 0 || x < 0 ) return 0.0;
+
+ return igamc((double)(k+1) ,x);
+ }
+
+ /**
+ * Returns the sum of the terms k+1 to infinity of the Poisson
+ * distribution.
+ * @param k start
+ * @param x double value
+ */
+
+ static public double poissonc(int k, double x)
+ throws ArithmeticException {
+
+
+ if( k < 0 || x < 0 ) return 0.0;
+
+ return igam((double)(k+1),x);
+ }
+
+
+
+ /**
+ * @param a double value
+ * @return The area under the Gaussian probability density
+ * function, integrated from minus infinity to x:
+ */
+
+ static public double normal( double a)
+ throws ArithmeticException {
+ double x, y, z;
+
+ x = a * SQRTH;
+ z = Math.abs(x);
+
+ if( z < SQRTH ) y = 0.5 + 0.5 * erf(x);
+ else {
+ y = 0.5 * erfc(z);
+ if( x > 0 ) y = 1.0 - y;
+ }
+
+ return y;
+ }
+
+
+ /**
+ * @param a double value
+ * @return The complementary Error function
+ *
+ *
+ * Converted to Java from
+ * Cephes Math Library Release 2.2: July, 1992
+ * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ */
+
+ static public double erfc(double a)
+ throws ArithmeticException {
+ double x,y,z,p,q;
+
+ double P[] = {
+ 2.46196981473530512524E-10,
+ 5.64189564831068821977E-1,
+ 7.46321056442269912687E0,
+ 4.86371970985681366614E1,
+ 1.96520832956077098242E2,
+ 5.26445194995477358631E2,
+ 9.34528527171957607540E2,
+ 1.02755188689515710272E3,
+ 5.57535335369399327526E2
+ };
+ double Q[] = {
+ //1.0
+ 1.32281951154744992508E1,
+ 8.67072140885989742329E1,
+ 3.54937778887819891062E2,
+ 9.75708501743205489753E2,
+ 1.82390916687909736289E3,
+ 2.24633760818710981792E3,
+ 1.65666309194161350182E3,
+ 5.57535340817727675546E2
+ };
+
+ double R[] = {
+ 5.64189583547755073984E-1,
+ 1.27536670759978104416E0,
+ 5.01905042251180477414E0,
+ 6.16021097993053585195E0,
+ 7.40974269950448939160E0,
+ 2.97886665372100240670E0
+ };
+ double S[] = {
+ //1.00000000000000000000E0,
+ 2.26052863220117276590E0,
+ 9.39603524938001434673E0,
+ 1.20489539808096656605E1,
+ 1.70814450747565897222E1,
+ 9.60896809063285878198E0,
+ 3.36907645100081516050E0
+ };
+
+ if( a < 0.0 ) x = -a;
+ else x = a;
+
+ if( x < 1.0 ) return 1.0 - erf(a);
+
+ z = -a * a;
+
+ if( z < -MAXLOG ) {
+ if( a < 0 ) return( 2.0 );
+ else return( 0.0 );
+ }
+
+ z = Math.exp(z);
+
+ if( x < 8.0 ) {
+ p = polevl( x, P, 8 );
+ q = p1evl( x, Q, 8 );
+ } else {
+ p = polevl( x, R, 5 );
+ q = p1evl( x, S, 6 );
+ }
+
+ y = (z * p)/q;
+
+ if( a < 0 ) y = 2.0 - y;
+
+ if( y == 0.0 ) {
+ if( a < 0 ) return 2.0;
+ else return( 0.0 );
+ }
+
+
+ return y;
+ }
+
+ /**
+ * @param a double value
+ * @return The Error function
+ *
+ *
+ * Converted to Java from
+ * Cephes Math Library Release 2.2: July, 1992
+ * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ */
+
+ static public double erf(double x)
+ throws ArithmeticException {
+ double y, z;
+ double T[] = {
+ 9.60497373987051638749E0,
+ 9.00260197203842689217E1,
+ 2.23200534594684319226E3,
+ 7.00332514112805075473E3,
+ 5.55923013010394962768E4
+ };
+ double U[] = {
+ //1.00000000000000000000E0,
+ 3.35617141647503099647E1,
+ 5.21357949780152679795E2,
+ 4.59432382970980127987E3,
+ 2.26290000613890934246E4,
+ 4.92673942608635921086E4
+ };
+
+ if( Math.abs(x) > 1.0 ) return( 1.0 - erfc(x) );
+ z = x * x;
+ y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
+ return y;
+}
+
+
+
+
+
+ static private double polevl( double x, double coef[], int N )
+ throws ArithmeticException {
+
+ double ans;
+
+ ans = coef[0];
+
+ for(int i=1; i<=N; i++) { ans = ans*x+coef[i]; }
+
+ return ans;
+ }
+
+ static private double p1evl( double x, double coef[], int N )
+ throws ArithmeticException {
+
+ double ans;
+
+ ans = x + coef[0];
+
+ for(int i=1; i 0.5 ) {
+ p += 1.0;
+ z = p - q;
+ }
+ z = q * Math.sin( Math.PI * z );
+ if( z == 0.0 ) throw new
+ ArithmeticException("lgamma: Overflow");
+ z = LOGPI - Math.log( z ) - w;
+ return z;
+ }
+
+ if( x < 13.0 ) {
+ z = 1.0;
+ while( x >= 3.0 ) {
+ x -= 1.0;
+ z *= x;
+ }
+ while( x < 2.0 ) {
+ if( x == 0.0 ) throw new
+ ArithmeticException("lgamma: Overflow");
+ z /= x;
+ x += 1.0;
+ }
+ if( z < 0.0 ) z = -z;
+ if( x == 2.0 ) return Math.log(z);
+ x -= 2.0;
+ p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
+ return( Math.log(z) + p );
+ }
+
+ if( x > 2.556348e305 ) throw new
+ ArithmeticException("lgamma: Overflow");
+
+ q = ( x - 0.5 ) * Math.log(x) - x + 0.91893853320467274178;
+ if( x > 1.0e8 ) return( q );
+
+ p = 1.0/(x*x);
+ if( x >= 1000.0 )
+ q += (( 7.9365079365079365079365e-4 * p
+ - 2.7777777777777777777778e-3) *p
+ + 0.0833333333333333333333) / x;
+ else
+ q += polevl( p, A, 4 ) / x;
+ return q;
+ }
+
+
+
+ /**
+ * @param aa double value
+ * @param bb double value
+ * @param xx double value
+ * @return The Incomplete Beta Function evaluated from zero to xx.
+ *
+ *
+ * Converted to Java from
+ * Cephes Math Library Release 2.3: July, 1995
+ * Copyright 1984, 1995 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ */
+
+ public static double ibeta( double aa, double bb, double xx )
+ throws ArithmeticException {
+ double a, b, t, x, xc, w, y;
+ boolean flag;
+
+ if( aa <= 0.0 || bb <= 0.0 ) throw new
+ ArithmeticException("ibeta: Domain error!");
+
+ if( (xx <= 0.0) || ( xx >= 1.0) ) {
+ if( xx == 0.0 ) return 0.0;
+ if( xx == 1.0 ) return 1.0;
+ throw new ArithmeticException("ibeta: Domain error!");
+ }
+
+ flag = false;
+ if( (bb * xx) <= 1.0 && xx <= 0.95) {
+ t = pseries(aa, bb, xx);
+ return t;
+ }
+
+ w = 1.0 - xx;
+
+ /* Reverse a and b if x is greater than the mean. */
+ if( xx > (aa/(aa+bb)) ) {
+ flag = true;
+ a = bb;
+ b = aa;
+ xc = xx;
+ x = w;
+ } else {
+ a = aa;
+ b = bb;
+ xc = w;
+ x = xx;
+ }
+
+ if( flag && (b * x) <= 1.0 && x <= 0.95) {
+ t = pseries(a, b, x);
+ if( t <= MACHEP ) t = 1.0 - MACHEP;
+ else t = 1.0 - t;
+ return t;
+ }
+
+ /* Choose expansion for better convergence. */
+ y = x * (a+b-2.0) - (a-1.0);
+ if( y < 0.0 )
+ w = incbcf( a, b, x );
+ else
+ w = incbd( a, b, x ) / xc;
+
+ /* Multiply w by the factor
+ a b _ _ _
+ x (1-x) | (a+b) / ( a | (a) | (b) ) . */
+
+ y = a * Math.log(x);
+ t = b * Math.log(xc);
+ if( (a+b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG ) {
+ t = Math.pow(xc,b);
+ t *= Math.pow(x,a);
+ t /= a;
+ t *= w;
+ t *= gamma(a+b) / (gamma(a) * gamma(b));
+ if( flag ) {
+ if( t <= MACHEP ) t = 1.0 - MACHEP;
+ else t = 1.0 - t;
+ }
+ return t;
+ }
+ /* Resort to logarithms. */
+ y += t + lgamma(a+b) - lgamma(a) - lgamma(b);
+ y += Math.log(w/a);
+ if( y < MINLOG )
+ t = 0.0;
+ else
+ t = Math.exp(y);
+
+ if( flag ) {
+ if( t <= MACHEP ) t = 1.0 - MACHEP;
+ else t = 1.0 - t;
+ }
+ return t;
+ }
+
+/* Continued fraction expansion #1
+ * for incomplete beta integral
+ */
+
+ private static double incbcf( double a, double b, double x )
+ throws ArithmeticException {
+ double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
+ double k1, k2, k3, k4, k5, k6, k7, k8;
+ double r, t, ans, thresh;
+ int n;
+ double big = 4.503599627370496e15;
+ double biginv = 2.22044604925031308085e-16;
+
+ k1 = a;
+ k2 = a + b;
+ k3 = a;
+ k4 = a + 1.0;
+ k5 = 1.0;
+ k6 = b - 1.0;
+ k7 = k4;
+ k8 = a + 2.0;
+
+ pkm2 = 0.0;
+ qkm2 = 1.0;
+ pkm1 = 1.0;
+ qkm1 = 1.0;
+ ans = 1.0;
+ r = 1.0;
+ n = 0;
+ thresh = 3.0 * MACHEP;
+ do {
+ xk = -( x * k1 * k2 )/( k3 * k4 );
+ pk = pkm1 + pkm2 * xk;
+ qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ xk = ( x * k5 * k6 )/( k7 * k8 );
+ pk = pkm1 + pkm2 * xk;
+ qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ if( qk != 0 ) r = pk/qk;
+ if( r != 0 ) {
+ t = Math.abs( (ans - r)/r );
+ ans = r;
+ } else
+ t = 1.0;
+
+ if( t < thresh ) return ans;
+
+ k1 += 1.0;
+ k2 += 1.0;
+ k3 += 2.0;
+ k4 += 2.0;
+ k5 += 1.0;
+ k6 -= 1.0;
+ k7 += 2.0;
+ k8 += 2.0;
+
+ if( (Math.abs(qk) + Math.abs(pk)) > big ) {
+ pkm2 *= biginv;
+ pkm1 *= biginv;
+ qkm2 *= biginv;
+ qkm1 *= biginv;
+ }
+ if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
+ pkm2 *= big;
+ pkm1 *= big;
+ qkm2 *= big;
+ qkm1 *= big;
+ }
+ } while( ++n < 300 );
+
+ return ans;
+ }
+/* Continued fraction expansion #2
+ * for incomplete beta integral
+ */
+
+ static private double incbd( double a, double b, double x )
+ throws ArithmeticException {
+ double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
+ double k1, k2, k3, k4, k5, k6, k7, k8;
+ double r, t, ans, z, thresh;
+ int n;
+ double big = 4.503599627370496e15;
+ double biginv = 2.22044604925031308085e-16;
+
+ k1 = a;
+ k2 = b - 1.0;
+ k3 = a;
+ k4 = a + 1.0;
+ k5 = 1.0;
+ k6 = a + b;
+ k7 = a + 1.0;
+ k8 = a + 2.0;
+
+ pkm2 = 0.0;
+ qkm2 = 1.0;
+ pkm1 = 1.0;
+ qkm1 = 1.0;
+ z = x / (1.0-x);
+ ans = 1.0;
+ r = 1.0;
+ n = 0;
+ thresh = 3.0 * MACHEP;
+ do {
+ xk = -( z * k1 * k2 )/( k3 * k4 );
+ pk = pkm1 + pkm2 * xk;
+ qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ xk = ( z * k5 * k6 )/( k7 * k8 );
+ pk = pkm1 + pkm2 * xk;
+ qk = qkm1 + qkm2 * xk;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ if( qk != 0 ) r = pk/qk;
+ if( r != 0 ) {
+ t = Math.abs( (ans - r)/r );
+ ans = r;
+ } else
+ t = 1.0;
+
+ if( t < thresh ) return ans;
+
+ k1 += 1.0;
+ k2 -= 1.0;
+ k3 += 2.0;
+ k4 += 2.0;
+ k5 += 1.0;
+ k6 += 1.0;
+ k7 += 2.0;
+ k8 += 2.0;
+
+ if( (Math.abs(qk) + Math.abs(pk)) > big ) {
+ pkm2 *= biginv;
+ pkm1 *= biginv;
+ qkm2 *= biginv;
+ qkm1 *= biginv;
+ }
+ if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
+ pkm2 *= big;
+ pkm1 *= big;
+ qkm2 *= big;
+ qkm1 *= big;
+ }
+ } while( ++n < 300 );
+
+ return ans;
+ }
+/* Power series for incomplete beta integral.
+ Use when b*x is small and x not too close to 1. */
+
+ static private double pseries( double a, double b, double x )
+ throws ArithmeticException {
+ double s, t, u, v, n, t1, z, ai;
+
+ ai = 1.0 / a;
+ u = (1.0 - b) * x;
+ v = u / (a + 1.0);
+ t1 = v;
+ t = u;
+ n = 2.0;
+ s = 0.0;
+ z = MACHEP * ai;
+ while( Math.abs(v) > z ) {
+ u = (n - b) * x / n;
+ t *= u;
+ v = t / (a + n);
+ s += v;
+ n += 1.0;
+ }
+ s += t1;
+ s += ai;
+
+ u = a * Math.log(x);
+ if( (a+b) < MAXGAM && Math.abs(u) < MAXLOG ) {
+ t = gamma(a+b)/(gamma(a)*gamma(b));
+ s = s * t * Math.pow(x,a);
+ } else {
+ t = lgamma(a+b) - lgamma(a) - lgamma(b) + u + Math.log(s);
+ if( t < MINLOG ) s = 0.0;
+ else s = Math.exp(t);
+ }
+ return s;
+ }
+ public static double getnormpdf(double x) {
+// k = find(sigma > 0);
+// if any(k)
+// xn = (x(k) - mu(k)) ./ sigma(k);
+// y(k) = exp(-0.5 * xn .^2) ./ (sqrt(2*pi) .* sigma(k));
+// end
+ double y = Math.exp(-0.5 *x *x) / Math.sqrt(2.0*Math.PI);
+ return y;
+ }
+ public static double getnormcdf(double x) {
+ //p = 0.5 * erfc(-(x-mu)./(sqrt(2)*sigma));
+ double p = 0.5 * erfc(-x/Math.sqrt(2));
+ return p;
+ }
+
+}
diff --git a/src/wsi/ra/math/interpolation/AbstractDataModifier.java b/src/wsi/ra/math/interpolation/AbstractDataModifier.java
new file mode 100644
index 00000000..d2316efe
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/AbstractDataModifier.java
@@ -0,0 +1,52 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: AbstractDataModifier.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:24:58 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+/*==========================================================================*
+* CLASS DECLARATION
+*==========================================================================*/
+
+/**
+ * The minimal set of functions which should implemented in a data modifier for
+ * AbstractDataSet.
+ */
+public abstract class AbstractDataModifier
+{
+
+ /*-------------------------------------------------------------------------*
+ * abstract methods
+ *-------------------------------------------------------------------------*/
+
+ /**
+ * Modifies the X data.
+ */
+ public abstract void modifyX(double[] setX);
+ /**
+ * Modifies the Y data.
+ */
+ public abstract void modifyY(double[] setY);
+ /**
+ * Modifies the data.
+ */
+ public abstract void modify(double[] setX, double[] setY);
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/math/interpolation/AbstractDataSet.java b/src/wsi/ra/math/interpolation/AbstractDataSet.java
new file mode 100644
index 00000000..633d9d2f
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/AbstractDataSet.java
@@ -0,0 +1,116 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: AbstractDataSet.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:04 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+/*==========================================================================*
+* CLASS DECLARATION
+*==========================================================================*/
+
+public abstract class AbstractDataSet
+{
+ /*-------------------------------------------------------------------------*
+ * public member variables
+ *-------------------------------------------------------------------------*/
+
+ /*--------------------------------------------------------------o-----------*
+ * protected member variables
+ *-------------------------------------------------------------------------*/
+ /**
+ * double array of X data.
+ *
+ * @see #yDoubleData
+ */
+ protected double[] xDoubleData = { -1, 1 };
+ /**
+ * double array of Y data.
+ *
+ * @see #xDoubleData
+ */
+ protected double[] yDoubleData = { 1, 1 };
+
+ /*-------------------------------------------------------------------------*
+ * abstract methods
+ *-------------------------------------------------------------------------*/
+
+ /**
+ * Returns the length of the data set
+ * @return the length of the data set
+ */
+ public int getLength()
+ {
+ return xDoubleData.length;
+ }
+ /**
+ * Returns an array of the X data
+ * @return the array of the X data
+ */
+ public double[] getXData()
+ {
+ return xDoubleData;
+ }
+ /**
+ * Returns an array of the Y data
+ * @return the array of the Y data
+ */
+ public double[] getYData()
+ {
+ return yDoubleData;
+ }
+ /**
+ * Returns the X label of the data set
+ * @return the X label of the data set
+ */
+ public abstract String getXLabel();
+ /**
+ * Returns the Y label of the data set
+ * @return the Y label of the data set
+ */
+ public abstract String getYLabel();
+
+ /**
+ * Modifies the X data.
+ *
+ * @param the data modifier
+ */
+ public void modifyXData(AbstractDataModifier modifier)
+ {
+ modifier.modifyX(xDoubleData);
+ }
+ /**
+ * Modifies the Y data.
+ *
+ * @param the data modifier
+ */
+ public void modifyYData(AbstractDataModifier modifier)
+ {
+ modifier.modifyY(yDoubleData);
+ }
+ /**
+ * Modifies the data.
+ *
+ * @param the data modifier
+ */
+ public void modifyData(AbstractDataModifier modifier)
+ {
+ modifier.modify(xDoubleData, yDoubleData);
+ }
+}
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
diff --git a/src/wsi/ra/math/interpolation/BasicDataSet.java b/src/wsi/ra/math/interpolation/BasicDataSet.java
new file mode 100644
index 00000000..3e5eed0a
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/BasicDataSet.java
@@ -0,0 +1,100 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: BasicDataSet.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:11 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+import wsi.ra.sort.XYDoubleArray;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+/**
+ * The minimum wrapper class for an AbstractDataSet.
+ */
+public class BasicDataSet extends AbstractDataSet
+{
+ /*-------------------------------------------------------------------------*
+ * protected member variables
+ *-------------------------------------------------------------------------*/
+ protected int dataType = -1;
+ protected String xLabel = null;
+ protected String yLabel = null;
+
+ /*------------------------------------------------------------------------*
+ * constructor
+ *------------------------------------------------------------------------*/
+ public BasicDataSet()
+ {
+ this(null, null, null, null);
+ }
+
+ public BasicDataSet(XYDoubleArray data)
+ {
+ this(data.x, data.y, null, null);
+ }
+
+ public BasicDataSet(XYDoubleArray data, String xLabel, String yLabel)
+ {
+ this(data.x, data.y, xLabel, yLabel);
+ }
+
+ public BasicDataSet(
+ double[] xDoubleData,
+ double[] yDoubleData,
+ int dataType)
+ {
+ this(xDoubleData, yDoubleData, null, null);
+ }
+
+ public BasicDataSet(
+ double[] xDoubleData,
+ double[] yDoubleData,
+ String xLabel,
+ String yLabel)
+ {
+ this.xDoubleData = xDoubleData;
+ this.yDoubleData = yDoubleData;
+ this.xLabel = xLabel;
+ this.yLabel = yLabel;
+ }
+
+ /*-------------------------------------------------------------------------*
+ * public methods
+ *-------------------------------------------------------------------------*/
+
+ public int getDataType()
+ {
+ return dataType;
+ }
+
+ public String getXLabel()
+ {
+ return xLabel;
+ }
+
+ public String getYLabel()
+ {
+ return yLabel;
+ }
+
+ public String getAdditionalInformation(String parm1)
+ {
+ return new String();
+ }
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/math/interpolation/InterpolationException.java b/src/wsi/ra/math/interpolation/InterpolationException.java
new file mode 100644
index 00000000..6f6015ef
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/InterpolationException.java
@@ -0,0 +1,44 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: InterpolationException.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:17 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+
+/*==========================================================================*
+ * CLASS DECLARATION
+ *==========================================================================*/
+
+/**
+ * Exception for interpolation error.
+ */
+public class InterpolationException extends Exception
+{
+
+ public InterpolationException()
+ {
+ super();
+ }
+ public InterpolationException(String s)
+ {
+ super(s);
+ }
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/math/interpolation/LinearInterpolation.java b/src/wsi/ra/math/interpolation/LinearInterpolation.java
new file mode 100644
index 00000000..ab04fadb
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/LinearInterpolation.java
@@ -0,0 +1,573 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: LinearInterpolation.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Original author: Charles S. Stanton
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:23 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+//import cern.jet.stat.*;
+
+/**
+ * Defines the routines for the spline interpolation of data.
+ */
+public class LinearInterpolation
+{
+ AbstractDataSet abstractDataSet = null;
+ private double[] x, y;
+ //vectors of x,y
+ private double sumX = 0;
+ private double sumY = 0;
+ private double sumXY = 0;
+ private double sumXsquared = 0;
+ private double sumYsquared = 0;
+ private double Sxx, Sxy, Syy, n;
+ private double a = 0, b = 0;
+ //coefficients of regression
+ private int dataLength;
+ private double[][] residual;
+ // residual[0][i] = x[i], residual[1][i]= residual
+ private double maxAbsoluteResidual = 0.0;
+ private double SSR = 0.0;
+ //regression sum of squares
+ private double SSE = 0.0;
+ //error sum of squares
+ private double minX = Double.POSITIVE_INFINITY;
+ private double maxX = Double.NEGATIVE_INFINITY;
+ private double minY = Double.POSITIVE_INFINITY;
+ private double maxY = Double.NEGATIVE_INFINITY;
+
+ //MISC
+ String xName, yName;
+ double aCILower, aCIUpper, bCILower, bCIUpper; //confidence interval
+ double t, bSE, aSE;
+ double MSE, F;
+ static double[] t025 =
+ {
+ Double.NaN,
+ 12.706,
+ 4.303,
+ 3.182,
+ 2.776,
+ 2.571,
+ 2.447,
+ 2.365,
+ 2.306,
+ 2.262,
+ 2.228,
+ 2.201,
+ 2.179,
+ 2.160,
+ 2.145,
+ 2.131,
+ 2.120,
+ 2.110,
+ 2.101,
+ 2.093,
+ 2.086,
+ 2.080,
+ 2.075,
+ 2.069,
+ 2.064,
+ 2.060,
+ 2.056,
+ 2.052,
+ 2.048,
+ 2.045,
+ 1.960 };
+
+ /*------------------------------------------------------------------------*
+ * constructor
+ *------------------------------------------------------------------------*/
+ /**
+ * Initializes this class.
+ */
+ public LinearInterpolation() throws InterpolationException
+ {
+ this.abstractDataSet = null;
+ }
+
+ /**
+ * Constructor for regression calculator.
+ *
+ * @param x is the array of x data
+ * @param y is the array of y data
+ */
+ public LinearInterpolation(double[] x, double[] y)
+ {
+ this.x = x;
+ this.y = y;
+ if (x.length != y.length)
+ {
+ System.out.println("x, y vectors must be of same length");
+ }
+ else
+ {
+ dataLength = x.length;
+ doStatistics();
+ }
+
+ }
+
+ /**
+ * Initializes this class and calculates the second derivative of the spline.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ */
+ public LinearInterpolation(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ this.setAbstractDataSet(abstractDataSet);
+ }
+
+ public void setAbstractDataSet(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ this.abstractDataSet = abstractDataSet;
+ x = abstractDataSet.getXData();
+ y = abstractDataSet.getYData();
+ if (x.length != y.length)
+ {
+ System.out.println("x, y vectors must be of same length");
+ }
+ else
+ {
+ dataLength = x.length;
+ doStatistics();
+ }
+ }
+
+ /**
+ * Find the p value for a given value of F.
+ * Requires the COLT high performance library:
+ * http://hoschek.home.cern.ch/hoschek/colt/
+ *
+ * @param fValue the value for the CDF
+ * @return The P value
+ */
+ // public double getP(double fValue) {
+ // double answer;
+ // double y1;
+ // double y2;
+ // //nu1 = 1;
+ // //x2 =1
+ // double nu2 = n - 2;
+ // y1 = nu2 / (nu2 + fValue);
+ // y2 = 0.0;
+ // answer = Gamma.incompleteBeta(nu2 / 2.0, 1 / 2.0, y1)
+ // - Gamma.incompleteBeta(nu2 / 2.0, 1 / 2.0, y2);
+ // return answer;
+ // }
+
+ /*
+ * Here are the accessor methods
+ *
+ */
+ /**
+ * Gets the intercept of the regression line.
+ *
+ * @return The intercept.
+ */
+ public double getIntercept()
+ {
+ return a;
+ }
+
+ /**
+ * Gets the Slope of the regression line.
+ *
+ * @return The slope.
+ */
+ public double getSlope()
+ {
+ return b;
+ }
+
+ /**
+ * Gets the residuals of the regression.
+ *
+ * @return The residuals.
+ */
+ public double[][] getResiduals()
+ {
+ return residual;
+ }
+
+ /**
+ * Gets the x data for the regression.
+ *
+ * @return The array of x values.
+ */
+ public double[] getDataX()
+ {
+ return x;
+ }
+
+ /**
+ * Gets the y data for the regression.
+ *
+ * @return The array of y values.
+ */
+ public double[] getDataY()
+ {
+ return y;
+ }
+
+ /**
+ * Gets the minimum of the x values.
+ *
+ * @return The minimum.
+ */
+ public double getMinX()
+ {
+ return minX;
+ }
+
+ /**
+ * Gets the maximum of the x values.
+ *
+ * @return The maximum.
+ */
+ public double getMaxX()
+ {
+ return maxX;
+ }
+
+ /**
+ * Gets the minimum of the y values.
+ *
+ * @return The minumum.
+ */
+ public double getMinY()
+ {
+ return minY;
+ }
+
+ /**
+ * Gets the maximum of the y values.
+ *
+ * @return The maximum.
+ */
+ public double getMaxY()
+ {
+ return maxY;
+ }
+
+ /**
+ * Gets the maximum absolute residual.
+ *
+ * @return The maximum.
+ */
+ public double getMaxAbsoluteResidual()
+ {
+ return maxAbsoluteResidual;
+ }
+
+ /**
+ * Gets the sum of the square x deviations from mean of x.
+ *
+ * @return The Sxx value
+ */
+ public double getSxx()
+ {
+ return Sxx;
+ }
+
+ /**
+ * Gets the sum of the square y deviations from mean of y.
+ *
+ * @return The Syy value
+ */
+ public double getSyy()
+ {
+ return Syy;
+ }
+
+ /**
+ * Gets SSR = Sxy * Sxy / Sxx;
+ *
+ * @return The SSR value
+ */
+ public double getSSR()
+ {
+ return SSR;
+ }
+
+ /**
+ * Gets SSE = Syy - SSR.
+ *
+ * @return The SSE value
+ */
+ public double getSSE()
+ {
+ return SSE;
+ }
+
+ /**
+ * Gets the mean square error MSE.
+ *
+ * @return The MSE value
+ */
+ public double getMSE()
+ {
+ return SSE / (n - 2);
+ }
+
+ /**
+ * Gets the mean XBar of x.
+ *
+ * @return The XBar value
+ */
+ public double getXBar()
+ {
+ return sumX / n;
+ }
+
+ /**
+ * Gets the mean YBar of y.
+ *
+ * @return The YBar value
+ */
+ public double getYBar()
+ {
+ return sumY / n;
+ }
+
+ /**
+ * Gets the sample size.
+ *
+ * @return The sample size.
+ */
+ public int getDataLength()
+ {
+ return x.length;
+ }
+
+ /**
+ * Gets the Pearson R statistic of the regression.
+ *
+ * @return The PearsonR value
+ */
+ public double getPearsonR()
+ {
+ return Sxy / Math.sqrt(Sxx * Syy);
+ }
+
+ /**
+ * Gets the sum of the x squared values.
+ *
+ * @return The sum of the x squared values.
+ */
+ public double getSumXSquared()
+ {
+ return sumXsquared;
+ }
+
+ /**
+ * reset data to 0
+ */
+ public void reset()
+ {
+ x = new double[0];
+ y = new double[0];
+ dataLength = 0;
+ n = 0.0;
+ residual = new double[0][0];
+
+ sumX = 0;
+ sumXsquared = 0;
+ sumY = 0;
+ sumYsquared = 0;
+ sumXY = 0;
+
+ }
+
+ /**
+ * Adds a new point to the regression (for interactive use).
+ *
+ * @param xValue The new x value
+ * @param yValue The new y value
+ */
+ public void addPoint(double xValue, double yValue)
+ {
+ dataLength++;
+ double[] xNew = new double[dataLength];
+ double[] yNew = new double[dataLength];
+ System.arraycopy(x, 0, xNew, 0, dataLength - 1);
+ System.arraycopy(y, 0, yNew, 0, dataLength - 1);
+ xNew[dataLength - 1] = xValue;
+ yNew[dataLength - 1] = yValue;
+ x = xNew;
+ y = yNew;
+ updateStatistics(xValue, yValue);
+ }
+
+ private void doStatistics()
+ {
+ //Find sum of squares for x,y and sum of xy
+ for (int i = 0; i < dataLength; i++)
+ {
+ minX = Math.min(minX, x[i]);
+ maxX = Math.max(maxX, x[i]);
+ minY = Math.min(minY, y[i]);
+ maxY = Math.max(maxY, y[i]);
+ sumX += x[i];
+ sumY += y[i];
+ sumXsquared += x[i] * x[i];
+ sumYsquared += y[i] * y[i];
+ sumXY += x[i] * y[i];
+ }
+ //Caculate regression coefficients
+ n = (double) dataLength;
+ Sxx = sumXsquared - sumX * sumX / n;
+ Syy = sumYsquared - sumY * sumY / n;
+ Sxy = sumXY - sumX * sumY / n;
+ b = Sxy / Sxx;
+ a = (sumY - b * sumX) / n;
+ SSR = Sxy * Sxy / Sxx;
+ SSE = Syy - SSR;
+ calculateResiduals();
+ }
+
+ private void calculateResiduals()
+ {
+ residual = new double[2][dataLength];
+ maxAbsoluteResidual = 0.0;
+ for (int i = 0; i < dataLength; i++)
+ {
+ residual[0][i] = x[i];
+ residual[1][i] = y[i] - (a + b * x[i]);
+ maxAbsoluteResidual =
+ Math.max(maxAbsoluteResidual, Math.abs(y[i] - (a + b * x[i])));
+ }
+ }
+
+ //update statistics for a single additional data point
+ private void updateStatistics(double xValue, double yValue)
+ {
+ //Find sum of squares for x,y and sum of xy
+ n++;
+ sumX += xValue;
+ sumY += yValue;
+ sumXsquared += xValue * xValue;
+ sumYsquared += yValue * yValue;
+ sumXY += xValue * yValue;
+ //Caculate regression coefficients
+ n = (double) dataLength;
+ Sxx = sumXsquared - sumX * sumX / n;
+ Syy = sumYsquared - sumY * sumY / n;
+ Sxy = sumXY - sumX * sumY / n;
+ b = Sxy / Sxx;
+ a = (sumY - b * sumX) / n;
+ SSR = Sxy * Sxy / Sxx;
+ SSE = Syy - SSR;
+ calculateResiduals();
+ }
+
+ /**
+ * regression line y = a + bx.
+ *
+ * @param x
+ * @return double
+ * @throws InterpolationException
+ */
+ public double getY(double x) throws InterpolationException
+ {
+ return a + b * x;
+ }
+
+ public String toString()
+ {
+ StringBuffer sb = new StringBuffer(1000);
+
+ sb.append("Regression Statistics for " + yName + " = a + b*" + xName);
+ sb.append("");
+ sb.append("Sample Statistics");
+ int n = this.getDataLength();
+ sb.append("Sample size n = " + n);
+ sb.append("Mean of " + yName + " Y bar = " + this.getYBar());
+ sb.append("s_Y");
+ sb.append("= " + Math.sqrt(this.getSyy() / ((float) (n - 1))));
+ sb.append("Pearson correlation R = " + this.getPearsonR());
+ sb.append("");
+ sb.append("Coefficient Estimates");
+ a = this.getIntercept();
+ b = this.getSlope();
+ sb.append("a = " + a);
+ sb.append("b = " + b);
+ sb.append("");
+
+ sb.append("95% Confidence Intervals");
+
+ if (n > 32)
+ {
+ t = t025[30];
+ }
+ else if (n > 2)
+ {
+ t = t025[n - 2];
+ }
+ else
+ {
+ t = Double.NaN;
+ }
+ MSE = this.getMSE();
+ if (n > 2)
+ {
+ bSE = Math.sqrt(MSE / this.getSxx());
+ }
+ else
+ {
+ bSE = Double.NaN;
+ }
+ aSE = bSE * Math.sqrt(this.getSumXSquared() / n);
+ aCILower = a - t * aSE;
+ aCIUpper = a + t * aSE;
+ sb.append("a : (" + aCILower + ", " + aCIUpper + ")");
+ bCILower = b - t * bSE;
+ bCIUpper = b + t * bSE;
+ sb.append("b : (" + bCILower + ", " + bCIUpper + ")");
+ sb.append("");
+ sb.append("Analysis of Variance");
+ sb.append("Source Degrees Freedom Sum of Squares");
+ sb.append("");
+ SSR = this.getSSR();
+ //allow one degree of freedom for mean
+ sb.append(
+ "model 1 "
+ + SSR);
+ sb.append(
+ "error "
+ + (n - 2)
+ + " "
+ + this.getSSE());
+ sb.append(
+ "total(corrected) "
+ + (n - 1)
+ + " "
+ + this.getSyy());
+ sb.append("");
+ sb.append("MSE =" + MSE);
+ F = SSR / MSE;
+ sb.append("F = " + F + " ");
+ //sb.append("p = " + this.getP(F));
+
+ return sb.toString();
+ }
+}
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/math/interpolation/PolyInterpolation.java b/src/wsi/ra/math/interpolation/PolyInterpolation.java
new file mode 100644
index 00000000..868d7947
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/PolyInterpolation.java
@@ -0,0 +1,453 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: PolyInterpolation.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:30 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+/**
+ * Defines the routines for the interpolation of data.
+ */
+public class PolyInterpolation
+{
+ AbstractDataSet abstractDataSet = null;
+ boolean sloppy = true;
+ double[] polynomialCoefficients = null;
+
+ /*------------------------------------------------------------------------*
+ * constructor
+ *------------------------------------------------------------------------*/
+ /**
+ * Initializes this class.
+ */
+ public PolyInterpolation() throws InterpolationException
+ {
+ this.abstractDataSet = null;
+ sloppy = true;
+ polynomialCoefficients = null;
+ }
+
+ /**
+ * Initializes this class and calculates the coefficients of the polynom.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ */
+ public PolyInterpolation(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ this.abstractDataSet = abstractDataSet;
+ sloppy = true;
+ this.polynomialCoefficients = calculatePolynomialCoefficients();
+ }
+
+ /**
+ * Initializes this class and calculates the coefficients of the polynom.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ * @param sloppy if true Neville's algorithm which is used in the
+ * polynomialInterpolation-routines does only print a
+ * warning message on the screen and does not throw an
+ * Exception if two x values are identical.
+ */
+ public PolyInterpolation(AbstractDataSet abstractDataSet, boolean sloppy)
+ throws InterpolationException
+ {
+ this.abstractDataSet = abstractDataSet;
+ this.sloppy = sloppy;
+ this.polynomialCoefficients = calculatePolynomialCoefficients();
+ }
+
+ /**
+ * Sets the new AbstractDataSet and calculates the coefficients
+ * of the polynom.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ */
+ public void setAbstractDataSet(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ this.abstractDataSet = abstractDataSet;
+ this.polynomialCoefficients = calculatePolynomialCoefficients();
+ }
+
+ /**
+ * Uses the polynom with the calculated coefficients to calculate the
+ * y value. This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 5, pages 173-176.
+ * The Neville's algorithm which is used in the polynomialInterpolation-
+ * routines returns also the error of this interpolated point.
+ *
+ * @param x the x value
+ * @return the interpolated y value
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #getYandDerivatives(double, int)
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ */
+
+ public double getY(double x)
+ {
+ int n = polynomialCoefficients.length - 1;
+ double y = polynomialCoefficients[n];
+ for (int j = n - 1; j >= 0; j--)
+ y = y * x + polynomialCoefficients[j];
+
+ return y;
+ }
+
+ /**
+ * Uses the polynom with the calculated coefficients to calculate the
+ * y value and the derivatives at the point x,
+ * y. This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 5, pages 173-176.
+ * The Neville's algorithm which is used in the polynomialInterpolation-
+ * routines returns also the error of this interpolated point.
+ *
+ * @param x the x value
+ * @param ndDerivateNumber the number of the calculated derivatives
+ * @return the interpolated y value at ...[0], the 1st derivativ value at
+ * ...[1], the 2nd derivativ at ...[2] and so on ...
+ * @see #getY(double)
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ */
+ public double[] getYandDerivatives(double x, int ndDerivateNumber)
+ throws InterpolationException
+ {
+ if (ndDerivateNumber < 0)
+ throw new InterpolationException("Negative derivative numbers make no sense.");
+ else if (ndDerivateNumber == 0)
+ {
+ double[] pd = new double[1];
+ pd[0] = getY(x);
+ return pd;
+ }
+
+ int nnd, j, i;
+ int nc = polynomialCoefficients.length - 1;
+ double[] pd = new double[ndDerivateNumber + 1];
+ double cnst = 1.0;
+
+ pd[0] = polynomialCoefficients[nc];
+ for (j = 1; j <= ndDerivateNumber; j++)
+ pd[j] = 0.0;
+ for (i = nc - 1; i >= 0; i--)
+ {
+ nnd = (ndDerivateNumber < (nc - i) ? ndDerivateNumber : nc - i);
+ for (j = nnd; j >= 1; j--)
+ pd[j] = pd[j] * x + pd[j - 1];
+ pd[0] = pd[0] * x + polynomialCoefficients[i];
+ }
+ for (i = 2; i <= ndDerivateNumber; i++)
+ {
+ cnst *= i;
+ pd[i] *= cnst;
+ }
+
+ return pd;
+ }
+
+ /**
+ * Neville's interpolation algorithm. This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 108-122.
+ *
+ * @param x the x value
+ * @return the interpolated y value and the interpolation error
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #getY(double)
+ * @see #getYandDerivatives(double, int)
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ */
+ public PolynomialInterpolationResult polynomialInterpolation(double x)
+ throws InterpolationException
+ {
+ if (abstractDataSet == null)
+ throw new InterpolationException(
+ "No data." + " The AbstractDataSet was not defined.");
+ return polynomialInterpolation(
+ abstractDataSet.getXData(),
+ abstractDataSet.getYData(),
+ x);
+ }
+
+ /**
+ * Neville's interpolation algorithm. This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 108-122.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ * @param x the x value
+ * @return the interpolated y value and the interpolation error
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #getY(double)
+ * @see #getYandDerivatives(double, int)
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ */
+ public PolynomialInterpolationResult polynomialInterpolation(
+ AbstractDataSet abstractDataSet,
+ double x)
+ throws InterpolationException
+ {
+ if (abstractDataSet == null)
+ throw new InterpolationException(
+ "No data." + " The AbstractDataSet was not defined.");
+ return polynomialInterpolation(
+ abstractDataSet.getXData(),
+ abstractDataSet.getYData(),
+ x);
+ }
+
+ /**
+ * Neville's interpolation algorithm. This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 108-122.
+ *
+ * @param xa the array of x values
+ * @param ya the array of y values
+ * @param x the x value
+ * @return the interpolated y value and the interpolation error
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #getY(double)
+ * @see #getYandDerivatives(double, int)
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ */
+ public PolynomialInterpolationResult polynomialInterpolation(
+ double[] xa,
+ double[] ya,
+ double x)
+ throws InterpolationException
+ {
+ if (xa == null || ya == null)
+ throw new InterpolationException("No data.");
+ int i, m, ns = 1;
+ double den, dif, dift, ho, hp, w;
+ double[] c = new double[xa.length + 1];
+ double[] d = new double[xa.length + 1];
+ PolynomialInterpolationResult result =
+ new PolynomialInterpolationResult();
+
+ dif = Math.abs(x - xa[1 - 1]);
+ for (i = 1; i <= xa.length; i++)
+ {
+ if ((dift = Math.abs(x - xa[i - 1])) < dif)
+ {
+ ns = i;
+ dif = dift;
+ }
+ c[i] = ya[i - 1];
+ d[i] = ya[i - 1];
+ //System.out.println("x"+xa[i-1]+" y"+ya[i-1]);
+ }
+ result.y = ya[ns - 1];
+ //System.out.println("y="+result.y+" ns="+ns);
+ ns--;
+ for (m = 1; m < xa.length; m++)
+ {
+ for (i = 1; i <= xa.length - m; i++)
+ {
+ ho = xa[i - 1] - x;
+ hp = xa[i + m - 1] - x;
+ w = c[i + 1] - d[i];
+ if ((den = ho - hp) == 0.0)
+ {
+ if (sloppy)
+ {
+ System.out.println(
+ "Two identical x values. The values must be distinct.");
+ den = 1.0;
+ }
+ else
+ throw new InterpolationException("Two identical x values.");
+ }
+ den = w / den;
+ d[i] = hp * den;
+ c[i] = ho * den;
+ }
+ result.y
+ += (result.yError =
+ (2 * ns < (xa.length - m) ? c[ns + 1] : d[ns--]));
+ }
+ return result;
+ }
+
+ /**
+ * Calculates the coefficients of a polynom of the grade N-1. This
+ * interpolation algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 108-122.
+ *
+ * @return the array with the polynomial coefficients y = ...[0] +
+ * ...[1]*x2 + ...[2]*x3 + ...
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #getY(double)
+ * @see #getYandDerivatives(double, int)
+ */
+ public double[] calculatePolynomialCoefficients()
+ throws InterpolationException
+ {
+ if (abstractDataSet == null)
+ throw new InterpolationException(
+ "No data." + " The AbstractDataSet was not defined.");
+ return calculatePolynomialCoefficients(
+ abstractDataSet.getXData(),
+ abstractDataSet.getYData());
+ }
+
+ /**
+ * Calculates the coefficients of a polynom of the grade N-1. This
+ * interpolation algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 108-122.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ * @return the array with the polynomial coefficients y = ...[0] +
+ * ...[1]*x2 + ...[2]*x3 + ...
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(double[], double[])
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #getY(double)
+ * @see #getYandDerivatives(double, int)
+ */
+ public double[] calculatePolynomialCoefficients(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ if (abstractDataSet == null)
+ throw new InterpolationException(
+ "No data." + " The AbstractDataSet was not defined.");
+ return calculatePolynomialCoefficients(
+ abstractDataSet.getXData(),
+ abstractDataSet.getYData());
+ }
+
+ /**
+ * Calculates the coefficients of a polynom of the grade N-1. This
+ * interpolation algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 108-122.
+ *
+ * @param x the array of x values
+ * @param y the array of y values
+ * @return the array with the polynomial coefficients y = ...[0] +
+ * ...[1]*x2 + ...[2]*x3 + ...
+ * @see #calculatePolynomialCoefficients()
+ * @see #calculatePolynomialCoefficients(AbstractDataSet)
+ * @see #polynomialInterpolation(double)
+ * @see #polynomialInterpolation(AbstractDataSet, double)
+ * @see #polynomialInterpolation(double[], double[], double)
+ * @see #getY(double)
+ * @see #getYandDerivatives(double, int)
+ */
+ public double[] calculatePolynomialCoefficients(double x[], double y[])
+ {
+ int k, j, i, n = x.length - 1;
+ double phi, ff, b;
+ double[] s = new double[n + 1];
+ double[] cof = new double[n + 1];
+
+ for (i = 0; i <= n; i++)
+ {
+ s[i] = cof[i] = 0.0;
+ }
+ s[n] = -x[0];
+
+ for (i = 1; i <= n; i++)
+ {
+ for (j = n - i; j <= n - 1; j++)
+ {
+ s[j] -= x[i] * s[j + 1];
+ }
+ s[n] -= x[i];
+ }
+
+ for (j = 0; j < n; j++)
+ {
+ phi = n + 1;
+ for (k = n; k >= 1; k--)
+ {
+ phi = k * s[k] + x[j] * phi;
+ }
+ ff = y[j] / phi;
+ b = 1.0;
+ for (k = n; k >= 0; k--)
+ {
+ cof[k] += b * ff;
+ b = s[k] + x[j] * b;
+ }
+ }
+ return cof;
+ }
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/math/interpolation/PolynomialInterpolationResult.java b/src/wsi/ra/math/interpolation/PolynomialInterpolationResult.java
new file mode 100644
index 00000000..d63dabb2
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/PolynomialInterpolationResult.java
@@ -0,0 +1,56 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: PolynomialInterpolationResult.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:36 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+/*==========================================================================*
+ * CLASS DECLARATION
+ *==========================================================================*/
+
+/**
+ * The result data for a polynomial interpolation.
+ */
+public class PolynomialInterpolationResult
+{
+ /*-------------------------------------------------------------------------*
+ * public member variables
+ *-------------------------------------------------------------------------*/
+
+ public double y = Double.NaN;
+ public double yError = Double.NaN;
+
+ /*------------------------------------------------------------------------*
+ * constructor
+ *------------------------------------------------------------------------*/
+
+ public PolynomialInterpolationResult()
+ {
+ y = Double.NaN;
+ yError = Double.NaN;
+ }
+
+ public PolynomialInterpolationResult(double y, double yError)
+ {
+ this.y = y;
+ this.yError = yError;
+ }
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/math/interpolation/SplineInterpolation.java b/src/wsi/ra/math/interpolation/SplineInterpolation.java
new file mode 100644
index 00000000..191de54d
--- /dev/null
+++ b/src/wsi/ra/math/interpolation/SplineInterpolation.java
@@ -0,0 +1,299 @@
+///////////////////////////////////////////////////////////////////////////////
+// Filename: $RCSfile: SplineInterpolation.java,v $
+// Purpose: Some interpolation stuff.
+// Language: Java
+// Compiler: JDK 1.4
+// Authors: Joerg K. Wegner
+// Version: $Revision: 1.1 $
+// $Date: 2003/07/22 19:25:42 $
+// $Author: wegnerj $
+//
+// Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+//
+///////////////////////////////////////////////////////////////////////////////
+
+package wsi.ra.math.interpolation;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+/**
+ * Defines the routines for the spline interpolation of data.
+ */
+public class SplineInterpolation
+{
+ AbstractDataSet abstractDataSet = null;
+ double[] secondDerivative = null;
+ double[] xArray = null;
+ double[] yArray = null;
+ boolean ascendingData = true;
+
+ /*------------------------------------------------------------------------*
+ * constructor
+ *------------------------------------------------------------------------*/
+ /**
+ * Initializes this class.
+ */
+ public SplineInterpolation() throws InterpolationException
+ {
+ this.abstractDataSet = null;
+ this.secondDerivative = null;
+ }
+
+ /**
+ * Initializes this class and calculates the second derivative of the spline.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ */
+ public SplineInterpolation(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ this.setAbstractDataSet(abstractDataSet);
+ }
+
+ /**
+ * Sets the new AbstractDataSet and calculates the second
+ * derivative of the spline.
+ *
+ * @param abstractDataSet the AbstractDataSet
+ */
+ public void setAbstractDataSet(AbstractDataSet abstractDataSet)
+ throws InterpolationException
+ {
+ this.abstractDataSet = abstractDataSet;
+ double[] x = abstractDataSet.getXData();
+ double[] y = abstractDataSet.getYData();
+ boolean ascending = false;
+ boolean descending = false;
+ int n = x.length;
+
+ xArray = new double[n];
+ yArray = new double[n];
+ xArray[n - 1] = x[0];
+ yArray[n - 1] = y[0];
+ for (int i = 0; i < n - 1; i++)
+ {
+ xArray[i] = x[n - i - 1];
+ yArray[i] = y[n - i - 1];
+ if (x[i] < x[i + 1])
+ {
+ //if(descending)throw new InterpolationException("The x values must be"+
+ // " in continous ascending/descending order.");
+ ascending = true;
+ }
+ else
+ {
+ //if(ascending)throw new InterpolationException("The x values must be"+
+ // " in continous ascending/descending order.");
+ descending = true;
+ }
+ }
+ ascendingData = ascending;
+
+ if (ascendingData)
+ {
+ xArray = null;
+ yArray = null;
+ xArray = abstractDataSet.getXData();
+ yArray = abstractDataSet.getYData();
+ }
+ this.secondDerivative =
+ spline(
+ xArray,
+ yArray,
+ (yArray[1] - yArray[0]) / (xArray[1] - xArray[0]),
+ (yArray[n - 1] - yArray[n - 2]) / (xArray[1] - xArray[n - 2]));
+ }
+
+ /**
+ * Uses the spline with the calculated second derivative values to calculate
+ * the y value. This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 5, pages 173-176.
+ */
+ public double getY(double x) throws InterpolationException
+ {
+ return splineInterpolation(xArray, yArray, secondDerivative, x);
+ }
+
+ public double getDerivative(double x) throws InterpolationException
+ {
+ return splineInterpolatedDerivative(
+ xArray,
+ yArray,
+ secondDerivative,
+ x);
+ }
+
+ /**
+ * Calculates the second derivative of the data. It's important that the
+ * xi values of the function yi=f(xi) are
+ * in ascending order, as x0<x1<x2<... .
+ * This algorithm was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 113-116.
+ */
+ public double[] spline(double[] x, double[] y, double yp0, double ypn)
+ throws InterpolationException
+ {
+ if (x[0] > x[1])
+ throw new InterpolationException(
+ "The x values must be" + " in ascending order.");
+ int n = x.length;
+ double[] y2 = new double[n];
+ double[] u = new double[n - 1];
+ int i, k;
+ double p, qn, sig, un;
+
+ if (yp0 > 0.99e30)
+ y2[0] = u[0] = 0.0;
+ else
+ {
+ y2[0] = -0.5;
+ u[0] =
+ (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp0);
+ }
+ for (i = 2; i <= n - 1; i++)
+ {
+ sig = (x[i - 1] - x[i - 2]) / (x[i] - x[i - 2]);
+ p = sig * y2[i - 2] + 2.0;
+ y2[i - 1] = (sig - 1.0) / p;
+ u[i - 1] =
+ (y[i] - y[i - 1]) / (x[i] - x[i - 1])
+ - (y[i - 1] - y[i - 2]) / (x[i - 1] - x[i - 2]);
+ u[i - 1] =
+ (6.0 * u[i - 1] / (x[i] - x[i - 2]) - sig * u[i - 2]) / p;
+ }
+ if (ypn > 0.99e30)
+ {
+ qn = un = 0.0;
+ }
+ else
+ {
+ qn = 0.5;
+ un =
+ (3.0 / (x[n - 1] - x[n - 2]))
+ * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
+ }
+ y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
+ for (k = n - 1; k >= 1; k--)
+ {
+ y2[k - 1] = y2[k - 1] * y2[k] + u[k - 1];
+ }
+
+ return y2;
+ }
+
+ /**
+ * Calculates the second derivative of the data. This algorithm
+ * was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 113-116.
+ */
+ public double splineInterpolation(
+ double[] xa,
+ double[] ya,
+ double[] y2a,
+ double x)
+ throws InterpolationException
+ {
+ int n = xa.length;
+ if (n != ya.length || n != y2a.length)
+ {
+ throw new InterpolationException("Arrays have different lengths.");
+ }
+ double y;
+ int klo, khi, k;
+ double h, b, a;
+
+ klo = 0;
+ khi = n - 1;
+ while (khi - klo > 1)
+ {
+ k = (khi + klo) >> 1;
+ if (xa[k] > x)
+ khi = k;
+ else
+ klo = k;
+ }
+ h = xa[khi] - xa[klo];
+ //System.out.println(""+x+" between "+xa[khi]+" "+xa[klo]);
+ if (h == 0.0)
+ throw new InterpolationException("Two identical x values. The values must be distinct.");
+ a = (xa[khi] - x) / h;
+ b = (x - xa[klo]) / h;
+ y =
+ a * ya[klo]
+ + b * ya[khi]
+ + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi])
+ * (h * h)
+ / 6.0;
+ return y;
+ }
+
+ /**
+ * Calculates the second derivative of the data. This algorithm
+ * was taken from:
+ *
+ * William H. Press, Saul A. Teukolosky, William T. Vetterling, Brian P. Flannery,
+ * "Numerical Recipes in C - The Art of Scientific Computing", Second Edition,
+ * Cambridge University Press,
+ * ISBN 0-521-43108-5,
+ * chapter 3, pages 113-116.
+ */
+ public double splineInterpolatedDerivative(
+ double[] xa,
+ double[] ya,
+ double[] y2a,
+ double x)
+ throws InterpolationException
+ {
+ int n = xa.length;
+ if (n != ya.length || n != y2a.length)
+ {
+ throw new InterpolationException("Arrays have different lengths.");
+ }
+ double dydx;
+ int klo, khi, k;
+ double h, b, a;
+
+ klo = 0;
+ khi = n - 1;
+ while (khi - klo > 1)
+ {
+ k = (khi + klo) >> 1;
+ if (xa[k] > x)
+ khi = k;
+ else
+ klo = k;
+ }
+ h = xa[khi] - xa[klo];
+ //System.out.println(""+x+" between "+xa[khi]+" "+xa[klo]);
+ if (h == 0.0)
+ throw new InterpolationException("Two identical x values. The values must be distinct.");
+ a = (xa[khi] - x) / h;
+ b = (x - xa[klo]) / h;
+ dydx =
+ (ya[khi] - ya[klo]) / h
+ - ((3 * (a * a) - 1) * h * y2a[klo]) / 6.0
+ + ((3 * (b * b) - 1) * h * y2a[khi]) / 6.0;
+ return dydx;
+ }
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
\ No newline at end of file
diff --git a/src/wsi/ra/print/PagePrinter.java b/src/wsi/ra/print/PagePrinter.java
new file mode 100644
index 00000000..9fd7746c
--- /dev/null
+++ b/src/wsi/ra/print/PagePrinter.java
@@ -0,0 +1,156 @@
+/**
+ * Filename: $RCSfile: PagePrinter.java,v $
+ * Purpose:
+ * Language: Java
+ * Compiler: JDK 1.3
+ * Authors: Fabian Hennecke
+ * Version: $Revision: 1.1.1.1 $
+ * $Date: 2003/07/03 14:59:40 $
+ * $Author: ulmerh $
+ * Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+ */
+
+package wsi.ra.print;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+import java.awt.* ;
+import java.awt.print.* ;
+import java.awt.image.BufferedImage ;
+import javax.swing.* ;
+
+/*==========================================================================*
+ * CLASS DECLARATION
+ *==========================================================================*/
+
+public class PagePrinter
+{
+ Component c;
+ Graphics g;
+ PageFormat pf;
+ public boolean fit_in_possible = true;
+
+ public PagePrinter( Component c, Graphics g, PageFormat pf ) {
+ this.c = c;
+ this.g = g;
+ this.pf = pf;
+ }
+
+ public int print(){
+ Dimension old = c.getSize();
+
+ int x = (int)pf.getImageableX(),
+ y = (int)pf.getImageableY();
+
+ g.translate( x, y );
+
+ double w = (int)pf.getImageableWidth(),
+ h = (int)pf.getImageableHeight();
+
+ if( old.width > w || old.height > h ){
+ boolean rec_turn = false, rec_fit_in = false;
+ if( ( old.width > old.height && h > w ) ||
+ ( old.width < old.height && h < w ) ) {
+ rec_turn = true;
+ if( old.width > h || old.height > w ) rec_fit_in = true;
+ }
+ else rec_fit_in = true;
+
+ JLabel[] text = new JLabel[4];
+ text[0] = new JLabel("The component which should be printed");
+ text[1] = new JLabel("is too large.");
+ text[2] = new JLabel("You can choose if the component should be");
+ JCheckBox cbFitIn = new JCheckBox("fitted-in", rec_fit_in),
+ cbTurn = new JCheckBox("turned", rec_turn );
+ text[3] = new JLabel("(Recommended choice is pre-selected)");
+
+ if( !fit_in_possible ){
+ cbFitIn.setEnabled( false );
+ cbFitIn.setSelected( false );
+ }
+
+ GridBagLayout gbl = new GridBagLayout();
+ JPanel panel = new JPanel( gbl );
+ GridBagConstraints gbc = new GridBagConstraints();
+ gbc.gridx = gbc.gridy = 0;
+ gbc.weightx = gbc.weighty = .5;
+
+ gbc.gridwidth = 2;
+ gbl.setConstraints( text[0], gbc );
+ panel.add( text[0] );
+
+ gbc.gridy++;
+ gbl.setConstraints( text[1], gbc );
+ panel.add( text[1] );
+
+ gbc.gridy++;
+ gbl.setConstraints( text[2], gbc );
+ panel.add( text[2] );
+
+ gbc.gridy++;
+ gbc.gridwidth = 1;
+ gbl.setConstraints( cbFitIn, gbc );
+ panel.add( cbFitIn );
+
+ gbc.gridx++;
+ gbl.setConstraints( cbTurn, gbc );
+ panel.add( cbTurn );
+ gbc.gridx = 0;
+ gbc.gridwidth = 2;
+ gbc.gridy++;
+ gbl.setConstraints( text[3], gbc);
+ panel.add( text[3] );
+
+ int choice = JOptionPane.showOptionDialog( c, panel, "Fit-in",
+ JOptionPane.OK_CANCEL_OPTION,
+ JOptionPane.QUESTION_MESSAGE,
+ null, null, null );
+
+ if( choice == JOptionPane.CANCEL_OPTION || choice == JOptionPane.CLOSED_OPTION )
+ return Printable.NO_SUCH_PAGE;
+
+ else if( choice == JOptionPane.OK_OPTION ){
+
+ if( cbTurn.isSelected() ){
+ BufferedImage img;
+ if( cbFitIn.isSelected() ){
+ double m = Math.min( h / (double)old.width, w / (double)old.height );
+ img = (BufferedImage)c.createImage( (int)( old.height * m ), (int)( old.width * m ) );
+ Graphics2D g2 = img.createGraphics();
+ g2.rotate( Math.toRadians( 90 ) );
+ g2.translate( 0, - old.height * m );
+ c.setSize( (int)( old.width * m ), (int)( old.height * m ) );
+ c.paint( g2 );
+ c.setSize( old );
+ }
+ else{
+ img = (BufferedImage)c.createImage( old.height, old.width );
+ Graphics2D g2 = img.createGraphics();
+ g2.rotate( Math.toRadians( 90 ) );
+ g2.translate( 0, - old.height );
+ c.paint( g2 );
+ }
+ g.drawImage( img, 0, 0, c.getBackground(), c );
+ }
+
+ else{
+ if( cbFitIn.isSelected() ){
+ double m = Math.min( w / (double)old.width, h / (double)old.height );
+ c.setSize( (int)( old.width * m ), (int)( old.height * m ) );
+ c.paint( g );
+ c.setSize( old );
+ }
+ else c.paint( g );
+ }
+ }
+ }
+ else c.paint( g );
+ return Printable.PAGE_EXISTS;
+ }
+}
+
+/****************************************************************************
+ * END OF FILE
+ ****************************************************************************/
diff --git a/src/wsi/ra/sort/XYDoubleArray.java b/src/wsi/ra/sort/XYDoubleArray.java
new file mode 100644
index 00000000..263ae6af
--- /dev/null
+++ b/src/wsi/ra/sort/XYDoubleArray.java
@@ -0,0 +1,114 @@
+/**
+ * Filename: $RCSfile: XYDoubleArray.java,v $
+ * Purpose: A description of the contents of this file.
+ * Language: Java
+ * Compiler: JDK 1.4
+ * Authors: Joerg K. Wegner
+ * Version: $Revision: 1.6 $
+ * $Date: 2001/10/17 11:37:42 $
+ * $Author: wegnerj $
+ * Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+ */
+
+/*==========================================================================*
+ * PACKAGE
+ *==========================================================================*/
+
+package wsi.ra.sort;
+
+/*==========================================================================*
+ * IMPORTS
+ *==========================================================================*/
+
+//import wsi.ra.jchart2d.*;
+import java.util.*;
+
+/*==========================================================================*
+ * CLASS DECLARATION
+ *==========================================================================*/
+
+/**
+ * Defines two double[] arrays.
+ *
+ * @version $Revision: 1.6 $, $Date: 2001/10/17 11:37:42 $
+ *
+ */
+public class XYDoubleArray implements Cloneable{
+ /*-------------------------------------------------------------------------*
+ * public member variables
+ *-------------------------------------------------------------------------*/
+ //public static final int TABLE_OUTPUT=JChart2DUtils.TABLE_OUTPUT;
+ //public static final int VECTOR_OUTPUT=JChart2DUtils.VECTOR_OUTPUT;
+ public double[] x=null;
+ public double[] y=null;
+
+ /*-------------------------------------------------------------------------*
+ * constructor
+ *-------------------------------------------------------------------------*/
+ public XYDoubleArray(int length)
+ {
+ this.x=new double[length];
+ this.y=new double[length];
+ }
+
+ public XYDoubleArray(double[] x, double[] y)
+ {
+ this.x=x;
+ this.y=y;
+ }
+
+ /*-------------------------------------------------------------------------*
+ * public methods
+ *-------------------------------------------------------------------------*/
+
+ public final void swap(int a, int b)
+ {
+ double xx = x[a]; x[a] = x[b]; x[b] = xx;
+ double yy = y[a]; y[a] = y[b]; y[b] = yy;
+ }
+
+ /*public final void sortX(){
+ QuickInsertSort quickInsertSort = new QuickInsertSort();
+ quickInsertSort.sortX(this);
+ }
+
+ public final void sortY(){
+ QuickInsertSort quickInsertSort = new QuickInsertSort();
+ quickInsertSort.sortY(this);
+ }*/
+
+ //
+ // uses jchart2d methods !;-(
+ //
+ /*public String toString(){
+ return toString(TABLE_OUTPUT);
+ }
+
+ private String toString(int output){
+ if(x==null || y==null)return null;
+ BasicDataSet data=new BasicDataSet(x, y,
+ AbstractDataSet.CONTINUOUS_DATA,
+ "x",
+ "y");
+ String s=JChart2DUtils.writeDataSetToString(data, output);
+ return s;
+ }*/
+
+ public Object clone()
+ {
+ double[] newX=new double[x.length];
+ double[] newY=new double[y.length];
+ XYDoubleArray newArray=new XYDoubleArray(newX, newY);
+
+ for(int i=0;i lineData = new ArrayList(100);
+
+ int lineCnt = 0;
+ try
+ {
+ while ((line = lnr.readLine()) != null)
+ {
+ line = line.trim();
+ if (strStartsWithPrefix(line, ignorePrefix) < 0) {
+ if (lineCnt >= lOffset) lineData.add(line);
+ lineCnt++;
+ if ((lCnt > 0) && (lineData.size() == lCnt)) break;
+ }
+// if (line.trim().length() > 0) {
+// if ((ignorePrefix == null) || (strStartsWithPrefix(line, ignorePrefix) < 0)) {
+// // if there are no prefixes given or none of them fits, add the line
+// lineData.add(line);
+// }
+// }
+ }
+ }
+ catch (IOException ex)
+ {
+ logger.error(ex.getMessage());
+ }
+
+ return lineData;
+ }
+
+ /**
+ * Parse columns of a data array containing double data. Columns may be selected by giving their
+ * indices in an int array. If selectedCols is null, all columns are selected. All selected columns
+ * are expected to contain double data and to be of same length. If rawData is null, null is returned.
+ *
+ * @param rawData Strings containing an array with double data columns
+ * @param colSplit String regexp for the splitting of a line
+ * @param selectedCols indices of the columns to retrieve, null for all.
+ * @return
+ */
+ public static double[][] parseDoubleArray(ArrayList rawData, String colSplit, int[] selectedCols) {
+ String[] entries;
+ double dat[][] = null;
+ if (rawData != null) {
+ try {
+ for (int i=0; i= start) && (data[i][col] <= end)) {
+ if (cnt == 0) startIndex = i;
+ cnt++;
+ } else if (cnt > 0) break;
+ }
+
+ double[][] selData = new double[cnt][data[0].length];
+ System.arraycopy(data, startIndex, selData, 0, cnt);
+ return selData;
+ }
+
+ /**
+ * Load double data from a text file. An ignore list of prefixes may be specified. The start line and number of lines
+ * to read may be specified, if lCnt is -1, as many lines as possible are read. The cols array may contain an integer
+ * list of columns to be read. If null, as many columns as possible are read.
+ * The data file is expected to be uniform, meaning that all lines which are not ignored, contain double data values
+ * in all columns.
+ *
+ * @param fname file name to read
+ * @param ignorePrefix lines starting with any of these Strings will be ignored
+ * @param colSplit String regexp for the splitting of a line
+ * @param lOffset start at a certain line (0 for top)
+ * @param lCnt read as many lines, -1 or 0 for all (from offset). Ignored lines are not counted!
+ * @param selectedCols indices of the columns to retrieve, null for all.
+ * @return
+ */
+ public static double[][] loadDoubleData(String fname, String[] ignorePrefix, String colSplit, int lOffset, int lCnt, int[] selectedCols) {
+ return parseDoubleArray((ArrayList)readLines(fname, ignorePrefix, lOffset, lCnt), colSplit, selectedCols);
+ }
+
+ /**
+ * Fill a line of an array with double values parsed from a String array. A subset of
+ * Columns may be selected by giving their indeces in an integer array cols. If cols
+ * is null, all are converted.
+ *
+ * @param dest
+ * @param lineCnt
+ * @param entries
+ * @param cols
+ */
+ public static void fillLine(double[][] dest, int lineCnt, String[] entries, int[] cols) {
+ if (((cols == null) && (dest[lineCnt].length != entries.length)) || (cols != null && (dest[lineCnt].length != cols.length))) {
+ System.err.println("error, array dimensions dont match! (SIFTComparisonMatrix");
+ }
+ if (cols == null) {
+ for (int i=0; i 1) &&
+ (resourceLocation.charAt(1) == ':')))
+ {
+ return getBytesFromFile(resourceLocation);
+ }
+
+ InputStream in = ClassLoader.getSystemResourceAsStream(resourceLocation);
+
+ if (in == null)
+ {
+ // try again for web start applications
+ in = this.getClass().getClassLoader().getResourceAsStream(
+ resourceLocation);
+ }
+
+ if (in == null)
+ {
+ return null;
+ }
+
+ if (logger.isDebugEnabled())
+ {
+ logger.debug("Stream opened for " + resourceLocation);
+ }
+
+ byte[] bytes = getBytesFromStream(in);
+
+ return bytes;
+ }
+
+ /**
+ * Gets the byte data from a file contained in a JAR or ZIP file.
+ *
+ * @param urlToZipArchive Description of the Parameter
+ * @param internalArchivePath Description of the Parameter
+ * @return the byte array of the file.
+ */
+ private byte[] getBytesFromArchive(String urlToZipArchive,
+ String internalArchivePath)
+ {
+ URL url = null;
+ int size = -1;
+ byte[] b = null;
+
+ try
+ {
+ url = new URL(urlToZipArchive);
+
+ // extracts just sizes only.
+ ZipFile zf = new ZipFile(url.getFile());
+ Enumeration e = zf.entries();
+
+ while (e.hasMoreElements())
+ {
+ ZipEntry ze = (ZipEntry) e.nextElement();
+
+ if (ze.getName().equals(internalArchivePath))
+ {
+ if (ze.isDirectory())
+ {
+ return null;
+ }
+
+ // only files with <65536 bytes are allowed
+ if (ze.getSize() > 65536)
+ {
+ System.out.println(
+ "Resource files should be smaller than 65536 bytes...");
+ }
+
+ size = (int) ze.getSize();
+ }
+ }
+
+ zf.close();
+
+ FileInputStream fis = new FileInputStream(url.getFile());
+ BufferedInputStream bis = new BufferedInputStream(fis);
+ ZipInputStream zis = new ZipInputStream(bis);
+ ZipEntry ze = null;
+
+ while ((ze = zis.getNextEntry()) != null)
+ {
+ if (ze.getName().equals(internalArchivePath))
+ {
+ b = new byte[(int) size];
+
+ int rb = 0;
+ int chunk = 0;
+
+ while (((int) size - rb) > 0)
+ {
+ chunk = zis.read(b, rb, (int) size - rb);
+
+ if (chunk == -1)
+ {
+ break;
+ }
+
+ rb += chunk;
+ }
+ }
+ }
+ }
+ catch (Exception e)
+ {
+ logger.error(e.getMessage());
+
+ return null;
+ }
+
+ return b;
+ }
+
+ /**
+ * Gets the byte data from a file.
+ *
+ * @param fileName Description of the Parameter
+ * @return the byte array of the file.
+ */
+ private byte[] getBytesFromFile(String fileName)
+ {
+ if (fileName.startsWith("/cygdrive/"))
+ {
+ int length = "/cygdrive/".length();
+ fileName = fileName.substring(length, length + 1) + ":" +
+ fileName.substring(length + 1);
+ }
+
+ if (logger.isDebugEnabled())
+ {
+ logger.debug("Trying to get file from " + fileName);
+ }
+
+ File file = new File(fileName);
+ FileInputStream fis = null;
+
+ try
+ {
+ fis = new FileInputStream(file);
+ }
+ catch (Exception e)
+ {
+ logger.error(e.getMessage());
+
+ return null;
+ }
+
+ BufferedInputStream bis = new BufferedInputStream(fis);
+
+ // only files with <65536 bytes are allowed
+ //if( file.length() > 65536 ) System.out.println("Resource files should be smaller than 65536 bytes...");
+ int size = (int) file.length();
+ byte[] b = new byte[size];
+ int rb = 0;
+ int chunk = 0;
+
+ try
+ {
+ while (((int) size - rb) > 0)
+ {
+ chunk = bis.read(b, rb, (int) size - rb);
+
+ if (chunk == -1)
+ {
+ break;
+ }
+
+ rb += chunk;
+ }
+ }
+ catch (Exception e)
+ {
+ logger.error(e.getMessage());
+
+ return null;
+ }
+
+ return b;
+ }
+
+ /**
+ * Gets the byte data from a file.
+ *
+ * @param fileName Description of the Parameter
+ * @return the byte array of the file.
+ */
+ private byte[] getBytesFromStream(InputStream stream)
+ {
+ if (stream == null)
+ {
+ return null;
+ }
+
+ if (logger.isDebugEnabled())
+ {
+ logger.debug("Trying to get file from stream.");
+ }
+
+ BufferedInputStream bis = new BufferedInputStream(stream);
+
+ try
+ {
+ int size = (int) bis.available();
+ byte[] b = new byte[size];
+ int rb = 0;
+ int chunk = 0;
+
+ while (((int) size - rb) > 0)
+ {
+ chunk = bis.read(b, rb, (int) size - rb);
+
+ if (chunk == -1)
+ {
+ break;
+ }
+
+ rb += chunk;
+ }
+
+ return b;
+ }
+ catch (Exception e)
+ {
+ logger.error(e.getMessage());
+
+ return null;
+ }
+ }
+
+ /**
+ *
+ */
+ public static Properties readProperties(String resourceName) throws Exception {
+// if (TRACE)
+// System.out.println("EvAClient.readProperties of " + resourceName);
+ Properties prop = new Properties();
+ BasicResourceLoader loader = BasicResourceLoader.instance();
+
+ byte bytes[] = loader.getBytesFromResourceLocation(resourceName);
+ if (bytes != null) {
+ ByteArrayInputStream bais = new ByteArrayInputStream(bytes);
+ prop.load(bais);
+ }
+ if (prop != null)
+ return prop;
+ /////////////
+
+ int slInd = resourceName.lastIndexOf('/');
+ if (slInd != -1)
+ resourceName = resourceName.substring(slInd + 1);
+ Properties userProps = new Properties();
+ File propFile = new File(File.separatorChar + "resources" + File.separatorChar + resourceName);
+ if (propFile.exists()) {
+ try {
+ userProps.load(new FileInputStream(propFile));
+ } catch (Exception ex) {
+ System.out.println("Problem reading user properties: " + propFile);
+ }
+ }
+ return userProps;
+ }
+
+
+}
diff --git a/src/wsi/ra/tool/DummyCategory.java b/src/wsi/ra/tool/DummyCategory.java
new file mode 100644
index 00000000..750cec67
--- /dev/null
+++ b/src/wsi/ra/tool/DummyCategory.java
@@ -0,0 +1,41 @@
+package wsi.ra.tool;
+
+/**
+ * Dummy class replacing the log4j Category because log4j couldnt be included in a clean
+ * way and seemingly wasnt used in the main classes anyways.
+ *
+ * @author mkron, mpaly
+ *
+ */
+public class DummyCategory {
+ static DummyCategory dummy = new DummyCategory();
+ static boolean bDebugEnabled = false;
+
+ public void error() {
+ System.err.println("Error");
+ }
+
+ public void error(String msg) {
+ System.err.println(msg);
+ }
+
+ public static DummyCategory getInstance(String str) {
+ return dummy;
+ }
+
+ public boolean isDebugEnabled() {
+ return bDebugEnabled;
+ }
+
+ public void debug(String str) {
+ System.err.println(str);
+ }
+
+ public void info(String str) {
+ System.err.println(str);
+ }
+
+ public void warn(String str) {
+ System.err.println(str);
+ }
+}
\ No newline at end of file
diff --git a/src/wsi/ra/tool/IntegerArrayList.java b/src/wsi/ra/tool/IntegerArrayList.java
new file mode 100644
index 00000000..f8d78540
--- /dev/null
+++ b/src/wsi/ra/tool/IntegerArrayList.java
@@ -0,0 +1,161 @@
+/**
+ * Filename: $RCSfile: IntegerArrayList.java,v $
+ * Purpose: This class is similar to 'java.util.ArrayList', except that it can
+ * only hold and manage integer values.
+ * Language: Java
+ * Compiler: JDK 1.2
+ * Authors: Fred Rapp
+ * Version: $Revision: 1.1.1.1 $
+ * $Date: 2003/07/03 14:59:40 $
+ * $Author: ulmerh $
+ * Copyright (c) Dept. Computer Architecture, University of Tuebingen, Germany
+ */
+
+package wsi.ra.tool;
+
+/*==========================================================================*
+ * CLASS DECLARATION
+ *==========================================================================*/
+
+/**
+ * This class is similar to 'java.util.ArrayList', except that it can
+ * only hold and manage integer values.
+ */
+public class IntegerArrayList
+{
+ /*-------------------------------------------------------------------------*
+ * private member variables
+ *-------------------------------------------------------------------------*/
+
+ private int[] integerArray = null;
+ private int integerCount = 0;
+
+ private int initialSize;
+ private int incrementFactor;
+
+ /*------------------------------------------------------------------------*
+ * constructor
+ *------------------------------------------------------------------------*/
+
+ public IntegerArrayList()
+ {
+ // call other constructor with default values
+ this(100, 3);
+ }
+
+ public IntegerArrayList(int initialSize, int incrementFactor)
+ {
+ this.initialSize = initialSize;
+ this.incrementFactor = incrementFactor;
+
+ // create new integer array of initial size
+ integerArray = new int[initialSize];
+ }
+
+ public void destroy() { integerArray = null; integerCount = 0; }
+
+ /*-------------------------------------------------------------------------*
+ * public methods
+ *-------------------------------------------------------------------------*/
+
+ /**
+ * Clears the contents of this integer list and sets it back to it's initial state.
+ */
+ public void clear() { integerArray = new int[initialSize]; integerCount = 0; }
+
+ /**
+ * Returns the number of components in this list.
+ */
+ public int size() { return integerCount; }
+
+ /**
+ * Parses given string to integer and adds it to the integer list.
+ * @return true if parsing was successful, false if an error occured.
+ */
+ public boolean add(String stringValue)
+ {
+ boolean success = false;
+ try {
+ int value = Integer.parseInt(stringValue);
+ add(value);
+ success = true;
+ }
+ catch (Exception e) {}
+ return success;
+ }
+
+ /**
+ * Adds a new integer value to the integer list.
+ */
+ public void add(int value)
+ {
+ // check integer array size
+ if (integerCount == integerArray.length)
+ {
+ // increase size of int array
+ int old_length = integerArray.length;
+ int new_length = incrementFactor * old_length;
+ int[] new_array = new int[new_length];
+ for (int i=0; i= integerCount)) { return Integer.MIN_VALUE; }
+
+ return integerArray[index];
+ }
+
+ /**
+ * Returns the contents of this list as an array of integer values.
+ * @param exactLength determines if the returned array should have the exactly right length, or if longer arrays are allowed
+ */
+ public int[] getIntegerArray(boolean exactLength)
+ {
+ // check if we can simply return our internal integer array
+ if (!exactLength || (integerArray.length == integerCount)) { return integerArray; }
+
+ // make a copy of the array with the exactly right length
+ int size = integerCount;
+ int[] new_array = new int[ size];
+ for (int i=0; i< size; i++) { new_array[i] = integerArray[i]; }
+ return new_array;
+ }
+
+ /**
+ * Returns the contents of this list as an array of double values.
+ */
+ public double[] getDoubleArray()
+ {
+ // make a copy of the array
+ int size = integerCount;
+ double[] new_array = new double[ size];
+ for (int i=0; i< size; i++) { new_array[i] = integerArray[i]; }
+ return new_array;
+ }
+
+ /**
+ * Returns the contents of this list as an array of strings.
+ */
+ public String[] getStringArray()
+ {
+ // make a copy of the array
+ int size = integerCount;
+ String[] new_array = new String[ size];
+ for (int i=0; i maximum)) {
+ maxIndex = i;
+ maximum = doubles[i];
+ }
+ }
+
+ return maxIndex;
+ }
+
+ /**
+ * Returns index of maximum element in a given
+ * array of integers. First maximum is returned.
+ *
+ * @param ints the array of integers
+ * @return the index of the maximum element
+ */
+ public static int maxIndex(int [] ints) {
+
+ int maximum = 0;
+ int maxIndex = 0;
+
+ for (int i = 0; i < ints.length; i++) {
+ if ((i == 0) || (ints[i] > maximum)) {
+ maxIndex = i;
+ maximum = ints[i];
+ }
+ }
+
+ return maxIndex;
+ }
+
+ /**
+ * Computes the mean for an array of doubles.
+ *
+ * @param vector the array
+ * @return the mean
+ */
+ public static double mean(double[] vector) {
+
+ double sum = 0;
+
+ if (vector.length == 0) {
+ return 0;
+ }
+ for (int i = 0; i < vector.length; i++) {
+ sum += vector[i];
+ }
+ return sum / (double) vector.length;
+ }
+
+ /**
+ * Returns index of minimum element in a given
+ * array of integers. First minimum is returned.
+ *
+ * @param ints the array of integers
+ * @return the index of the minimum element
+ */
+ public static int minIndex(int [] ints) {
+
+ int minimum = 0;
+ int minIndex = 0;
+
+ for (int i = 0; i < ints.length; i++) {
+ if ((i == 0) || (ints[i] < minimum)) {
+ minIndex = i;
+ minimum = ints[i];
+ }
+ }
+
+ return minIndex;
+ }
+
+ /**
+ * Returns index of minimum element in a given
+ * array of doubles. First minimum is returned.
+ *
+ * @param doubles the array of doubles
+ * @return the index of the minimum element
+ */
+ public static int minIndex(double [] doubles) {
+
+ double minimum = 0;
+ int minIndex = 0;
+
+ for (int i = 0; i < doubles.length; i++) {
+ if ((i == 0) || (doubles[i] < minimum)) {
+ minIndex = i;
+ minimum = doubles[i];
+ }
+ }
+
+ return minIndex;
+ }
+
+ /**
+ * Normalizes the doubles in the array by their sum.
+ *
+ * @param doubles the array of double
+ * @exception IllegalArgumentException if sum is Zero or NaN
+ */
+ public static void normalize(double[] doubles) {
+
+ double sum = 0;
+ for (int i = 0; i < doubles.length; i++) {
+ sum += doubles[i];
+ }
+ normalize(doubles, sum);
+ }
+
+ /**
+ * Normalizes the doubles in the array using the given value.
+ *
+ * @param doubles the array of double
+ * @param sum the value by which the doubles are to be normalized
+ * @exception IllegalArgumentException if sum is zero or NaN
+ */
+ public static void normalize(double[] doubles, double sum) {
+
+ if (Double.isNaN(sum)) {
+ throw new IllegalArgumentException("Can't normalize array. Sum is NaN.");
+ }
+ if (sum == 0) {
+ // Maybe this should just be a return.
+ throw new IllegalArgumentException("Can't normalize array. Sum is zero.");
+ }
+ for (int i = 0; i < doubles.length; i++) {
+ doubles[i] /= sum;
+ }
+ }
+
+
+ /**
+ * Computes the variance for an array of doubles.
+ *
+ * @param vector the array
+ * @return the variance
+ */
+ public static double variance(double[] vector) {
+
+ double sum = 0, sumSquared = 0;
+
+ if (vector.length <= 1) {
+ return 0;
+ }
+ for (int i = 0; i < vector.length; i++) {
+ sum += vector[i];
+ sumSquared += (vector[i] * vector[i]);
+ }
+ double result = (sumSquared - (sum * sum / (double) vector.length)) /
+ (double) (vector.length - 1);
+
+ // We don't like negative variance
+ if (result < 0) {
+ return 0;
+ } else {
+ return result;
+ }
+ }
+
+ /**
+ * Computes the sum of the elements of an array of doubles.
+ *
+ * @param doubles the array of double
+ * @return the sum of the elements
+ */
+ public static double sum(double[] doubles) {
+
+ double sum = 0;
+
+ for (int i = 0; i < doubles.length; i++) {
+ sum += doubles[i];
+ }
+ return sum;
+ }
+
+ /**
+ * Computes the sum of the elements of an array of integers.
+ *
+ * @param ints the array of integers
+ * @return the sum of the elements
+ */
+ public static int sum(int[] ints) {
+
+ int sum = 0;
+
+ for (int i = 0; i < ints.length; i++) {
+ sum += ints[i];
+ }
+ return sum;
+ }
+
+ /**
+ * Returns c*log2(c) for a given integer value c.
+ *
+ * @param c an integer value
+ * @return c*log2(c) (but is careful to return 0 if c is 0)
+ */
+ public static final double xlogx(int c) {
+
+ if (c == 0) {
+ return 0.0;
+ }
+ return c * StatisticUtils.log2((double) c);
+ }
+
+ /**
+ * Returns c*log2(c) for a given value c.
+ *
+ * @param c an integer value
+ * @return c*log2(c) (but is careful to return 0 if c is 0)
+ */
+ public static final double xlogx(double c) {
+
+ if (c == 0) {
+ return 0.0;
+ }
+ return c * StatisticUtils.log2( c);
+ }
+
+ /**
+ * Tests if a is equal to b.
+ *
+ * @param a a double
+ * @param b a double
+ */
+ public static final boolean eq(double a, double b){
+
+ return (a - b < SMALL) && (b - a < SMALL);
+ }
+
+ /**
+ * Tests if a is smaller or equal to b.
+ *
+ * @param a a double
+ * @param b a double
+ */
+ public static final boolean smOrEq(double a,double b) {
+
+ return (a-b < SMALL);
+ }
+
+ /**
+ * Tests if a is greater or equal to b.
+ *
+ * @param a a double
+ * @param b a double
+ */
+ public static final boolean grOrEq(double a,double b) {
+
+ return (b-a < SMALL);
+ }
+
+ /**
+ * Tests if a is smaller than b.
+ *
+ * @param a a double
+ * @param b a double
+ */
+ public static final boolean sm(double a,double b) {
+
+ return (b-a > SMALL);
+ }
+
+ /**
+ * Tests if a is greater than b.
+ *
+ * @param a a double
+ * @param b a double
+ */
+ public static final boolean gr(double a,double b) {
+
+ return (a-b > SMALL);
+ }
+
+ /**
+ * returns root mean square error.
+ */
+ public static final double rmsError(double array1[], double array2[])
+ {
+ if ((array1 == null) || (array2 == null)) { return -1.0; }
+
+ double errorValueRMS = 0;
+ for (int i=0; i
+#include "JMatLink.h"
+#include
+
+
+#include
+#include
+#include "engine.h"
+
+#define V5_COMPAT
+#define BUFSIZE 2560
+
+ mxArray *T = NULL;
+ mxArray *result = NULL;
+ char buffer[BUFSIZE];
+ double *TP;
+ mxArray *arrayP;
+ jdouble scalar;
+ int engOpenMarkerI = 0; // Remember the engine that was opened with
+ // engOpen-command
+
+ #define enginePMax 10 // Maximum number of engines opened at the
+ // same time.
+ Engine *engineP[ enginePMax ]; // Array of pointers to engines
+ // Element 0 is reserved
+
+ jboolean debugB = JNI_FALSE; // No debug messages from start
+
+
+/******************************************************************************/
+int getUnusedEnginePointer()
+{
+ /* Find unsed entry in array of pointers to engine */
+ /* j=0 is reserved!! */
+
+ int j;
+
+ for (j=1; j= enginePMax )
+ )
+ {
+ return; // pointer out of allowed region
+ }
+
+ // free entry in array
+ engineP[p] = 0;
+
+}
+
+
+JNIEXPORT void JNICALL Java_wsi_ra_tool_matlab_JMatLink_setDebugNATIVE
+ (JNIEnv *env, jobject obj, jboolean d)
+{
+ if (d) debugB = JNI_TRUE;
+ else debugB = JNI_FALSE;
+
+}
+
+
+/********************************************************************************/
+
+/********************************************************************************/
+/******************** int engOpenNATIVE( startcmd ) *********************/
+/********************************************************************************/
+
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engOpenNATIVE__Ljava_lang_String_2
+ (JNIEnv *env, jobject obj, jstring startCmdS_JNI)
+{
+ int engineI;
+ const char *openS;
+
+ if (engOpenMarkerI != 0) return engOpenMarkerI; // engOpen already used before
+ // return handle valid connection
+
+ openS = (*env)->GetStringUTFChars(env, startCmdS_JNI, 0);
+
+ /* find unused entry in engine array */
+ engineI = getUnusedEnginePointer();
+
+ if (engineI==0) return 0; // no more pointers available
+
+ if (!(engineP[engineI] = engOpen(openS)) )
+ {
+ if (debugB) fprintf(stderr, "\nCan't start MATLAB engine\n");
+ (*env)->ReleaseStringUTFChars(env, startCmdS_JNI, openS); // free memory
+ delEnginePointer(engineI);
+ return 0;
+ }
+
+ engOpenMarkerI = engineI; // Remember engine that was opened with engOpen()
+
+ (*env)->ReleaseStringUTFChars(env, startCmdS_JNI, openS); // free memory
+
+ return engineI;
+}
+
+/********************************************************************************/
+/**************** int engOpenSingleUseNATIVE( startcmd ) ****************/
+/********************************************************************************/
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engOpenSingleUseNATIVE
+ (JNIEnv *env, jobject obj, jstring startCmdS_JNI)
+{
+ const char *openS = (*env)->GetStringUTFChars(env, startCmdS_JNI, 0);
+ int retStatus = 0;
+ int engineI;
+
+ /* find unused entry in engine array */
+ engineI = getUnusedEnginePointer();
+
+ if (engineI==0) return 0; // no more pointers avaiblable
+
+ if (!(engineP[engineI] = engOpenSingleUse(openS, NULL, &retStatus)))
+ {
+ if (debugB) fprintf(stderr, "\nCan't start MATLAB engine\n");
+ (*env)->ReleaseStringUTFChars(env, startCmdS_JNI, openS); // free memory
+
+ delEnginePointer(engineI);
+ return 0;
+ }
+
+ (*env)->ReleaseStringUTFChars(env, startCmdS_JNI, openS); // free memory
+
+ return engineI;
+}
+
+
+/********************************************************************************/
+/**************** int engCloseNATIVE( int epI ) ****************/
+/********************************************************************************/
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engCloseNATIVE__I
+ (JNIEnv *env, jobject obj, jint engine)
+{
+ int retValI;
+
+ // Check if engine pointer is within allowed region
+ if (( engine < 1 ) || ( engine >= enginePMax ))
+ {
+ return 0; // Pointer is out of allowed region
+ }
+
+
+ if ( engineP[ engine ] != NULL )
+ {
+ retValI = engClose(engineP[ engine ]);
+ delEnginePointer( engine );
+ if (engine == engOpenMarkerI)
+ {
+ // This engine was opened with engOpen() before
+ engOpenMarkerI = 0;
+ }
+ if (debugB) printf("\n engClose \n");
+ }
+ else
+ {
+ return 0;
+ }
+}
+
+
+
+/********************************************************************************/
+/*********** int engEvalStringNATIVE( int epI, String evalS ) *************/
+/********************************************************************************/
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engEvalStringNATIVE__ILjava_lang_String_2
+ (JNIEnv *env, jobject obj, jint engine, jstring evalS_JNI)
+{
+ int retValI = 0;
+ const char *evalS = (*env)->GetStringUTFChars(env, evalS_JNI, 0);
+
+
+ // Check if engine pointer is within allowed region
+ if (( engine < 1 ) || ( engine >= enginePMax ))
+ {
+ return 0; // Pointer is out of allowed region
+ }
+
+ retValI = engEvalString(engineP[ engine ], evalS);
+
+ //printf("evalString %i",OpenB);
+
+ (*env)->ReleaseStringUTFChars(env, evalS_JNI, evalS); // free memory
+
+ return retValI;
+}
+
+
+
+/********************************************************************************/
+// public native void engOutputBufferNATIVE(int epI, int buflenI );
+/********************************************************************************/
+JNIEXPORT jstring JNICALL Java_wsi_ra_tool_matlab_JMatLink_engOutputBufferNATIVE
+ (JNIEnv *env, jobject obj, jint engine, jint buflen)
+{
+ // !!!!! buflen not implemented yet
+
+ // Check if engine pointer is within allowed region
+ if (( engine < 1 ) || ( engine >= enginePMax ))
+ {
+ return 0; // Pointer is out of allowed region
+ }
+
+ engOutputBuffer(engineP[ engine ], buffer, BUFSIZE);
+
+ if (debugB) printf("JMatLink %s", buffer);
+ return (*env)->NewStringUTF(env, buffer);
+}
+
+
+/********************************************************************************/
+// public native engPutArrayNATIVE(int epI, String nameS, double[][] array );
+/********************************************************************************/
+JNIEXPORT void JNICALL Java_wsi_ra_tool_matlab_JMatLink_engPutArrayNATIVE
+ (JNIEnv *env, jobject obj, jint engine, jstring arrayS_JNI, jobjectArray valueDD_JNI)
+{
+
+ int i, j;
+
+ const char *arrayS;
+ int rowCount;
+ jobject colPtr;
+ int colCount;
+ double *tPtrR;
+ jdouble *arrayElements;
+
+ // Check if engine pointer is within allowed region
+ if (( engine < 1 ) || ( engine >= enginePMax ))
+ {
+ return; // Pointer is out of allowed region
+ }
+
+
+ arrayS = (*env)->GetStringUTFChars(env, arrayS_JNI, 0);
+ rowCount = (*env)->GetArrayLength(env, valueDD_JNI);
+ colPtr = (*env)->GetObjectArrayElement(env, valueDD_JNI, 0);
+ colCount = (*env)->GetArrayLength(env, colPtr);
+ //jboolean *isCopy;
+
+ //double *arrayElements = (*env)->GetDoubleArrayElements(env, colPtr, isCopy);
+
+ if (debugB) printf("engPutArray [][] %s %i %i\n", arrayS, rowCount, colCount);
+
+ T = mxCreateDoubleMatrix(rowCount, colCount, mxREAL);
+ mxSetName(T, arrayS);
+ //printf("matrix created and name set\n");
+ tPtrR = mxGetPr(T);
+
+ for (i=0; iGetObjectArrayElement(env, valueDD_JNI, i);
+ //printf("got colPtr %i\n",i);
+ arrayElements = (*env)->GetDoubleArrayElements(env, colPtr, 0);
+ //printf("got array Elements\n");
+ for (j=0; jReleaseStringUTFChars(env, arrayS_JNI, arrayS); // free memory
+}
+
+
+
+/********************************************************************************/
+// public native double engGetScalarNATIVE(int epI, String nameS);
+/********************************************************************************/
+JNIEXPORT jdouble JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetScalarNATIVE
+ (JNIEnv *env, jobject obj, jint engine, jstring scalarS_JNI)
+{
+
+
+ const char *scalarS= (*env)->GetStringUTFChars(env, scalarS_JNI, 0);
+
+ // Check if engine pointer is within allowed region
+ if (( engine < 1 ) || ( engine >= enginePMax ))
+ {
+ return 0.0; // Pointer is out of allowed region
+ }
+
+
+ if (debugB) printf("native engGetScalar %s \n",scalarS);
+
+ //if (OpenB==0) return 0.0; // engine is not open
+
+ arrayP = engGetArray( engineP[ engine ], scalarS);
+
+ if (arrayP == NULL) {
+ printf("Could not get scalar from MATLAB workspace.\n");
+ (*env)->ReleaseStringUTFChars(env, scalarS_JNI, scalarS); // free memory
+ return 0.0;
+ }
+ else {
+ scalar = mxGetScalar(arrayP);
+ }
+
+ mxDestroyArray(arrayP); // free memory
+ (*env)->ReleaseStringUTFChars(env, scalarS_JNI, scalarS); // free memory
+
+ return scalar;
+}
+
+
+/********************************************************************************/
+/********** double[][] engGetArrayNATIVE( int epI, String nameS ) *********/
+/********************************************************************************/
+JNIEXPORT jobjectArray JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetArrayNATIVE
+ (JNIEnv *env, jobject obj, jint engine, jstring arrayS_JNI)
+{
+ // NOTE: in java there are only 1-dimensional (array[]) arrays.
+ // higher dimensional arrays are arrays of arrays.
+ jarray rowA;
+ jclass rowClass;
+
+ jarray columnA;
+ int in,im;
+ int m = 0;
+ int n = 0;
+
+ jdouble *rowElements;
+ double *TP;
+
+ // convert array name to c-string
+ const char *arrayS = (*env)->GetStringUTFChars(env, arrayS_JNI, 0);
+
+ // Check if engine pointer is within allowed region
+ //if (( engine < 1 ) || ( engine >= enginePMax ))
+ //{
+ // return NULL; // Pointer is out of allowed region
+ //}
+
+
+ if (engineP[ engine ] != NULL)
+ {
+
+ if (debugB) printf("native engGetArray %s \n",arrayS);
+
+ arrayP = engGetArray( engineP[ engine ] , arrayS);
+ if (arrayP == NULL)
+ {
+ printf("Could not get array %s from MATLAB workspace.\n", arrayS);
+ (*env)->ReleaseStringUTFChars(env, arrayS_JNI, arrayS); // free memory
+ return NULL;
+ }
+
+ m = mxGetM( arrayP ); // rows
+ n = mxGetN( arrayP ); // columns
+ TP = mxGetPr( arrayP ); // get pointer to values
+
+
+ /* create an array of double and get its class */
+ rowA = (*env)->NewDoubleArray( env, n); // row vector
+ rowClass = (*env)->GetObjectClass( env, rowA); // row class
+
+ /* create an array of object with the rowClass as the
+ the default element */
+ columnA = (*env)->NewObjectArray( env, m, rowClass, NULL); // column vector
+
+
+ for (im=0; imNewDoubleArray( env, n); // row vector
+ rowClass = (*env)->GetObjectClass( env, rowA); // row class
+ rowElements = (*env)->GetDoubleArrayElements(env, rowA, 0); // row elements
+
+ for (in=0; inSetObjectArrayElement(env, columnA, im, rowA);
+
+ (*env)->ReleaseDoubleArrayElements(env, rowA, rowElements, 0);
+ }
+ } // engineP[ ]
+
+ mxDestroyArray(arrayP); // free memory
+ (*env)->ReleaseStringUTFChars(env, arrayS_JNI, arrayS); // free memory
+
+
+ // are the following two line ok? return NULL ?!?!?
+ if (engineP[ engine ] != NULL) return columnA;
+ else return NULL;
+}
+
+/********************************************************************************/
+// public String[] engGetCharArrayNATIVE(String name)
+/********************************************************************************/
+JNIEXPORT jobjectArray JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetCharArrayNATIVE
+ (JNIEnv *env, jobject jthis, jint engine, jstring arrayS_JNI)
+{
+
+
+ // Check if engine pointer is within allowed region
+ //if (( engine < 1 ) || ( engine >= enginePMax ))
+ //{
+ // return NULL; // Pointer is out of allowed region
+ //}
+
+ const char *arrayS = (*env)->GetStringUTFChars(env, arrayS_JNI, 0);
+
+
+ if (debugB) printf("native engGetArray %s \n",arrayS);
+
+ arrayP = engGetArray( engineP[ engine ], arrayS);
+ if (arrayP == NULL) {
+ printf("Could not get array %s from MATLAB workspace.\n", arrayS);
+ (*env)->ReleaseStringUTFChars(env, arrayS_JNI, arrayS); // free memory
+ return NULL;
+ }
+ if (mxIsChar(arrayP) != 0) {
+ printf("The array %s is not of type char.\n", arrayS);
+ (*env)->ReleaseStringUTFChars(env, arrayS_JNI, arrayS); // free memory
+ return NULL;
+ }
+
+
+ mxDestroyArray(arrayP); // free memory
+ (*env)->ReleaseStringUTFChars(env, arrayS_JNI, arrayS); // free memory
+} // end engGetCharArrayNATIVE
+
+
+
diff --git a/src/wsi/ra/tool/matlab/JMatLink.h b/src/wsi/ra/tool/matlab/JMatLink.h
new file mode 100644
index 00000000..7f61a6bf
--- /dev/null
+++ b/src/wsi/ra/tool/matlab/JMatLink.h
@@ -0,0 +1,149 @@
+/* DO NOT EDIT THIS FILE - it is machine generated */
+#include
+/* Header for class wsi_ra_tool_matlab_JMatLink */
+
+#ifndef _Included_wsi_ra_tool_matlab_JMatLink
+#define _Included_wsi_ra_tool_matlab_JMatLink
+#ifdef __cplusplus
+extern "C" {
+#endif
+/* Inaccessible static: threadInitNumber */
+/* Inaccessible static: stopThreadPermission */
+#undef wsi_ra_tool_matlab_JMatLink_MIN_PRIORITY
+#define wsi_ra_tool_matlab_JMatLink_MIN_PRIORITY 1L
+#undef wsi_ra_tool_matlab_JMatLink_NORM_PRIORITY
+#define wsi_ra_tool_matlab_JMatLink_NORM_PRIORITY 5L
+#undef wsi_ra_tool_matlab_JMatLink_MAX_PRIORITY
+#define wsi_ra_tool_matlab_JMatLink_MAX_PRIORITY 10L
+#undef wsi_ra_tool_matlab_JMatLink_idleI
+#define wsi_ra_tool_matlab_JMatLink_idleI 0L
+#undef wsi_ra_tool_matlab_JMatLink_engOpenI
+#define wsi_ra_tool_matlab_JMatLink_engOpenI 1L
+#undef wsi_ra_tool_matlab_JMatLink_engCloseI
+#define wsi_ra_tool_matlab_JMatLink_engCloseI 2L
+#undef wsi_ra_tool_matlab_JMatLink_engEvalStringI
+#define wsi_ra_tool_matlab_JMatLink_engEvalStringI 3L
+#undef wsi_ra_tool_matlab_JMatLink_engGetScalarI
+#define wsi_ra_tool_matlab_JMatLink_engGetScalarI 4L
+#undef wsi_ra_tool_matlab_JMatLink_engGetVectorI
+#define wsi_ra_tool_matlab_JMatLink_engGetVectorI 5L
+#undef wsi_ra_tool_matlab_JMatLink_engGetArrayI
+#define wsi_ra_tool_matlab_JMatLink_engGetArrayI 6L
+#undef wsi_ra_tool_matlab_JMatLink_engPutArray2dI
+#define wsi_ra_tool_matlab_JMatLink_engPutArray2dI 9L
+#undef wsi_ra_tool_matlab_JMatLink_engOutputBufferI
+#define wsi_ra_tool_matlab_JMatLink_engOutputBufferI 10L
+#undef wsi_ra_tool_matlab_JMatLink_engGetCharArrayI
+#define wsi_ra_tool_matlab_JMatLink_engGetCharArrayI 11L
+#undef wsi_ra_tool_matlab_JMatLink_destroyJMatLinkI
+#define wsi_ra_tool_matlab_JMatLink_destroyJMatLinkI 12L
+#undef wsi_ra_tool_matlab_JMatLink_engOpenSingleUseI
+#define wsi_ra_tool_matlab_JMatLink_engOpenSingleUseI 13L
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: displayHelloWorld
+ * Signature: ()V
+ */
+JNIEXPORT void JNICALL Java_wsi_ra_tool_matlab_JMatLink_displayHelloWorld
+ (JNIEnv *, jobject);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engTestNATIVE
+ * Signature: ()V
+ */
+JNIEXPORT void JNICALL Java_wsi_ra_tool_matlab_JMatLink_engTestNATIVE
+ (JNIEnv *, jobject);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engOpenNATIVE
+ * Signature: (Ljava/lang/String;)I
+ */
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engOpenNATIVE
+ (JNIEnv *, jobject, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engOpenSingleUseNATIVE
+ * Signature: (Ljava/lang/String;)I
+ */
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engOpenSingleUseNATIVE
+ (JNIEnv *, jobject, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engCloseNATIVE
+ * Signature: (I)I
+ */
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engCloseNATIVE
+ (JNIEnv *, jobject, jint);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engEvalStringNATIVE
+ * Signature: (ILjava/lang/String;)I
+ */
+JNIEXPORT jint JNICALL Java_wsi_ra_tool_matlab_JMatLink_engEvalStringNATIVE
+ (JNIEnv *, jobject, jint, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engGetScalarNATIVE
+ * Signature: (ILjava/lang/String;)D
+ */
+JNIEXPORT jdouble JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetScalarNATIVE
+ (JNIEnv *, jobject, jint, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engGetVectorNATIVE
+ * Signature: (ILjava/lang/String;)[D
+ */
+JNIEXPORT jdoubleArray JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetVectorNATIVE
+ (JNIEnv *, jobject, jint, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engGetArrayNATIVE
+ * Signature: (ILjava/lang/String;)[[D
+ */
+JNIEXPORT jobjectArray JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetArrayNATIVE
+ (JNIEnv *, jobject, jint, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engGetCharArrayNATIVE
+ * Signature: (ILjava/lang/String;)[Ljava/lang/String;
+ */
+JNIEXPORT jobjectArray JNICALL Java_wsi_ra_tool_matlab_JMatLink_engGetCharArrayNATIVE
+ (JNIEnv *, jobject, jint, jstring);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engPutArrayNATIVE
+ * Signature: (ILjava/lang/String;[[D)V
+ */
+JNIEXPORT void JNICALL Java_wsi_ra_tool_matlab_JMatLink_engPutArrayNATIVE
+ (JNIEnv *, jobject, jint, jstring, jobjectArray);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: engOutputBufferNATIVE
+ * Signature: (II)Ljava/lang/String;
+ */
+JNIEXPORT jstring JNICALL Java_wsi_ra_tool_matlab_JMatLink_engOutputBufferNATIVE
+ (JNIEnv *, jobject, jint, jint);
+
+/*
+ * Class: wsi_ra_tool_matlab_JMatLink
+ * Method: setDebugNATIVE
+ * Signature: (Z)V
+ */
+JNIEXPORT void JNICALL Java_wsi_ra_tool_matlab_JMatLink_setDebugNATIVE
+ (JNIEnv *, jobject, jboolean);
+
+#ifdef __cplusplus
+}
+#endif
+#endif
diff --git a/src/wsi/ra/tool/matlab/JMatLink.java b/src/wsi/ra/tool/matlab/JMatLink.java
new file mode 100644
index 00000000..d9d48988
--- /dev/null
+++ b/src/wsi/ra/tool/matlab/JMatLink.java
@@ -0,0 +1,2014 @@
+/*****************************************************************************
+
+* JMatLink *
+
+******************************************************************************
+
+* (c) 1999-2001 Stefan Mueller (email: stefan@held-mueller.de) *
+
+******************************************************************************
+
+* 19.01.1999 beginning (reuse of HelloWorld example, Java Tutorial) (0.01) *
+
+* 06.02.1999 separation from gui (0.02) *
+
+* 27.03.1999 engGetScalar, engGetVector, engGetArray (0.03) *
+
+* 02.05.1999 engGetScalar already works *
+
+* 20.05.1999 deleted engPutArray (float..., int...) (0.04) *
+
+* 13.08.1999 thread implementation due to ActiveX thread problem (0.05) *
+
+* 15.08.1999 intruduced 2 locking mechanisms to lock the engine and *
+
+* wait for data to be return to caller. (0.06) *
+
+* 16.08.1999 wait for all engine calls to return savely (0.06a)*
+
+* 29.08.1999 made all engine routines "synchronized" (0.07) *
+
+* 11.09.1999 get char array from workspace (0.08) *
+
+* 21.09.1999 engGetCharArray working (conversion to double matrix first) *
+
+* (transfer as double throught engine and conversion *
+
+* back to double) engine doesn't support char arrays (0.09) *
+
+* 10.10.1999 engPutArray(String name, double[][] values) (0.10) *
+
+* 04/07/2001 restart work on JMatLink (0.11) *
+
+* 05/24/2001 engOpenSingleUse (0.11) *
+
+* 05/24/2001 engClose(int pointer) (0.11) *
+
+* 06/04/2001 engEvalString(ep, evalString) *
+
+* 06/10/2001 use of maximum number of arguments of all engine functions *
+
+* 07/31/2001 setDebug() to suppress messages to the console window *
+
+* 08/09/2001 extensive documentation *
+
+* 08/09/2001 final version. upgrade version number (1.00) *
+
+******************************************************************************/
+
+
+
+//****************************************************************************
+
+// The matlab-engine must only be called from the SAME thread *
+
+// all the time!! On Windows systems the ActiveX implementation is *
+
+// supposed to be the reason for this. I don't know if that happens *
+
+// on other platforms, too. *
+
+// To achieve this. All native methods of this class, all accesses to *
+
+// the engine, are called from the SAME thread. Since I don't know *
+
+// how to call them directly I set up a mechanism to send messages *
+
+// to that thread and ask it to process all requests. Some locking *
+
+// mechanism locks up the engine for each single call in order *
+
+// to stay out of concurrent situations / accesses. *
+
+//****************************************************************************
+
+
+
+
+
+/** ToDo list
+
+* engOpen is supposed to throw an exception if not successfull
+
+* I'm not sure what happens if large matrices are passed back and
+
+ force between this class and other applications. Maybe I
+
+ should always clear all global arrays engGetArray1dD...
+
+* I'm not sure if the locking mechanism is good or bad. There
+
+ may be a problem, if two subsequent calls to two different
+
+ engine routines occur at the same time.
+
+* what should I do with the return values of the engine functions? Throw
+
+ an exception?
+
+* make engPutArray also for int,
+
+* make engGetArray/Scalar also for int
+
+* make engGetCharArray also for different engine pointers
+
+* make something like engGetOutputBuffer
+
+*/
+
+
+
+package wsi.ra.tool.matlab;
+
+
+
+// **** Thread programming. Notify and start and stop ****
+
+// Look at
+
+// file:///e%7C/jdk1.2/docs/guide/misc/threadPrimitiveDeprecation.html
+
+// on how to suspend, start, stop, interrupt a thread savely in JDK1.2.x
+
+// The problem is ActiveX: one only can make calls to the engine library
+
+// (especially engEvalString) from ONE single thread.
+import java.io.*;
+import wsi.ra.tool.*;
+public class JMatLink extends Thread {
+ private static JMatLink m_Instance = null;
+ /**
+ *
+ */
+ public static JMatLink getInstance () {
+ //System.out.println("JMatLink getInstance () {");
+ if (m_Instance==null)
+ m_Instance = new JMatLink();
+ return m_Instance;
+ }
+ // static declarations
+ // the variable "status" is used to tell the main
+ // thread what to do.
+ private final static int idleI = 0;
+ private final static int engOpenI = 1;
+ private final static int engCloseI = 2;
+ private final static int engEvalStringI = 3;
+ private final static int engGetScalarI = 4;
+ private final static int engGetVectorI = 5;
+ private final static int engGetArrayI = 6;
+ private final static int engPutArray2dI = 9;
+ private final static int engOutputBufferI = 10;
+ private final static int engGetCharArrayI = 11;
+ private final static int destroyJMatLinkI = 12;
+ private final static int engOpenSingleUseI = 13;
+ // All variables are global to allow all methods
+
+ // and the main thread to share all data
+
+ private int status = idleI;
+
+ private String arrayS;
+
+ private String engEvalStringS;
+
+ private String engOutputBufferS;
+
+ private double engGetScalarD;
+
+ private double[] engGetVectorD;
+
+ private double[][] engGetArrayD;
+
+ private double engPutArrayD;
+
+ private double[] engPutArray1dD;
+
+ private double[][] engPutArray2dD;
+
+ private String[] engGetCharArrayS;
+
+ private int epI; /* Engine pointer */
+
+ private int retValI; /* return Value of eng-methods */
+
+ private String startCmdS; /* start command for engOpen... */
+
+ private int buflenI; /* output buffer length */
+
+ private boolean debugB = false;
+
+
+
+ // Locks
+
+ private boolean lockEngineB = false;
+
+ private boolean lockThreadB = false;
+
+ private boolean lockWaitForValueB = false;
+
+ private boolean destroyJMatLinkB = false;
+
+
+
+ private Thread runner;
+
+
+
+ // *********************** native declarations ****************************
+
+ // NEVER call native methods directly, like
+
+ // JMatLink.engEvalStringNATIVE("a=1"). Matlab's engine has quite some
+
+ // thread problems.
+
+ private native void displayHelloWorld();
+
+ private native void engTestNATIVE();
+
+
+
+ private native int engOpenNATIVE (String startCmdS );
+
+ private native int engOpenSingleUseNATIVE (String startCmdS );
+
+
+
+ private native int engCloseNATIVE (int epI);
+
+
+
+ private native int engEvalStringNATIVE (int epI, String evalS );
+
+
+
+ private native double engGetScalarNATIVE (int epI, String nameS );
+
+ private native double[] engGetVectorNATIVE (int epI, String nameS );
+
+ private native double[][] engGetArrayNATIVE (int epI, String nameS );
+
+ private native String[] engGetCharArrayNATIVE (int epI, String nameS );
+
+
+
+ private native void engPutArrayNATIVE (int epI, String matrixS, double[][] valuesDD);
+
+
+
+ private native String engOutputBufferNATIVE (int epI, int buflenI );
+
+
+
+ private native void setDebugNATIVE (boolean debugB );
+
+
+
+ // *******************************************************************************
+
+ // ************** load JMatLink library into memory **********************
+
+ static {
+ //System.out.println("loading !!");
+ try { //System.out.println("loading");
+
+ System.loadLibrary("JMatLink");
+
+
+// String path="../lib/"+SystemHelper.getOperationSystemName()+"/JMatLink.dll";
+
+// System.out.println("load: "+path);
+
+// System.loadLibrary(path);
+
+
+
+// String libPath = "." + System.getProperty("file.separator", "/");
+
+// //libPath += SystemHelper.getOperationSystemName() + System.getProperty("file.separator", "/");
+
+// libPath += System.mapLibraryName("JMatLink");
+
+//
+
+// // make sure that we have the absolute path
+
+// libPath = new File(libPath).getAbsolutePath();
+
+//
+
+// System.out.println("load: "+"d:/workingAt/JCompChem/lib/windows/JMatLink.dll");
+
+// System.loadLibrary("../lib/windows/JMatLink_old.dll");
+
+ //System.out.println("loaded");
+
+ }
+
+ catch (UnsatisfiedLinkError e) {
+
+ System.err.println("ERROR: Could not load the JMatLink library");
+
+ e.printStackTrace();
+
+ }
+ //System.out.println("loading !! end");
+
+ }
+
+
+
+
+
+ // ************************ JMatLink constructor ***************************
+
+ /** This is the constructor for the JMatLink library.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ private JMatLink() {
+
+ if (debugB) System.out.println("JMatLink constructor");
+
+ runner = new Thread(this);
+
+ runner.start();
+
+ }
+
+
+
+ // **** terminate running thread ****
+
+ public void destroy() {
+
+ destroyJMatLinkB = true;
+
+ notifyAll();
+
+ }
+
+
+
+ public void kill() {
+
+ destroyJMatLinkB = true;
+
+ callThread(destroyJMatLinkI);
+
+ }
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+ // ***************************** engOpen *****************************
+
+ /** Open engine. This command is used to open a single connection
+
+ * to matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized int engOpen()
+
+ {
+ return engOpen( "" ); // return value is pointer to engine
+
+ }
+
+
+
+
+
+ // ***************************** engOpen *******************************
+
+ /** Open engine. This command is used to open a single connection
+
+ * to matlab. This command is only useful on unix systems. On windows
+
+ * the optional parameter must be NULL.
+
+ *
+
+ *
E.g.:
+
+ *
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen("commands to start matlab");
+
+ * engine.engEvalString("surf(peaks)");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized int engOpen(String startCmdS)
+
+ {
+
+ // startup MATLAB and set up connection
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.startCmdS = startCmdS;
+
+
+
+ callThread( engOpenI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ return this.epI;
+
+ }
+
+
+
+
+
+ // ************************** engOpenSingleUse *****************************
+
+ /** Open engine for single use. This command is used to open
+
+ * multiple connections to matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int a,b;
+
+ * JMatLink engine = new JMatLink();
+
+ * a = engine.engOpenSingleUse(); // start first matlab session
+
+ * b = engine.engOpenSingleUse(); // start second matlab session
+
+ * engine.engEvalString(a, "surf(peaks)");
+
+ * engine.engEvalString(b, "foo=ones(10,0)");
+
+ * engine.engClose(a);
+
+ * engine.engClose(b);
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized int engOpenSingleUse()
+
+ {
+
+ return engOpenSingleUse("");
+
+ }
+
+
+
+ // ************************** engOpenSingleUse *****************************
+
+ /** Open engine for single use. This command is used to open
+
+ * multiple connections to matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int a,b;
+
+ * JMatLink engine = new JMatLink();
+
+ * a = engine.engOpenSingleUse("start matlab"); // start first matlab session
+
+ * b = engine.engOpenSingleUse("start matlab"); // start second matlab session
+
+ * engine.engEvalString(a, "surf(peaks)");
+
+ * engine.engEvalString(b, "foo=ones(10,0)");
+
+ * engine.engClose(a);
+
+ * engine.engClose(b);
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized int engOpenSingleUse(String startCmdS)
+
+ {
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.startCmdS = startCmdS;
+
+
+
+ callThread( engOpenSingleUseI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ return this.epI;
+
+ }
+
+
+
+
+
+ // ***************************** engClose *****************************
+
+ /** Close the connection to matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engClose()
+
+ {
+
+ engClose( this.epI );
+
+ }
+
+
+
+
+
+ // ***************************** engClose *****************************
+
+ /** Close a specified connection to an instance of matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int a,b;
+
+ * JMatLink engine = new JMatLink();
+
+ * a = engine.engOpenSingleUse(); // start first matlab session
+
+ * b = engine.engOpenSingleUse(); // start second matlab session
+
+ * engine.engEvalString(b, "surf(peaks)");
+
+ * engine.engEvalString(a, "array = randn(23)");
+
+ * engine.engClose(a); // Close the first connection to matlab
+
+ * engine.engClose(b); // Close the second connection to matlab
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engClose( int epI)
+
+ {
+
+ // close connection and terminate MATLAB
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.epI = epI;
+
+
+
+ callThread( engCloseI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ // return retValI; Return value indicates success
+
+ }
+
+
+
+
+
+ // ***************************** engEvalString *****************************
+
+ /** Evaluate an expression in matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engEvalString(String evalS)
+
+ {
+
+ engEvalString( this.epI, evalS);
+
+ }
+
+
+
+
+
+ // ***************************** engEvalString *************************
+
+ /** Evaluate an expression in a specified workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int a,b;
+
+ * JMatLink engine = new JMatLink();
+
+ * a = engine.engOpenSingleUse();
+
+ * engine.engEvalString(a, "surf(peaks)");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engEvalString(int epI, String evalS)
+
+ {
+
+ // evaluate expression "evalS" in specified engine Ep
+
+ if (debugB) System.out.println("eval(ep,String) in " + epI + " "+evalS);
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.epI = epI;
+
+ engEvalStringS = evalS;
+
+
+
+ callThread( engEvalStringI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ if (debugB) System.out.println("eval(ep,String) out "+epI+" "+evalS);
+
+ // return retValI; Return value indicates success
+
+ }
+
+
+
+
+
+ // ***************************** engGetScalar **************************
+
+ /** Get a scalar value from matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * double a;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("foo = sin( 3 )");
+
+ * a = engine.engGetScalarValue("foo");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized double engGetScalar(String arrayS)
+
+ {
+
+ return engGetScalar( this.epI, arrayS);
+
+ }
+
+
+
+
+
+ // ***************************** engGetScalar **************************
+
+ /** Get a scalar value from a specified workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * double a;
+
+ * int b;
+
+ * JMatLink engine = new JMatLink();
+
+ * b = engine.engOpenSigleUse();
+
+ * engine.engEvalString(b, "foo = sin( 3 )");
+
+ * a = engine.engGetScalarValue(b, "foo");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized double engGetScalar(int epI, String arrayS)
+
+ {
+
+ // Get scalar value or element (1,1) of an array from
+
+ // MATLAB's workspace
+
+ // Only real values are supported right now
+
+
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ /* copy parameters to global variables */
+
+ this.epI = epI;
+
+ this.arrayS = arrayS;
+
+
+
+ callThread( engGetScalarI );
+
+
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ return engGetScalarD;
+
+ }
+
+
+
+
+
+ // ***************************** engGetVector **************************
+
+ /** Get an array (1 * n) from matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * double[] array;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("array = randn(10,1);");
+
+ * array = engine.engGetVector("array");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized double[] engGetVector( String arrayS )
+
+ {
+
+ return engGetVector( this.epI, arrayS );
+
+ }
+
+
+
+ // ***************************** engGetVector **************************
+
+ /** Get an array (1 * n) from a specified workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int b;
+
+ * double[] array;
+
+ * JMatLink engine = new JMatLink();
+
+ * b = engine.engOpenSingleUse();
+
+ * engine.engEvalString(b, "array = randn(10,1);");
+
+ * array = engine.engGetVector(b, "array");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized double[] engGetVector( int epI, String arrayS )
+
+ {
+
+ // only real values are supported so far
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.epI = epI;
+
+ this.arrayS = arrayS;
+
+
+
+ callThread( engGetVectorI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ return engGetVectorD;
+
+ }
+
+
+
+
+
+ // ***************************** engGetArray ***************************
+
+ /** Get an array from matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int b;
+
+ * double[][] array;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("array = randn(10);");
+
+ * array = engine.engGetArray("array");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized double[][] engGetArray( String arrayS )
+
+ {
+
+ return engGetArray( this.epI, arrayS );
+
+ }
+
+
+
+
+
+ // ***************************** engGetArray ***************************
+
+ /** Get an array from a specified instance/workspace of matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int b;
+
+ * double[][] array;
+
+ * JMatLink engine = new JMatLink();
+
+ * b = engine.engOpenSingleUse();
+
+ * engine.engEvalString(b, "array = randn(10);");
+
+ * array = engine.engGetArray(b, "array");
+
+ * engine.engClose(b);
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized double[][] engGetArray( int epI, String arrayS )
+
+ {
+
+ // only real values are supported so far
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.epI = epI;
+
+ this.arrayS = arrayS;
+
+
+
+ callThread( engGetArrayI );
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+
+
+ return engGetArrayD;
+
+ }
+
+
+
+ // ************************** engGetCharArray *****************************
+
+ /** Get an 'char' array (string) from matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * String array;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("array = 'hello world';");
+
+ * array = engine.engCharArray("array");
+
+ * System.out.println("output = "+ array);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized String[] engGetCharArray(String arrayS)
+
+ {
+
+ // convert to double array
+
+ engEvalString( "engGetCharArrayD=double(" + arrayS +")" );
+
+
+
+ // get double array
+
+ double[][] arrayD = engGetArray("engGetCharArrayD");
+
+
+
+ // delete temporary double array
+
+ engEvalString("clear engGetCharArrayD");
+
+
+
+ // convert double back to char
+
+ return double2String( arrayD );
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray ***************************
+
+ /** Put an array into a specified workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int array = 1;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engPutArray("array", array);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray( String arrayS, int valueI )
+
+ {
+
+ engPutArray( this.epI, arrayS, new Integer(valueI).doubleValue());
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray ***************************
+
+ // public synchronized void engPutArray( String arrayS, int[] valuesI )
+
+ // {
+
+ // engPutArray( this.epI, arrayS, (double[])valuesI );
+
+ // }
+
+
+
+
+
+ // ***************************** engPutArray ***************************
+
+ /** Put an array into matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * double array = 1;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engPutArray("array", array);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray( String arrayS, double valueD )
+
+ {
+
+ engPutArray( this.epI, arrayS, valueD);
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray *****************************
+
+ /** Put an array into a specified instance/workspace of matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int b;
+
+ * double array = 1;
+
+ * JMatLink engine = new JMatLink();
+
+ * b = engine.engOpenSingleUse();
+
+ * engine.engPutArray(b, "array", array);
+
+ * engine.engClose(b);
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray( int epI, String arrayS, double valueD )
+
+ {
+
+
+
+ double vDD[][] = {{0.0}};
+
+ vDD[0][0] = valueD;
+
+ engPutArray( epI, arrayS, vDD ); // nxn dimensional
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray ***************************
+
+ /** Put an array (1 dimensional) into a specified instance/workspace of
+
+ * matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * double[] array = {1.0 , 2.0 , 3.0};
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engPutArray("array", array);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray( String arrayS, double[] valuesD )
+
+ {
+
+ engPutArray( this.epI, arrayS, valuesD );
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray *****************************
+
+ /** Put an array (1 dimensional) into a specified instance/workspace of
+
+ * matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int b;
+
+ * double[] array = {1.0 , 2.0 , 3.0};
+
+ * JMatLink engine = new JMatLink();
+
+ * b = engine.engOpenSingleUse();
+
+ * engine.engPutArray(b, "array", array);
+
+ * engine.engClose(b);
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray(int epI, String arrayS, double[] valuesD)
+
+ {
+
+ double[][] vDD = new double[1][valuesD.length]; // 1xn array
+
+
+
+ if (debugB) System.out.println("length = "+valuesD.length);
+
+
+
+ vDD[0] = valuesD; // copy row
+
+
+
+ engPutArray( epI, arrayS, vDD );
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray ***************************
+
+ /** Put an array (2 dimensional) into matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * double[][] array={{1.0 , 2.0 , 3.0},
+
+ * {4.0 , 5.0 , 6.0}};
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpenSingleUse();
+
+ * engine.engPutArray("array", array);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray( String arrayS, double[][] valuesDD )
+
+ {
+
+ engPutArray( this.epI, arrayS, valuesDD );
+
+ }
+
+
+
+
+
+ // ***************************** engPutArray ***************************
+
+ /** Put an array (2 dimensional) into a specified instance/workspace of
+
+ * matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * int b;
+
+ * double[][] array={{1.0 , 2.0 , 3.0},
+
+ * {4.0 , 5.0 , 6.0}};
+
+ * JMatLink engine = new JMatLink();
+
+ * b = engine.engOpenSingleUse();
+
+ * engine.engPutArray(b, "array", array);
+
+ * engine.engClose(b);
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized void engPutArray(int epI, String arrayS, double[][] valuesDD)
+
+ {
+
+ // send an array to MATLAB
+
+ // only real values are supported so far
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.epI = epI;
+
+ this.arrayS = arrayS;
+
+ this.engPutArray2dD = valuesDD;
+
+
+
+ callThread( engPutArray2dI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+ }
+
+
+
+
+
+ // ***************************** engOutputBuffer ***********************
+
+ /** Return the outputs of previous commands from matlab's workspace.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * String buffer;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * buffer = engine.engOutputBuffer();
+
+ * System.out.println("workspace " + buffer);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized String engOutputBuffer( )
+
+ {
+
+ return engOutputBuffer( this.epI, this.buflenI );
+
+ }
+
+
+
+ // ***************************** engOutputBuffer ***********************
+
+ /** Return the outputs of previous commands from a specified instance/
+
+ * workspace form matlab.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * String buffer;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * buffer = engine.engOutputBuffer();
+
+ * System.out.println("workspace " + buffer);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized String engOutputBuffer( int epI )
+
+ {
+
+ return engOutputBuffer( epI, this.buflenI );
+
+ }
+
+
+
+ // ***************************** engOutputBuffer ***********************
+
+ /** Return the ouputs of previous commands in matlab's workspace.
+
+ *
+
+ * Right now the parameter buflen is not supported.
+
+ *
+
+ * E.g.:
+
+ *
+
+ * String buffer;
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpen();
+
+ * engine.engEvalString("surf(peaks)");
+
+ * buffer = engine.engOutputBuffer();
+
+ * System.out.println("workspace " + buffer);
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public synchronized String engOutputBuffer( int epI, int buflenI )
+
+ {
+
+ // get the output buffer from MATLAB
+
+ if (debugB) System.out.println("Thread in: "+Thread.currentThread().getName());
+
+
+
+ lockEngineLock();
+
+ lockWaitForValue();
+
+
+
+ this.epI = epI;
+
+ this.buflenI = buflenI;
+
+
+
+ callThread( engOutputBufferI );
+
+
+
+ WaitForValue();
+
+ releaseEngineLock();
+
+ if (debugB) System.out.println("Thread out: "+Thread.currentThread().getName());
+
+
+
+ return engOutputBufferS;
+
+ }
+
+
+
+
+
+ // ***************************** setDebug *******************************
+
+ /* Switch on or disable debug information printed to standard output.
+
+ *
+
+ * Default setting is debug info disabled.
+
+ *
E.g.:
+
+ *
+
+ * JMatLink engine = new JMatLink();
+
+ * engine.engOpenSingleUse();
+
+ * engine.setDebug(true);
+
+ * engine.engEvalString("a=ones(10,5);");
+
+ * engine.engClose();
+
+ *
+
+ ***************************************************************************/
+
+ public void setDebug( boolean debugB )
+
+ {
+
+ this.debugB = debugB;
+
+ setDebugNATIVE( debugB );
+
+ }
+
+
+
+
+
+ ////////////////////////////////////////////////////////////////////////////////
+
+ // This method notifys the main thread to call matlab's engine
+
+ // Since threads don't have methods, we set a variable which
+
+ // contains the necessary information about what to do.
+
+ private synchronized void callThread(int status)
+
+ {
+
+ this.status = status;
+
+ lockThreadB = false;
+
+ notifyAll();
+
+ }
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+// The run methods does ALL calls to the native methods
+
+// The keyword "synchronized" is neccessary to block the run()
+
+// method as long as one command needs to get executed.
+
+public synchronized void run()
+
+{
+
+ int tempRetVal;
+
+
+
+ if (debugB) System.out.println("JMatLink: thread is running");
+
+ while (true) {
+
+ // System.out.println("Number of Java-Threads: "+Thread.activeCount()+"");
+
+ // Thread thread = Thread.currentThread();
+
+ // System.out.println("Name of active Java-Threads: "+thread.getName()+"");
+
+ // System.out.println("active Java-Thread is Daemon: "+thread.isDaemon();+"");
+
+
+
+
+
+ switch (status) {
+
+ case engOpenI: epI = engOpenNATIVE( startCmdS );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engOpenSingleUseI: epI = engOpenSingleUseNATIVE( startCmdS );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engCloseI: retValI = engCloseNATIVE( epI );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engEvalStringI: retValI = engEvalStringNATIVE(epI, engEvalStringS);
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engGetScalarI: engGetScalarD = engGetScalarNATIVE(epI, arrayS );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engGetVectorI: engGetVectorD = engGetVectorNATIVE(epI, arrayS );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engGetArrayI: engGetArrayD = engGetArrayNATIVE(epI, arrayS );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engGetCharArrayI: engGetCharArrayS = engGetCharArrayNATIVE(epI, arrayS );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engPutArray2dI: engPutArrayNATIVE( epI, arrayS, engPutArray2dD );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+ case engOutputBufferI: engOutputBufferS = engOutputBufferNATIVE( epI, buflenI );
+
+ releaseWaitForValue();
+
+ break;
+
+
+
+
+
+ default: //System.out.println("thread default switch statem.");
+
+ }
+
+ status=0;
+
+
+
+ lockThreadB = true;
+
+ while (lockThreadB == true) {
+
+ synchronized(this) {
+
+ try { wait();} // wait until next command is available
+
+ catch (InterruptedException e) { }
+
+ }
+
+ }
+
+ //System.out.println("JMatLink: thread awoke and passed lock");
+
+ if (destroyJMatLinkB == true) break;
+
+ } // end while
+
+ if (debugB) System.out.println("JMatLink: thread terminated");
+
+} // end run
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+// The MATLAB engine is served by a thread. Threads don't have methods
+
+// which can be called. So we need to send messages to that thread
+
+// by using notifyAll. In the meantime NO OTHER methods is allowed to
+
+// access our thread (engine) so we lock everything up.
+
+ private void lockEngineLock(){
+
+ synchronized(this){
+
+ while (lockEngineB==true){
+
+ try { //System.out.println("lockEngineLock locked");
+
+ wait();} // wait until last command is finished
+
+ catch (InterruptedException e) { }
+
+ }
+
+ //now lockEngineB is false
+
+ lockEngineB = true;
+
+ }
+
+ } // end lockEngine
+
+
+
+ private synchronized void releaseEngineLock(){
+
+ lockEngineB = false;
+
+ notifyAll();
+
+ }
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+// The MATLAB engine is served by a thread. Threads don't have methods
+
+// which can be called directly. If we send a command that returns data
+
+// back to the calling function e.g. engGetArray("array"), we'll notify
+
+// the main thread to get the data from matlab. Since the data is collected
+
+// in another thread, we don't know exactly when the data is available, since
+
+// this is a concurrent situation.
+
+// The solution is simple: I always use a locking-mechanism to wait for the
+
+// data. The main thread will release the lock and the calling method can
+
+// return the data.
+
+//
+
+// Steps:
+
+// 1. a method that returns data calls the locking method
+
+// 2. notify the thread to call matlab
+
+// 3. wait for the returned data
+
+// 4. after the thread itself got the data it releases the locks method
+
+// 5. return data
+
+
+
+ private synchronized void lockWaitForValue(){
+
+ lockWaitForValueB = true;
+
+ }
+
+
+
+ private void WaitForValue(){
+
+ synchronized(this){
+
+ while (lockWaitForValueB==true){
+
+ try { //System.out.println("lockWaitForValue locked");
+
+ wait();} // wait for return value
+
+ catch (InterruptedException e) { }
+
+ }
+
+ }
+
+ //System.out.println("WaitForValue released");
+
+ } // end waitForValue
+
+
+
+ private synchronized void releaseWaitForValue(){
+
+ lockWaitForValueB = false;
+
+ notifyAll();
+
+ }
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+//// Utility methods ////
+
+
+
+
+
+ // Convert an n*n double array to n*1 String vector
+
+ private String[] double2String(double[][] d)
+
+ {
+
+ String encodeS[]=new String[d.length]; // String vector
+
+
+
+ // for all rows
+
+ for (int n=0; n
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+