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

Importing release version 322 from old repos

Browse Source
This commit is contained in:
Marcel Kronfeld
2007-12-11 16:38:11 +00:00
parent 8cecdb016d
commit 7ae15be788
668 changed files with 109288 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

182
src/wsi/ra/math/Jama/CholeskyDecomposition.java Normal file
View File

@@ -0,0 +1,182 @@
package wsi.ra.math.Jama;
/** Cholesky Decomposition.
<P>
For a symmetric, positive definite matrix A, the Cholesky decomposition
is an lower triangular matrix L so that A = L*L'.
<P>
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.
// <P>
// 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);
}
}

956
src/wsi/ra/math/Jama/EigenvalueDecomposition.java Normal file
View File

@@ -0,0 +1,956 @@
package wsi.ra.math.Jama;
import wsi.ra.math.Jama.util.Maths;
/** Eigenvalues and eigenvectors of a real matrix.
<P>
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.
<P>
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;
}
}

318
src/wsi/ra/math/Jama/LUDecomposition.java Normal file
View File

@@ -0,0 +1,318 @@
package wsi.ra.math.Jama;
/** LU Decomposition.
<P>
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.
<P>
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.
<P>
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.
<P>
@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;
}
}

1111
src/wsi/ra/math/Jama/Matrix.java Normal file
View File

File diff suppressed because it is too large Load Diff

218
src/wsi/ra/math/Jama/QRDecomposition.java Normal file
View File

@@ -0,0 +1,218 @@
package wsi.ra.math.Jama;
import wsi.ra.math.Jama.util.Maths;
/** QR Decomposition.
<P>
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.
<P>
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));
}
}

540
src/wsi/ra/math/Jama/SingularValueDecomposition.java Normal file
View File

@@ -0,0 +1,540 @@
package wsi.ra.math.Jama;
import wsi.ra.math.Jama.util.*;
/** Singular Value Decomposition.
<P>
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'.
<P>
The singular values, sigma[k] = S[k][k], are ordered so that
sigma[0] >= sigma[1] >= ... >= sigma[n-1].
<P>
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<p
// kase = 2 if s(k) is negligible and k<p
// kase = 3 if e[k-1] is negligible, k<p, and
// s(k), ..., s(p) are not negligible (qr step).
// kase = 4 if e(p-1) is negligible (convergence).
for (k = p-2; 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;
}
}

21
src/wsi/ra/math/Jama/util/Maths.java Normal file
View File

@@ -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;
}
}

184
src/wsi/ra/math/RNG.java Normal file
View File

@@ -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<lo.length;i++)
xin[i] = (hi[i]-lo[i])*random.nextDouble()+lo[i];
return xin;
}
/**
*
*/
public static double[] randomDoubleArray(double lo,double hi,int size) {
double[] xin = new double[size];
for (int i=0;i<size;i++)
xin[i] = (hi-lo)*random.nextDouble()+lo;
return xin;
}
/**
*
*/
public static double[] randomDoubleArray(double[] lo,double[] hi,double[] xin) {
for (int i=0;i<lo.length;i++)
xin[i] = (hi[i]-lo[i])*random.nextDouble()+lo[i];
return xin;
}
/**
*
*/
public static boolean randomBoolean() {
return (randomInt()==1);
}
/**
*
*/
public static int randomBit() {
return randomInt();
}
/**
*
*/
public static boolean flipCoin(double p) {
return (randomDouble()<p ? true : false);
}
/**
*
*/
public static float gaussianFloat(float dev) {
return (float)random.nextGaussian()*dev;
}
/**
*
*/
public static double gaussianDouble(double dev) {
return random.nextGaussian()*dev;
}
/**
*
*/
public static float exponentialFloat(float mean) {
return (float)(-mean*Math.log(randomDouble()));
}
/**
*
*/
public static double exponentialDouble(double mean) {
return -mean*Math.log(randomDouble());
}
}

1364
src/wsi/ra/math/SpecialFunction.java Normal file
View File

File diff suppressed because it is too large Load Diff

