net.algart.math.functions.LinearOperator Maven / Gradle / Ivy
/*
* The MIT License (MIT)
*
* Copyright (c) 2007-2024 Daniel Alievsky, AlgART Laboratory (http://algart.net)
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
package net.algart.math.functions;
import net.algart.math.Point;
import java.util.Objects;
/**
* Linear operator (affine transformation):
* O f(x) = f(Ax + b),
* where the numeric n x n matrix A
* and the n-dimensional vector b are parameters of the transformation.
* (It is a particular case of the {@link ProjectiveOperator projective transformation},
* when c vector is zero and d number is 1.0.)
*
* This class is immutable and thread-safe:
* there are no ways to modify settings of the created instance.
*
* @author Daniel Alievsky
*/
public final class LinearOperator extends ProjectiveOperator {
LinearOperator(double[] a, double[] diagonal, double[] b) {
super(a, diagonal, b, null, 1.0);
}
/**
* Returns an instance of this class, describing the linear operator with the specified matrix A
* and vector b:
* O f(x) = f(Ax + b).
* The coordinates of the vector b must be listed in b array.
* The elements of the matrix A must be listed, row by row, in the a array:
* A={aij},
* aij=a[i*n+j],
* i is the index of the row (0..n-1),
* j is the index of the column (0..n-1),
* n=b.length.
* The length a.length of the a array must be equal to the square n2
* of the length n=b.length of the b array.
* Empty arrays (n=0) are not allowed.
*
* The passed a and b Java arrays are cloned by this method: no references to them
* are maintained by the created instance.
*
* @param a the elements of A matrix.
* @param b the coordinates of b vector.
* @return the linear operator O f(x) = f(Ax + b).
* @throws NullPointerException if one of the arguments of the method is {@code null}.
* @throws IllegalArgumentException if b.length==0 or a.length!=b.length2.
*/
public static LinearOperator getInstance(double[] a, double[] b) {
Objects.requireNonNull(a, "Null A matrix");
Objects.requireNonNull(b, "Null b vector");
return new LinearOperator(a.clone(), null, b.clone());
}
/**
* Returns an instance of this class, describing the linear operator with the diagonal matrix A
* and vector b:
* O f(x) = f(Ax + b),
* where A={aij},
* aij=0.0 if i!=j,
* aii=diagonal[i].
* The coordinates of the vector b must be listed in b array.
* Empty arrays (diagonal.length=b.length=0) are not allowed.
*
*
This linear operator performs resizing and shift along coordinate axes.
*
*
The passed diagonal and b Java arrays are cloned by this method: no references to them
* are maintained by the created instance.
*
* @param diagonal the diagonal elements of A matrix (all other elements are supposed to be zero).
* @param b the coordinates of b vector.
* @return the linear operator O f(x) = f(Ax + b).
* @throws NullPointerException if one of the arguments of the method is {@code null}.
* @throws IllegalArgumentException if diagonal.length==0 or diagonal.length!=b.length.
*/
public static LinearOperator getDiagonalInstance(double[] diagonal, double[] b) {
Objects.requireNonNull(diagonal, "Null diagonal array");
Objects.requireNonNull(b, "Null b vector");
return new LinearOperator(null, diagonal.clone(), b.clone());
}
/**
* Equivalent to
* {@link #getDiagonalInstance(double[], double[]) getDiagonalInstance(diagonal, new double[diagonal.length])}
* (the case of zero b vector).
* This linear operator performs resizing along coordinate axes.
*
* @param diagonal the diagonal elements of A matrix (all other elements are supposed to be zero).
* @return the linear operator O f(x) = f(Ax).
* @throws NullPointerException if the argument of the method is {@code null}.
* @throws IllegalArgumentException if diagonal.length==0.
*/
public static LinearOperator getDiagonalInstance(double... diagonal) {
Objects.requireNonNull(diagonal, "Null diagonal array");
return new LinearOperator(null, diagonal.clone(), new double[diagonal.length]);
}
/**
* Equivalent to
* {@link #getDiagonalInstance(double[], double[]) getDiagonalInstance(diagonal, b)},
* where diagonal is an array consisting of b.length unit values (1.0).
* This linear operator performs shifting along coordinate axes.
*
* @param b the coordinates of b vector.
* @return the linear operator O f(x) = f(x + b).
* @throws NullPointerException if the argument of the method is {@code null}.
* @throws IllegalArgumentException if b.length==0.
*/
public static LinearOperator getShiftInstance(double... b) {
Objects.requireNonNull(b, "Null b array");
return new LinearOperator(null, null, b);
}
/**
* Returns an instance of this class, describing rotation in 2D plane by the specified angle (in radians)
* around the specified center.
* Almost equivalent to {@link #getInstance(double[], double[]) getInstance}(a,b),
* where:
*
* a = {cos,sin,-sin,cos} — matrix A
* (cos=StrictMath.cos(angle)),
* sin=StrictMath.sin(angle));b = {centerX-a[0]*centerX-a[1]*centerY,
* centerY-a[2]*centerX-a[3]*centerY} — vector
* b=c−Ac (c={centerX, centerY}).
*
* The only difference from these formulas is special processing some cases, when the angle is kπ/2
* with integer k (more precisely, k/2.0*StrictMath.PI):
* StrictMath.cos and StrictMath.sin methods can return inexact results here,
* but this method uses precise values ±1 in these cases.
*
* @param centerX the x-coordinate of the rotation center.
* @param centerY the y-coordinate of the rotation center.
* @param angle the rotation angle (in radians; positive values correspond to clockwise rotation,
* if the x axis is directed rightwards and the y axis
* is directed downwards, according traditions of computer image processing).
* @return 2-dimensional linear operator, describing this rotation.
*/
public static LinearOperator getRotation2D(double centerX, double centerY, double angle) {
double cos = StrictMath.cos(angle);
double sin = StrictMath.sin(angle);
for (int k = -4; k <= 4; k++) {
if (angle == k / 2.0 * StrictMath.PI) { // use precise cosine and sine in these cases
assert StrictMath.abs(StrictMath.round(cos) - cos) < 1e-6;
assert StrictMath.abs(StrictMath.round(sin) - sin) < 1e-6;
cos = StrictMath.round(cos);
sin = StrictMath.round(sin);
}
}
// y=A(x-c)+c=Ax+c-Ac
double[] a = {
cos, sin,
-sin, cos,
};
double[] b = {
centerX - a[0] * centerX - a[1] * centerY,
centerY - a[2] * centerX - a[3] * centerY,
};
return new LinearOperator(a, null, b);
}
/**
* Returns the n-dimensional linear operator, that transforms (maps)
* the given n+1 points p0, p1, ..., pn to
* the given another n+1 points q0, q1, ..., qn
* of the n-dimensional space.
* In other words, the matrix A and the vector b in the returned operator
* fulfil the following n+1 conditions:
*
*
* Ap0 + b = q0,
* Ap1 + b = q1,
* ...,
* Apn + b = qn
*
*
* It is possible that there is no such operator
* or there are many different solutions (degenerated cases).
* In this case, this method still returns some operator, but some coefficients of A matrix
* and b vector in the returned operator will probably be Double.NaN,
* Double.POSITIVE_INFINITY or Double.NEGATIVE_INFINITY.
*
*
All passed points must be n-dimensional,
* where n+1=p.length=q.length.
*
* @param q the destination points.
* @param p the source points.
* @return the n-dimensional linear operator, which maps pi to
* qi for all i=0,1,2,...,n.
* @throws NullPointerException if one of arguments of this method or one of elements of
* p and q arrays {@code null}.
* @throws IllegalArgumentException if the lengths of the passed p and
* q arrays are not equal,
* or if for some k
* p[k].{@link Point#coordCount() coordCount()}!=p.length-1 or
* q[k].{@link Point#coordCount() coordCount()}!=p.length-1.
* @throws OutOfMemoryError if there is not enough Java memory for storing two Java arrays
* double[n*n] and double[(n+1)*(n+1)],
* where n+1=p.length,
* or if (n+1)*(n+1)>Integer.MAX_VALUE.
*/
public static LinearOperator getInstanceByPoints(Point[] q, Point[] p) {
Objects.requireNonNull(p, "Null p argument");
Objects.requireNonNull(q, "Null q argument");
if (p.length != q.length)
throw new IllegalArgumentException("p and q point arrays lengths mismatch: p.length="
+ p.length + ", q.length=" + q.length);
if (p.length == 0)
throw new IllegalArgumentException("Empty p and q arrays");
final int n = p.length - 1;
if (((long) n + 1) * ((long) n + 1) > Integer.MAX_VALUE)
throw new OutOfMemoryError("Too large necessary matrix (more than Integer.MAX_VALUE elements)");
for (int k = 0; k < p.length; k++) {
if (p[k].coordCount() != n)
throw new IllegalArgumentException("n+1 n-dimensional points are necessary to "
+ "find the linear operator, but we have " + (n + 1) + " points, "
+ "and the source point #" + k + " is " + p[k].coordCount() + "-dimensional");
if (q[k].coordCount() != n)
throw new IllegalArgumentException("n+1 n-dimensional points are necessary to "
+ "find the linear operator, but we have " + (n + 1) + " points, "
+ "and the destination point #" + k + " is " + q[k].coordCount() + "-dimensional");
}
// A * (px0,py0,pz0,...) + b = (qx0,qy0,qz0,...)
// A * (px1,py1,pz1,...) + b = (qx1,qy1,qz1,...)
// ..., i.e.
//
// px0*a00 + py0*a01 + pz0*a02 + ... + bx = qx0
// px1*a00 + py1*a01 + pz1*a02 + ... + bx = qx1
// px2*a00 + py2*a01 + pz2*a02 + ... + bx = qx2
// ...
// px0*a10 + py0*a11 + pz0*a12 + ... + by = qy0
// px1*a10 + py1*a11 + pz1*a12 + ... + by = qy1
// px2*a10 + py2*a11 + pz2*a13 + ... + by = qy2
// etc.
// In other words, here are n independent equation systems Sv=t
double[] s = new double[(n + 1) * (n + 1)];
double[] t = new double[n + 1];
double[] v = new double[n + 1];
for (int i = 0, sOfs = 0; i <= n; i++) {
for (int j = 0; j < n; j++, sOfs++) {
s[sOfs] = p[i].coord(j);
}
s[sOfs++] = 1.0;
}
// S is the matrix of the equations set (same for all equations sets)
double[] a = new double[n * n];
double[] b = new double[n];
for (int i = 0, aOfs = 0; i < n; i++, aOfs += n) {
// we shall solve here equations set #i
// and find the line #i of the matrix A and the element #i of the vector b
for (int j = 0; j <= n; j++) {
t[j] = q[j].coord(i);
}
// t is the right side of the equations set
solveLinearEquationsSet(v, s.clone(), t);
// solveLinearEquationsSetByCramerMethod(v, s, t);
System.arraycopy(v, 0, a, aOfs, n);
b[i] = v[n];
}
return new LinearOperator(a, null, b);
}
/**
* Returns superposition of this and the passed linear operators.
* More precisely, if this operator corresponds to the affine transformation Ax + b,
* and the passed one corresponds to the affine transformation A'x + b',
* then the returned operator corresponds to the affine transformation
* A''x + b'' = A'(Ax + b) + b',
* i.e. in the returned operator
* A'' = A'A, b'' = A'b + b'.
*
* @param operator the second operator, that should be applied after this one.
* @return superposition of this and the passed operator.
* @throws NullPointerException if the argument of the method is {@code null}.
* @throws IllegalArgumentException if operator.{@link #n() n()}!=this.{@link #n() n()}.
*/
public LinearOperator superposition(LinearOperator operator) {
Objects.requireNonNull(operator, "Null operator argument");
if (operator.n != n)
throw new IllegalArgumentException("Passed and this operators dimensions mismatch: operator.n()="
+ operator.n + ", this.n()=" + n);
if (operator.isShift()) {
double[] b = this.b.clone();
for (int i = 0; i < n; i++) {
b[i] += operator.b[i];
}
return new LinearOperator(this.a, this.diagonal, b);
}
double[] a = this.a != null ? this.a : this.a(); // A
double[] aOther = operator.a != null ? operator.a : operator.a(); // A'
double[] aNew = new double[a.length]; // A'' (result)
double[] bNew = new double[n]; // b'' (result)
// y'' = A'y' + b' = A'(Ax+b) + b' = A'Ax + (A'b+b'), so A''=A'A, b''=A'b+b'
for (int in = 0; in < a.length; in += n) { // in = i * n
for (int j = 0; j < n; j++) {
double sum = 0.0;
for (int k = 0, disp = j; k < n; k++, disp += n) {
sum += aOther[in + k] * a[disp]; // aOther[i * n + k] * a[k * n + j];
}
aNew[in + j] = sum;
}
}
for (int i = 0, in = 0; i < n; i++, in += n) { // in = i * n
double sum = operator.b[i];
for (int j = 0; j < n; j++) {
sum += aOther[in + j] * b[j];
}
bNew[i] = sum;
}
return new LinearOperator(aNew, null, bNew);
}
/**
* Returns an instance of this class, identical to this one execpting that the new instance has the specified
* vector b.
*
*
The passed b Java array is cloned by this method: no references to it
* are maintained by the created instance.
*
* @param b the new coordinates of b vector.
* @return the linear operator with changed b vector.
* @throws NullPointerException if the argument of the method is {@code null}.
* @throws IllegalArgumentException if b.length!=this.{@link #n() n()}.
*/
public LinearOperator changeB(double... b) {
Objects.requireNonNull(b, "Null b array");
if (b.length != n)
throw new IllegalArgumentException("Passed b and this.b vector lengths mismatch: b.length="
+ b.length + ", this b.length=" + n);
return new LinearOperator(this.a, this.diagonal, b);
}
/**
* This implementation calculates destPoint by multiplication
* the srcPoint by the matrix A and adding the vector b.
* to the coordinates destPoint of the destination point.
*
* @param destPoint the coordinates of the destination point y, filled by this method.
* @param srcPoint the coordinates of the source point x.
* @throws NullPointerException if one of the arguments is {@code null}.
* @throws IllegalArgumentException if destPoint.length or srcPoint.length
* is not equal to the {@link #n() number of dimensions}.
*/
public void map(double[] destPoint, double[] srcPoint) {
calculateAxPlusB(destPoint, srcPoint);
}
/**
* Transforms the coordinates destPoint of the destination point in n-dimensional space
* back to the coordinates srcPoint of the original point.
*
*
To find the srcPoint, this method solves the system of linear equations Ax=y−b,
* where the matrix A and the vector b are the parameters of this transformation,
* y is destPoint, x is srcPoint.
* This method uses
* Gauss elimination algorithm
* with partial (column) pivoting.
* This algorithm requires O(n3) operations.
*
*
It is possible that there is no required srcPoint
* or there are many different solutions (degenerated cases).
* In this case, this method still returns some point, but some found srcPoint coordinates
* will probably be Double.NaN, Double.POSITIVE_INFINITY
* or Double.NEGATIVE_INFINITY.
*
*
Note: this method allocates some additional memory if the matrix A
* is not {@link #isDiagonal() diagonal}.
* If you don't want to occupy additional memory, you can directly use
* {@link #solveLinearEquationsSet(double[], double[], double[])} method.
*
*
Note: this method works correctly even if destPoint and srcPoint
* is the same Java array.
*
* @param srcPoint the coordinates of the source point x, filled by this method.
* @param destPoint the coordinates of the destinated point y.
* @throws NullPointerException if one of the arguments is {@code null}.
* @throws IllegalArgumentException if destPoint.length or srcPoint.length
* is not equal to the {@link #n() number of dimensions}.
*/
public void inverseMap(double[] srcPoint, double[] destPoint) {
Objects.requireNonNull(srcPoint, "Null srcPoint");
Objects.requireNonNull(destPoint, "Null destPoint");
if (srcPoint.length != n)
throw new IllegalArgumentException("Illegal length of srcPoint array: "
+ srcPoint.length + " for " + this);
if (destPoint.length != n)
throw new IllegalArgumentException("Illegal length of destPoint array: "
+ destPoint.length + " for " + this);
if (a != null) {
double[] yMinusB = new double[n];
for (int i = 0; i < n; i++) {
yMinusB[i] = destPoint[i] - b[i];
}
double[] aClone = this.a.clone();
solveLinearEquationsSet(srcPoint, aClone, yMinusB);
} else if (diagonal != null) {
for (int i = 0; i < n; i++) {
srcPoint[i] = (destPoint[i] - b[i]) / diagonal[i];
}
} else {
for (int i = 0; i < n; i++) {
srcPoint[i] = destPoint[i] - b[i];
}
}
}
/**
* Solves the system of linear equations Ax=y by
* Gauss elimination algorithm
* with partial (column) pivoting.
*
*
The coordinates of the vector y must be listed in y array.
* The elements of the matrix A must be listed, row by row, in the a array:
* A={aij},
* aij=a[i*n+j],
* i is the index of the row (0..n-1),
* j is the index of the column (0..n-1),
* n=b.length.
* The length a.length of the a array must be equal to the square n2
* of the length n=b.length of the b array.
* Empty arrays (n=0) are not allowed.
*
*
It is possible that there is no required x vector
* or there are many different solutions (degenerated cases).
* In this case, this method still find some x vector, but some found coordinates
* in the x array will probably be
* Double.NaN, Double.POSITIVE_INFINITY or Double.NEGATIVE_INFINITY.
*
*
This method is called in the {@link #inverseMap(double[], double[])} method,
* if the matrix A is not {@link #isDiagonal() diagonal}.
*
*
Warning: this method destroys the content of the passed a and y arrays!
* But this method does not occupy any additional memory, unlike {@link #inverseMap(double[], double[])} method.
*
*
Warning: this method will not work correctly if x and y
* is the same Java array.
*
* @param x the coordinates of x vector, filled by this method.
* @param a the elements of A matrix (row by row).
* @param y the coordinates of y vector.
* @throws NullPointerException if one of the arguments of the method is {@code null}.
* @throws IllegalArgumentException if the length of one of the passed arrays is 0,
* or if x.length!=y.length,
* or if a.length!=x.length2.
*/
public static void solveLinearEquationsSet(double[] x, double[] a, double[] y) {
Objects.requireNonNull(x, "Null x");
Objects.requireNonNull(y, "Null y");
Objects.requireNonNull(a, "Null a");
final int n = x.length;
if (y.length != n)
throw new IllegalArgumentException("x and y vector lengths mismatch: x.length="
+ n + ", y.length=" + y.length);
if (a.length != (long) n * (long) n)
throw new IllegalArgumentException("Illegal size of A matrix: a.length=" + a.length
+ " must be equal to x.length^2=" + (long) n * (long) n);
// System.out.print("Gauss solving " + n + "x" + n + "... ");
// long t1 = System.nanoTime();
// printEquationsSet(null, a, y);
// Elimination algorithm
for (int k = 0; k < n - 1; k++) {
// Let, for example, k=2. Now we have the following equations set (below m=this.n-1):
// a00*x0 + a01*x1 + a02*x2 + a03*x3 + ... + a0m*xm = y0
// a11*x1 + a12*x2 + a13*x3 + ... + a1m*xm = y1
// a22*x2 + a23*x3 + ... + a2m*xm = y2
// a32*x2 + a33*x3 + ... + a3m*xm = y3
// a42*x2 + a43*x3 + ... + a4m*xm = y4
// . . .
// am2*x2 + am3*x3 + ... + amm*xm = ym
// Finding the pivot - the maximal |aik|, i=k,k+1,...,n-1:
final int aOfsDiagonal = k * n + k;
int pivotIndex = k;
double pivot = a[aOfsDiagonal], pivotAbs = StrictMath.abs(pivot);
for (int i = k + 1, aOfs = aOfsDiagonal + n; i < n; i++, aOfs += n) {
double v = a[aOfs], vAbs = StrictMath.abs(v);
if (vAbs > pivotAbs) {
pivot = v;
pivotAbs = vAbs;
pivotIndex = i;
}
}
if (pivotAbs == 0.0) {
// Little optimization: avoiding useless futher calculations
for (int j = 0; j < n; j++) {
x[j] = Double.NaN;
}
return;
}
// Exchanging the line #k and line #i: this operation does not change the order of the unknowns
if (pivotIndex != k) {
double temp = y[k];
y[k] = y[pivotIndex];
y[pivotIndex] = temp;
for (int j = k, aOfsK = aOfsDiagonal, aOfs = pivotIndex * n + k; j < n; j++, aOfsK++, aOfs++) {
// Don't exchange first k elements in the lines: they are already zero
temp = a[aOfsK];
a[aOfsK] = a[aOfs];
a[aOfs] = temp;
}
}
assert a[aOfsDiagonal] == pivot : "Pivot element is not placed to the correct place (k = " + k + ")";
// Now |aik|<=|akk| for all i>k.
// Eliminating - subtracting the line #k, multiplied by aik/akk, from all lines #i
for (int i = k + 1, aOfs = k * n + n; i < n; i++) { // aOfs refers to the line begin!
double q = a[aOfs + k] / pivot;
y[i] -= y[k] * q;
for (int j = 0; j < k; j++, aOfs++) {
assert a[aOfs] == 0.0 : "The line was not filled by zero in previous iterations 0,...," + (k - 1);
}
a[aOfs++] = 0.0; // not necessary, but the resulting matrix will be "better" (triangular)
int aOfsK = aOfsDiagonal + 1;
assert aOfs == i * n + k + 1;
assert aOfsK == k * n + k + 1;
for (int j = k + 1; j < n; j++, aOfsK++, aOfs++) {
a[aOfs] -= a[aOfsK] * q;
}
}
}
// Now A is triangular:
// a00*x0 + a01*x1 + a02*x2 + a03*x3 + ... + a0m*xm = y0
// a11*x1 + a12*x2 + a13*x3 + ... + a1m*xm = y1
// a22*x2 + a23*x3 + ... + a2m*xm = y2
// a33*x3 + ... + a3m*xm = y3
// . . .
// amm*xm = ym
// Finding the unknowns (x):
for (int i = n - 1, aOfs = n * n - 1; i >= 0; i--) {
assert aOfs == i * n + n - 1;
double v = y[i];
for (int j = n - 1; j > i; j--, aOfs--) {
v -= a[aOfs] * x[j];
}
x[i] = v / a[aOfs];
aOfs -= i + 1;
}
// printEquationsSet(x, a, y);
// long t2 = System.nanoTime();
// System.out.println("done in " + (t2 - t1) * 1e-3 + " mcs = n^3*" + (t2 - t1) / Math.pow(n, 3) + " ns");
}
// Not used in current version. Very slow solution for matrices >3x3.
// The private "determinant" method is moved into Matrices class.
// static void solveLinearEquationsSetByCramerMethod(double[] x, double[] a, double[] y) {
// final int n = y.length;
// long t1 = System.nanoTime();
// System.out.print("Cramer solving " + n + "x" + n + "...");
// double det = determinant(n, a);
// double[] m = new double[a.length];
// for (int k = 0; k < n; k++) {
// for (int i = 0, aOfs = 0, mOfs = 0; i < n; i++) {
// for (int j = 0; j < n; j++, aOfs++, mOfs++) {
// m[mOfs] = j == k ? y[i] : a[aOfs];
// }
// }
// x[k] = determinant(n, m) / det;
// }
// long t2 = System.nanoTime();
// System.out.println(" done in " + (t2 - t1) * 1e-3 + " mcs = n^3*" + (t2 - t1) / Math.pow(n, 3) + " ns");
// }
/**
* Returns a brief string description of this object.
*
* @return a brief string description of this object.
*/
public String toString() {
boolean shift = isShift();
StringBuilder sA = new StringBuilder();
if (isDiagonal()) {
if (diagonal != null) {
sA.append("diag[");
for (int k = 0; k < diagonal.length; k++) {
if (k > 0)
sA.append(",");
sA.append(LinearFunc.goodFormat(diagonal[k]));
}
sA.append("]");
}
} else {
sA.append("A");
}
StringBuilder sB = new StringBuilder();
for (int k = 0; k < n; k++) {
if (k > 0)
sB.append(",");
sB.append(LinearFunc.goodFormat(b[k]));
}
return "linear (affine) " + n + "-dimensional operator "
+ sA + "x+b, b=(" + sB + ")" + (shift ? " (shift)" : "");
}
static void printEquationsSet(double[] x, double[] a, double[] y) {
final int n = y.length;
System.out.println();
for (int i = 0; i < n; i++) {
double v = 0.0;
for (int j = 0; j < n; j++) {
if (x != null) {
v += a[i * n + j] * x[j];
}
System.out.printf("%10.6f*x[%d] " + (j < n - 1 ? "+ " : "= "), a[i * n + j], j);
}
System.out.printf("%.7f", y[i]);
if (x != null)
System.out.printf(" (error %.7f)", v - y[i]);
System.out.println();
}
}
/**
* The simplest test for this class: finds a linear operator by pairs of points.
*/
static class Test extends ProjectiveOperator.Test {
@Override
int numberOfPoints(int dimCount) {
return dimCount + 1;
}
@Override
CoordinateTransformationOperator getOperator() {
return LinearOperator.getInstanceByPoints(q, p);
}
void badLinearEquationsSetTest() {
double factorial = 1.0;
for (int k = 2; k <= 16; k++) {
factorial *= k;
}
double[] a = new double[dimCount * dimCount];
for (int i = 1, aOfs = 0; i <= dimCount; i++) {
for (int j = 1; j <= dimCount; j++, aOfs++) {
a[aOfs] = factorial / (i + j);
}
}
LinearOperator lo = LinearOperator.getInstance(a, new double[dimCount]);
newRndPoints(10, dimCount);
for (int k = 0; k < p.length; k++) {
double[] srcPoint = p[k].coordinates();
double[] destPoint = new double[srcPoint.length];
q[k] = p[k];
lo.map(destPoint, srcPoint);
lo.inverseMap(srcPoint, destPoint);
r[k] = Point.valueOf(srcPoint);
}
double error = maxDiff();
System.out.println("Difference for \"bad\" equation set " + dimCount + "*" + dimCount + " = " + error);
LinearOperator lo2 = LinearOperator.getInstance(a, new double[dimCount]);
if (!lo.equals(lo2))
throw new AssertionError("Error in equals");
ProjectiveOperator po = ProjectiveOperator.getInstance(a, lo.b(), new double[dimCount], 1.0);
if (!lo.equals(po))
throw new AssertionError("Error in equals");
}
public static void main(String[] args) {
Test test = new Test();
test.init(args);
test.badLinearEquationsSetTest();
test.mainTest();
}
}
}