use of javax.vecmath.GMatrix in project bb4-common by bb4.
the class ConjugateGradientSolver method solve.
/**
* Find a solution or return the initial guess if something goes wrong.
* @param initialGuess the initial guess for the solution x, x0
* @return solution vector
*/
public GVector solve(GVector initialGuess) {
GVector x = new GVector(initialGuess);
GVector tempv = new GVector(initialGuess);
tempv.mul(matrix, initialGuess);
GVector bb = new GVector(b);
bb.sub(tempv);
GVector r = new GVector(bb);
GVector p = new GVector(r);
GVector xnew = new GVector(p);
GVector rnew = new GVector(p);
GVector pnew = new GVector(p);
GVector matrixMultp = new GVector(p);
GMatrix matrixInverse = new GMatrix(matrix);
matrixInverse.invert();
double error, norm;
int iteration = 0;
do {
matrixMultp.mul(matrix, p);
double lambda = (r.dot(p) / p.dot(matrixMultp));
xnew.scaleAdd(lambda, p, x);
rnew.scaleAdd(-lambda, matrixMultp, r);
double alpha = -(rnew.dot(matrixMultp) / p.dot(matrixMultp));
pnew.scaleAdd(alpha, p, rnew);
p.set(pnew);
r.set(rnew);
// System.out.println("the residual = "+r.toString());
x.set(xnew);
// error = Math.abs(r.dot(r)); // wrong way to compute norm
rnew.mul(r, matrixInverse);
norm = rnew.dot(r);
error = norm * norm;
// System.out.println("xi = "+x.toString());
iteration++;
// System.out.println("The error for iteration " + iteration + " is : " + error );
} while (error > eps && iteration < maxIterations);
if (error > eps || Double.isNaN(error) || Double.isInfinite(error)) {
// something went wrong
throw new IllegalStateException("Unable to converge on a solution. Error = " + error);
// return initialGuess;
}
return xnew;
}
use of javax.vecmath.GMatrix in project bb4-common by bb4.
the class ConjugateGradientSolverTest method solve3by4System.
/* didn't work. typo?
@Test
public void solveSimple3by3System() {
GMatrix A = new GMatrix(3, 3, new double[] {
5, 4, -1,
0, 10, -3,
0, 0, 1
});
GVector b = new GVector(new double[] {0, 11, 3});
solver = new ConjugateGradientSolver(A, b);
GVector solution = solver.solve();
System.out.println("Solution = "+ solution);
GVector expectedSolution = new GVector(new double[] {-1, 2, 3});
assertTrue("Unexpected solution:" + solution,
LinearUtil.appxVectorsEqual(expectedSolution, solution, 0.0000001));
}*/
// Conjugant gradient can only be applied to square matrices
@Test(expected = MismatchedSizeException.class)
public void solve3by4System() {
GMatrix A = new GMatrix(3, 4, new double[] { 2, 2, -1, 1, 1, 1, 1, 1, 2, -2, 2, 3 });
GVector b = new GVector(new double[] { 11, 4, 7 });
solver = new ConjugateGradientSolver(A, b);
solver.solve();
}
use of javax.vecmath.GMatrix in project imageio-ext by geosolutions-it.
the class NetCDFCFExperiment method Resampler.
private WritableRaster Resampler(final Array latData, final Array lonData, final int imageWidth, final int imageHeight, final int polyDegree, final Array data, final float fillValue) {
final Index latIndex = latData.getIndex();
final Index lonIndex = lonData.getIndex();
final int numCoeffs = (polyDegree + 1) * (polyDegree + 2) / 2;
final int XOFFSET = 0;
final int YOFFSET = 1;
final int stepX = 2;
final int stepY = 2;
int numNeededPoints = 0;
for (int xi = 0; xi < imageWidth; xi += stepX) {
for (int yi = 0; yi < imageHeight; yi += stepY) {
numNeededPoints++;
}
}
computeMatrixExtremes(latData, lonData, imageWidth, imageHeight, latIndex, lonIndex);
float[] destCoords = new float[2 * numNeededPoints];
float[] srcCoords = new float[2 * numNeededPoints];
/*
* Copy source and destination coordinates into float arrays. The
* destination coordinates are scaled in order to gets values similar to
* source coordinates (values will be identical if all "real world"
* coordinates are grid indices multiplied by a constant).
*/
int offset = 0;
for (int yi = 0; yi < imageHeight; yi += stepY) {
for (int xi = 0; xi < imageWidth; xi += stepX) {
srcCoords[offset] = xi;
srcCoords[offset + 1] = yi;
destCoords[offset] = (float) ((lonData.getFloat(lonIndex.set(xi)) - this.xmin) / this.periodX);
destCoords[offset + 1] = (float) ((this.ymax - latData.getFloat(latIndex.set(yi))) / this.periodY);
// destCoords[offset + 1] = ((latData.getFloat(latIndex.set(yi)) - this.ymin) / this.periodY);
offset += 2;
}
}
GMatrix A = new GMatrix(numNeededPoints, numCoeffs);
for (int coord = 0; coord < numNeededPoints; coord++) {
int var = 0;
for (int i = 0; i <= polyDegree; i++) {
for (int j = 0; j <= i; j++) {
double value = Math.pow(destCoords[2 * coord + XOFFSET], (double) (i - j)) * Math.pow(destCoords[2 * coord + YOFFSET], (double) j);
A.setElement(coord, var++, value);
}
}
}
GMatrix AtAi = new GMatrix(numCoeffs, numCoeffs);
GMatrix Ap = new GMatrix(numCoeffs, numNeededPoints);
AtAi.mulTransposeLeft(A, A);
AtAi.invert();
Ap.mulTransposeRight(AtAi, A);
GMatrix xVector = new GMatrix(numNeededPoints, 1);
GMatrix yVector = new GMatrix(numNeededPoints, 1);
for (int idx = 0; idx < numNeededPoints; idx++) {
xVector.setElement(idx, 0, srcCoords[2 * idx + XOFFSET]);
yVector.setElement(idx, 0, srcCoords[2 * idx + YOFFSET]);
}
GMatrix xCoeffsG = new GMatrix(numCoeffs, 1);
GMatrix yCoeffsG = new GMatrix(numCoeffs, 1);
xCoeffsG.mul(Ap, xVector);
yCoeffsG.mul(Ap, yVector);
float[] xCoeffs = new float[numCoeffs];
float[] yCoeffs = new float[numCoeffs];
for (int ii = 0; ii < numCoeffs; ii++) {
xCoeffs[ii] = new Double(xCoeffsG.getElement(ii, 0)).floatValue();
yCoeffs[ii] = new Double(yCoeffsG.getElement(ii, 0)).floatValue();
}
WritableRaster outDataCube;
WritableRandomIter iteratorDataCube;
SampleModel outSampleModel = RasterFactory.createBandedSampleModel(// data type
DataBuffer.TYPE_FLOAT, // width
imageWidth, // height
imageHeight, // num bands
1);
outDataCube = Raster.createWritableRaster(outSampleModel, null);
iteratorDataCube = RandomIterFactory.createWritable(outDataCube, null);
// Transfering data in the WritableRaster structure
Index indexInputVar = data.getIndex();
for (int jj = 0; jj < outDataCube.getNumBands(); jj++) {
for (int kk = 0; kk < outDataCube.getWidth(); kk++) {
for (int ll = 0; ll < outDataCube.getHeight(); ll++) {
iteratorDataCube.setSample(kk, ll, jj, data.getFloat(indexInputVar.set(ll, kk)));
}
}
}
WritableRaster target = RasterFactory.createWritableRaster(outSampleModel, null);
for (int bi = 0; bi < outDataCube.getNumBands(); bi++) {
for (int yi = 0; yi < imageHeight; yi++) {
for (int xi = 0; xi < imageWidth; xi++) {
float[] dstCoords = new float[2];
GMatrix regressionVec = new GMatrix(numCoeffs, 1);
int var = 0;
for (int i = 0; i <= polyDegree; i++) {
for (int j = 0; j <= i; j++) {
double value = Math.pow(xi, (double) (i - j)) * Math.pow(yi, (double) j);
regressionVec.setElement(var++, 0, value);
}
}
GMatrix xG = new GMatrix(1, 1);
GMatrix yG = new GMatrix(1, 1);
xG.mulTransposeLeft(regressionVec, xCoeffsG);
yG.mulTransposeLeft(regressionVec, yCoeffsG);
int X = (int) Math.round(xG.getElement(0, 0));
int Y = (int) Math.round(yG.getElement(0, 0));
if (X >= 0 && Y >= 0 && X < imageWidth && Y < imageHeight) {
target.setSample(xi, yi, bi, outDataCube.getSampleFloat(X, Y, bi));
} else {
// TODO: Change with fillvalue
// target.setSample(xi, yi, bi, Float.NaN);
target.setSample(xi, yi, bi, fillValue);
}
}
}
}
return target;
}
use of javax.vecmath.GMatrix in project bb4-common by bb4.
the class ConjugateGradientSolverTest method solveSimple2by2System.
@Test
public void solveSimple2by2System() {
GMatrix A = new GMatrix(2, 2, new double[] { 4, 1, 1, 3 });
GVector b = new GVector(new double[] { 1, 2 });
solver = new ConjugateGradientSolver(A, b);
GVector solution = solver.solve();
System.out.println("Solution = " + solution);
GVector expectedSolution = new GVector(new double[] { 1.0 / 11.0, 7.0 / 11.0 });
assertTrue("Unexpected solution:" + solution, LinearUtil.appxVectorsEqual(expectedSolution, solution, 0.000001));
}
use of javax.vecmath.GMatrix in project imageio-ext by geosolutions-it.
the class HOPSConverter method Resampler.
private WritableRaster Resampler(final Array latData, final Array lonData, final int imageWidth, final int imageHeight, final int polyDegree, final Array data, final float fillValue) {
final Index latIndex = latData.getIndex();
final Index lonIndex = lonData.getIndex();
final int numCoeffs = (polyDegree + 1) * (polyDegree + 2) / 2;
final int XOFFSET = 0;
final int YOFFSET = 1;
final int stepX = 2;
final int stepY = 2;
int numNeededPoints = 0;
for (int xi = 0; xi < imageWidth; xi += stepX) {
for (int yi = 0; yi < imageHeight; yi += stepY) {
numNeededPoints++;
}
}
computeMatrixExtremes(latData, lonData, imageWidth, imageHeight, latIndex, lonIndex);
float[] destCoords = new float[2 * numNeededPoints];
float[] srcCoords = new float[2 * numNeededPoints];
/*
* Copy source and destination coordinates into float arrays. The
* destination coordinates are scaled in order to gets values similar to
* source coordinates (values will be identical if all "real world"
* coordinates are grid indices multiplied by a constant).
*/
int offset = 0;
for (int yi = 0; yi < imageHeight; yi += stepY) {
for (int xi = 0; xi < imageWidth; xi += stepX) {
srcCoords[offset] = xi;
srcCoords[offset + 1] = yi;
destCoords[offset] = (float) ((lonData.getFloat(lonIndex.set(xi)) - this.xmin) / this.periodX);
destCoords[offset + 1] = (float) ((this.ymax - latData.getFloat(latIndex.set(yi))) / this.periodY);
// destCoords[offset + 1] = ((latData.getFloat(latIndex.set(yi)) - this.ymin) / this.periodY);
offset += 2;
}
}
GMatrix A = new GMatrix(numNeededPoints, numCoeffs);
for (int coord = 0; coord < numNeededPoints; coord++) {
int var = 0;
for (int i = 0; i <= polyDegree; i++) {
for (int j = 0; j <= i; j++) {
double value = Math.pow(destCoords[2 * coord + XOFFSET], (double) (i - j)) * Math.pow(destCoords[2 * coord + YOFFSET], (double) j);
A.setElement(coord, var++, value);
}
}
}
GMatrix AtAi = new GMatrix(numCoeffs, numCoeffs);
GMatrix Ap = new GMatrix(numCoeffs, numNeededPoints);
AtAi.mulTransposeLeft(A, A);
AtAi.invert();
Ap.mulTransposeRight(AtAi, A);
GMatrix xVector = new GMatrix(numNeededPoints, 1);
GMatrix yVector = new GMatrix(numNeededPoints, 1);
for (int idx = 0; idx < numNeededPoints; idx++) {
xVector.setElement(idx, 0, srcCoords[2 * idx + XOFFSET]);
yVector.setElement(idx, 0, srcCoords[2 * idx + YOFFSET]);
}
GMatrix xCoeffsG = new GMatrix(numCoeffs, 1);
GMatrix yCoeffsG = new GMatrix(numCoeffs, 1);
xCoeffsG.mul(Ap, xVector);
yCoeffsG.mul(Ap, yVector);
float[] xCoeffs = new float[numCoeffs];
float[] yCoeffs = new float[numCoeffs];
for (int ii = 0; ii < numCoeffs; ii++) {
xCoeffs[ii] = new Double(xCoeffsG.getElement(ii, 0)).floatValue();
yCoeffs[ii] = new Double(yCoeffsG.getElement(ii, 0)).floatValue();
}
WritableRaster outDataCube;
WritableRandomIter iteratorDataCube;
SampleModel outSampleModel = RasterFactory.createBandedSampleModel(// data type
DataBuffer.TYPE_FLOAT, // width
imageWidth, // height
imageHeight, // num bands
1);
outDataCube = Raster.createWritableRaster(outSampleModel, null);
iteratorDataCube = RandomIterFactory.createWritable(outDataCube, null);
// Transfering data in the WritableRaster structure
Index indexInputVar = data.getIndex();
for (int jj = 0; jj < outDataCube.getNumBands(); jj++) {
for (int kk = 0; kk < outDataCube.getWidth(); kk++) {
for (int ll = 0; ll < outDataCube.getHeight(); ll++) {
iteratorDataCube.setSample(kk, ll, jj, data.getFloat(indexInputVar.set(ll, kk)));
}
}
}
WritableRaster target = RasterFactory.createWritableRaster(outSampleModel, null);
for (int bi = 0; bi < outDataCube.getNumBands(); bi++) {
for (int yi = 0; yi < imageHeight; yi++) {
for (int xi = 0; xi < imageWidth; xi++) {
float[] dstCoords = new float[2];
GMatrix regressionVec = new GMatrix(numCoeffs, 1);
int var = 0;
for (int i = 0; i <= polyDegree; i++) {
for (int j = 0; j <= i; j++) {
double value = Math.pow(xi, (double) (i - j)) * Math.pow(yi, (double) j);
regressionVec.setElement(var++, 0, value);
}
}
GMatrix xG = new GMatrix(1, 1);
GMatrix yG = new GMatrix(1, 1);
xG.mulTransposeLeft(regressionVec, xCoeffsG);
yG.mulTransposeLeft(regressionVec, yCoeffsG);
int X = (int) Math.round(xG.getElement(0, 0));
int Y = (int) Math.round(yG.getElement(0, 0));
if (X >= 0 && Y >= 0 && X < imageWidth && Y < imageHeight) {
target.setSample(xi, yi, bi, outDataCube.getSampleFloat(X, Y, bi));
} else {
// TODO: Change with fillvalue
// target.setSample(xi, yi, bi, Float.NaN);
target.setSample(xi, yi, bi, fillValue);
}
}
}
}
return target;
}
Aggregations