52
src/wsi/ra/math/interpolation/AbstractDataModifier.java Normal file
View File

@@ -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
* <code>AbstractDataSet</code>.
*/
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
****************************************************************************/

116
src/wsi/ra/math/interpolation/AbstractDataSet.java Normal file
View File

@@ -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
****************************************************************************/

100
src/wsi/ra/math/interpolation/BasicDataSet.java Normal file
View File

@@ -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 <code>AbstractDataSet</code>.
*/
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
****************************************************************************/

44
src/wsi/ra/math/interpolation/InterpolationException.java Normal file
View File

@@ -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
****************************************************************************/

573
src/wsi/ra/math/interpolation/LinearInterpolation.java Normal file
View File

@@ -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 <code>AbstractDataSet</code>
*/
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
****************************************************************************/

453
src/wsi/ra/math/interpolation/PolyInterpolation.java Normal file
View File

@@ -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 <code>AbstractDataSet</code>
*/
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 <code>AbstractDataSet</code>
* @param sloppy if <code>true</code> Neville's algorithm which is used in the
* <code>polynomialInterpolation</code>-routines does only print a
* warning message on the screen and does not throw an
* <code>Exception</code> 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 <code>AbstractDataSet</code> and calculates the coefficients
* of the polynom.
*
* @param abstractDataSet the <code>AbstractDataSet</code>
*/
public void setAbstractDataSet(AbstractDataSet abstractDataSet)
throws InterpolationException
{
this.abstractDataSet = abstractDataSet;
this.polynomialCoefficients = calculatePolynomialCoefficients();
}
/**
* Uses the polynom with the calculated coefficients to calculate the
* <code>y</code> value. This algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
* The Neville's algorithm which is used in the <code>polynomialInterpolation</code>-
* 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
* <code>y</code> value and the derivatives at the point <code>x</code>,
* <code>y</code>. This algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
* The Neville's algorithm which is used in the <code>polynomialInterpolation</code>-
* 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:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*
* @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:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*
* @param abstractDataSet the <code>AbstractDataSet</code>
* @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:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*
* @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 <code>N-1</code>. This
* interpolation algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*
* @return the array with the polynomial coefficients y = ...[0] +
* ...[1]*x<SUP>2</SUP> + ...[2]*x<SUP>3</SUP> + ...
* @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 <code>N-1</code>. This
* interpolation algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*
* @param abstractDataSet the <code>AbstractDataSet</code>
* @return the array with the polynomial coefficients y = ...[0] +
* ...[1]*x<SUP>2</SUP> + ...[2]*x<SUP>3</SUP> + ...
* @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 <code>N-1</code>. This
* interpolation algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*
* @param x the array of x values
* @param y the array of y values
* @return the array with the polynomial coefficients y = ...[0] +
* ...[1]*x<SUP>2</SUP> + ...[2]*x<SUP>3</SUP> + ...
* @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
****************************************************************************/

56
src/wsi/ra/math/interpolation/PolynomialInterpolationResult.java Normal file
View File

@@ -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
****************************************************************************/

299
src/wsi/ra/math/interpolation/SplineInterpolation.java Normal file
View File

@@ -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 <code>AbstractDataSet</code>
*/
public SplineInterpolation(AbstractDataSet abstractDataSet)
throws InterpolationException
{
this.setAbstractDataSet(abstractDataSet);
}
/**
* Sets the new <code>AbstractDataSet</code> and calculates the second
* derivative of the spline.
*
* @param abstractDataSet the <code>AbstractDataSet</code>
*/
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 <code>y</code> value. This algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*/
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
* x<sub>i</sub> values of the function y<sub>i</sub>=f(x<sub>i</sub>) are
* in ascending order, as x<sub>0</sub>&lt;x<sub>1</sub>&lt;x<sub>2</sub>&lt;... .
* This algorithm was taken from:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*/
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:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*/
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:<br>
* <a href="http://www.ulib.org/webRoot/Books/Numerical_Recipes/" target="_top">
* 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.</a><br>
*/
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
****************************************************************************/