Skip to content

I
imagej-elphel
Commit 09a0cb0e authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

Moved tests to CholeskyBlockTest.java

parent b09025fb
Hide whitespace changes
Inline Side-by-side
Showing with 195 additions and 37 deletions
src/main/java/com/elphel/imagej/common/CholeskyBlock.java
View file @ 09a0cb0e
package com.elphel.imagej.common;
/**
** CholeskyBlock - multithreaded Chelesky decomposition
** CholeskyBlock - multithreaded Cholesky decomposition and solution
**
** Copyright (C) 2025 Elphel, Inc.
**
......@@ -30,16 +30,29 @@ public class CholeskyBlock {
public static int dflt_m = 70; // tile size if not specified during initialization
public static int dflt_threads = 100; // maximal number of threads to use for decomposition
public static int dflt_solve_threads =16; // 4; // maximal number of threads to use for solve()
public static boolean debug = false; // true; // false;
public final int m; // tile size
public final int np; // number of elements in row/col
public final int n; // number of tiles in row/col
public final int nf; // number of full rows/cols
public int decomp_threads;
public final double [] A;
public final double [] L;
public final double [] B;
public static boolean debug = false; // true;
private final int m; // tile size
private final int np; // number of elements in row/col
private final int n; // number of tiles in row/col
private final int nf; // number of full rows/cols
private int decomp_threads;
private final double [] A;
private final double [] L;
private final double [] B;
//wjtjlambda_copy1.getArray()
public CholeskyBlock (Matrix mA) {
m = dflt_m;
decomp_threads = dflt_threads;
np = mA.getRowDimension();
nf = np / m;
n = ((nf * m) < np)? (nf + 1) : nf;
A = new double [np*np]; // [n*n*m*m];
L = new double [A.length];
B = new double [np];
setup_ATriangle(mA.getArray());
choleskyBlockMulti();
}
public CholeskyBlock (
double [][] A_in) {
m = dflt_m;
......@@ -152,7 +165,7 @@ public class CholeskyBlock {
return new Matrix (B, np);
}
public double [][] get_LTriangle() {
private double [][] get_LTriangle() {
double [][] L_out = new double[np][np];
for (int tile_row = 0; tile_row < nf; tile_row++) {
for (int tile_col = 0; tile_col < tile_row; tile_col++) {
......@@ -213,11 +226,11 @@ public class CholeskyBlock {
* @param j tile column
* @return index in A and L arrays
*/
public int indx_IJ(int i, int j) {
private int indx_IJ(int i, int j) {
return j * (m * np) + ((j >= nf) ? (np-nf*m): m) * m * i;
}
public void setL21(
private void setL21(
int i, // i > j,
int j) { // j <nf
int indx_diag = indx_IJ(j,j);
......@@ -238,7 +251,7 @@ public class CholeskyBlock {
return;
}
public void setA22(
private void setA22(
int diag,// < col, < row
int row, // >= col
int col) {
......@@ -266,7 +279,7 @@ public class CholeskyBlock {
return;
}
public void choleskyBlockMulti() {
private void choleskyBlockMulti() {
final Thread[] threads = ImageDtt.newThreadArray(decomp_threads); // 1);
final AtomicInteger ai = new AtomicInteger(0);
cholesky_single(0);
......@@ -350,7 +363,7 @@ public class CholeskyBlock {
// Single-threaded, used for a single-tile Cholesky deconstruction
public void cholesky_single(int diag) {
private void cholesky_single(int diag) {
int h = (diag < nf) ? m : (np-nf*m);
int indx = indx_IJ(diag, diag);
Arrays.fill(L, indx, indx+h*h-1, 0);
......@@ -378,7 +391,7 @@ public class CholeskyBlock {
* B has data, will be modified
* @param diag number of the diagonal tile (last may be smaller)
*/
public void solveY(int diag) {
private void solveY(int diag) {
int h = (diag < nf) ? m : (np-nf*m);
int lindx = indx_IJ(diag, diag);
int xindx = diag * m; // start of B/x
......@@ -399,7 +412,7 @@ public class CholeskyBlock {
* @param diag number of the diagonal tile from the bottom-right,
* last may be smaller.
*/
public void solveX(int diag) {
private void solveX(int diag) {
int h = (diag < nf) ? m : (np-nf*m);
int lindx = indx_IJ(diag, diag);
int xindx = diag * m; // start of B/x
......@@ -420,7 +433,7 @@ public class CholeskyBlock {
* @param i tile row to subtract from in B
* @param j tile row corresponding to diagonal L tiles
*/
public void subYCol(
private void subYCol(
int i, // i > j,
int j) { // j < i
int lindx = indx_IJ(i,j);
......@@ -444,7 +457,7 @@ public class CholeskyBlock {
* @param j tile row corresponding to diagonal L tiles, from the bottom
*/
public void subXCol(
private void subXCol(
int i, // i < j,
int j) { // j > i
int lindx = indx_IJ(j,i); // transposed
......@@ -601,16 +614,17 @@ public class CholeskyBlock {
ImageDtt.startAndJoin(threads);
starts[4] = (((double) System.nanoTime()) * 1E-9) - start_time;
if (debug) {
System.out.println("\ncholeskyBlock.solve1(): ==== number of threads = "+threads.length+ "====");
System.out.println(String.format("choleskyBlock.solve1(): solveY(0) %12.9f sec",(starts[0])));
System.out.println(String.format("choleskyBlock.solve1(): solveY_other %12.9f sec- including 1-st column",(starts[2] - starts[0])));
System.out.println(String.format("choleskyBlock.solve1(): solveX(n-1) %12.9f sec",(starts[3] - starts[2])));
System.out.println(String.format("choleskyBlock.solve1(): solveX_other %12.9f sec",(starts[4] - starts[3])));
System.out.println("\ncholeskyBlock.solve(): ==== number of threads = "+threads.length+ "====");
System.out.println(String.format("choleskyBlock.solve(): solveY(0) %12.9f sec",(starts[0])));
System.out.println(String.format("choleskyBlock.solve(): solveY_other %12.9f sec- including 1-st column",(starts[2] - starts[0])));
System.out.println(String.format("choleskyBlock.solve(): solveX(n-1) %12.9f sec",(starts[3] - starts[2])));
System.out.println(String.format("choleskyBlock.solve(): solveX_other %12.9f sec",(starts[4] - starts[3])));
}
return;
}
public static Matrix solve (Matrix B, double [][] L) {
// for comparison/verification
public static Matrix solve_single (Matrix B, double [][] L) {
int n = L.length;
if (B.getRowDimension() != n) {
throw new IllegalArgumentException("Matrix row dimensions must agree.");
......@@ -633,5 +647,29 @@ public class CholeskyBlock {
}
return new Matrix(x,n);
}
/* Results
First run, next MULTITHREADED ones are faster. Ran on an 8-core x86, 16 threads
choleskyBlockMulti() last pass n= 33, tile_diag=32
testCholesky(): matrix size= 2258
testCholesky(): block_size= 70
testCholesky(): choleskyDecomposition: 1.7590559809468687 sec Jama CholeskyDecomposition, single-threaded
testCholesky(): CholeskyBlock(): 0.34954120498150587 sec FIRST run is faster than single-threaded, but slower than next
testCholesky(): cholesky.solve(): 0.030882437014952302 sec Jama CholeskyDecomposition.solve()
testCholesky(): block.solve(): 0.051133116940036416 sec FIRST run is slower than single-threaded!
testCholesky(): title= spd_a_2.tif
testCholesky(): dbg_title= spd_a_2.tifch_diff_choleskyBlock-choleskyDecomposition-70
DEBUG_LEVEL = 1, CLT_PARAMETERS.lwir.getDebugLevel() = 0 LOG_LEVEL=ERROR LOG_LEVEL_SET=false
choleskyBlockMulti() last pass n= 33, tile_diag=32
testCholesky(): matrix size= 2258
testCholesky(): block_size= 70
testCholesky(): choleskyDecomposition: 1.7251908951438963 sec
testCholesky(): CholeskyBlock(): 0.19312086096033454 sec
testCholesky(): cholesky.solve(): 0.03081041993573308 sec
testCholesky(): block.solve(): 0.005443983944132924 sec
testCholesky(): title= spd_a_2.tif
testCholesky(): dbg_title= spd_a_2.tifch_diff_choleskyBlock-choleskyDecomposition-70
*/
}
src/main/java/com/elphel/imagej/common/CholeskyBlockTest.java 0 → 100644
View file @ 09a0cb0e
package com.elphel.imagej.common;
import Jama.CholeskyDecomposition;
import Jama.Matrix;
import ij.ImagePlus;
public class CholeskyBlockTest {
public static void testCholesky(ImagePlus imp_src) {
float [] fpixels = (float[]) imp_src.getProcessor().getPixels();
int width = imp_src.getWidth();
int height = imp_src.getHeight();
int n = height;
double [][] a = new double [n][n];
double [][] b = new double [n][1];
for (int i = 0; i < n; i++) {
for (int j= 0; j < n; j++) {
a[i][j] = fpixels[i*width+j];
}
b[i][0] = fpixels[i*width+n];
}
testCholesky(
new Matrix(a), // Matrix wjtjlambda,
new Matrix(b), // Matrix jty)
imp_src.getTitle()); // String title);
}
public static Matrix[] testCholesky(
Matrix wjtjlambda_in,
Matrix jty_in,
String title) {
int block_size = 70; // 100;//70~best // 64; // 60; // 70; // 80; // 120; // 150; // 200; // 100; // 4; // 10; // 100; // 28;
boolean truncate = false;
int trunc_size = 199; // 0 to use full size
Matrix wjtjlambda,jty;
if (truncate) {
int tr_size = (trunc_size==0)?(((wjtjlambda_in.getRowDimension())/block_size) * block_size):trunc_size;
wjtjlambda=wjtjlambda_in.getMatrix (
0, // int i0,
tr_size-1, // int i1,
0, // int j0,
tr_size-1); //int j1)
jty=jty_in.getMatrix (
0, // int i0,
tr_size-1, // int i1,
0, // int j0,
0); //int j1)
} else {
wjtjlambda=wjtjlambda_in;
jty = jty_in;
}
String dbg_title=title+"ch_diff_choleskyBlock-choleskyDecomposition-"+block_size+(truncate?("-truncate"+trunc_size):"");
Matrix wjtjlambda_copy1 = new Matrix(wjtjlambda.getArrayCopy());
Matrix wjtjlambda_copy = new Matrix(wjtjlambda.getArrayCopy());
int matrix_size = wjtjlambda.getRowDimension();
double [] starts = new double[4];
double start_time = (((double) System.nanoTime()) * 1E-9);
CholeskyDecomposition choleskyDecomposition = new CholeskyDecomposition(wjtjlambda_copy);
starts[0] = (((double) System.nanoTime()) * 1E-9) - start_time;
CholeskyBlock choleskyBlock = new CholeskyBlock(wjtjlambda_copy1.getArray(),block_size);
starts[1] = (((double) System.nanoTime()) * 1E-9) - start_time;
Matrix mdelta_cholesky = choleskyDecomposition.solve(jty);
starts[2] = (((double) System.nanoTime()) * 1E-9) - start_time;
Matrix mdelta_cholesky_block = choleskyBlock.solve(jty);
starts[3] = (((double) System.nanoTime()) * 1E-9) - start_time;
System.out.println("testCholesky(): matrix size= "+matrix_size);
System.out.println("testCholesky(): block_size= "+block_size);
System.out.println("testCholesky(): choleskyDecomposition: "+(starts[0])+" sec");
System.out.println("testCholesky(): CholeskyBlock(): "+(starts[1]-starts[0])+" sec");
System.out.println("testCholesky(): cholesky.solve(): "+(starts[2]-starts[1])+" sec");
System.out.println("testCholesky(): block.solve(): "+(starts[3]-starts[2])+" sec");
System.out.println("testCholesky(): title= "+title);
System.out.println("testCholesky(): dbg_title= "+dbg_title);
Matrix ch_diff = choleskyBlock.getL().minus(choleskyDecomposition.getL());
double [][] dbg_img = {
choleskyBlock.getL().getRowPackedCopy(),
choleskyDecomposition.getL().getRowPackedCopy(),
ch_diff.getRowPackedCopy()};
String[] dbg_titles = {"choleskyBlock","choleskyDecomposition",
"choleskyBlock-choleskyDecomposition"};
ShowDoubleFloatArrays.showArrays(
dbg_img, // double[] pixels,
ch_diff.getRowDimension(), // int width,
ch_diff.getRowDimension(), // int height,
true,
dbg_title, // String title)
dbg_titles);
return new Matrix[] {mdelta_cholesky, mdelta_cholesky_block};
/*
First run, next MULTITHREADED ones are faster. Ran on an 8-core x86, 16 threads
choleskyBlockMulti() last pass n= 33, tile_diag=32
testCholesky(): matrix size= 2258
testCholesky(): block_size= 70
testCholesky(): choleskyDecomposition: 1.7590559809468687 sec Jama CholeskyDecomposition, single-threaded
testCholesky(): CholeskyBlock(): 0.34954120498150587 sec FIRST run is faster (5x) than single-threaded, but slower than next
testCholesky(): cholesky.solve(): 0.030882437014952302 sec Jama CholeskyDecomposition.solve()
testCholesky(): block.solve(): 0.051133116940036416 sec FIRST run is slower than single-threaded!
testCholesky(): title= spd_a_2.tif
testCholesky(): dbg_title= spd_a_2.tifch_diff_choleskyBlock-choleskyDecomposition-70
DEBUG_LEVEL = 1, CLT_PARAMETERS.lwir.getDebugLevel() = 0 LOG_LEVEL=ERROR LOG_LEVEL_SET=false
choleskyBlockMulti() last pass n= 33, tile_diag=32
testCholesky(): matrix size= 2258
testCholesky(): block_size= 70
testCholesky(): choleskyDecomposition: 1.7251908951438963 sec
testCholesky(): CholeskyBlock(): 0.19312086096033454 sec Second run, 9x acceleration
testCholesky(): cholesky.solve(): 0.03081041993573308 sec
testCholesky(): block.solve(): 0.005443983944132924 sec Second run, 5.7x acceleration
testCholesky(): title= spd_a_2.tif
testCholesky(): dbg_title= spd_a_2.tifch_diff_choleskyBlock-choleskyDecomposition-70
*/
}
}
src/main/java/com/elphel/imagej/common/CholeskyLLTMulti.java
View file @ 09a0cb0e
......@@ -304,7 +304,7 @@ public class CholeskyLLTMulti {
// starts[2] = (((double) System.nanoTime()) * 1E-9) - start_time;
// choleskyBlock.choleskyBlockMulti();
starts[3] = (((double) System.nanoTime()) * 1E-9) - start_time;
double [][] LTriangle = choleskyBlock.get_LTriangle();
double [][] LTriangle = choleskyBlock.getL().getArray(); // get_LTriangle();
Matrix mdelta_cholesky = choleskyDecomposition.solve(jty);
starts[4] = (((double) System.nanoTime()) * 1E-9) - start_time;
Matrix mdelta_cholesky_single = CholeskyLLTMulti.solve_single(jty, choleskyLLTMulti_single.getL());
......
src/main/java/com/elphel/imagej/correction/Eyesis_Correction.java
View file @ 09a0cb0e
......@@ -79,9 +79,9 @@ import javax.imageio.metadata.IIOInvalidTreeException;
import javax.imageio.metadata.IIOMetadata;
import javax.imageio.metadata.IIOMetadataFormatImpl;
import javax.imageio.metadata.IIOMetadataNode;
import javax.imageio.plugins.tiff.ExifGPSTagSet;
//import javax.imageio.plugins.tiff.ExifGPSTagSet;
import javax.imageio.plugins.tiff.TIFFDirectory;
import javax.imageio.plugins.tiff.TIFFTagSet;
//import javax.imageio.plugins.tiff.TIFFTagSet;
import javax.swing.JFileChooser;
import javax.swing.filechooser.FileFilter;
......@@ -96,8 +96,9 @@ import com.elphel.imagej.cameras.CLTParameters;
import com.elphel.imagej.cameras.ColorProcParameters;
import com.elphel.imagej.cameras.EyesisCorrectionParameters;
import com.elphel.imagej.cameras.ThermalColor;
import com.elphel.imagej.common.CholeskyBlockTest;
import com.elphel.imagej.common.CholeskyLDLTMulti;
import com.elphel.imagej.common.CholeskyLLTMulti;
//import com.elphel.imagej.common.CholeskyLLTMulti;
import com.elphel.imagej.common.DoubleFHT;
import com.elphel.imagej.common.DoubleGaussianBlur;
import com.elphel.imagej.common.GenericJTabbedDialog;
......@@ -108,7 +109,7 @@ import com.elphel.imagej.gpu.GPUTileProcessor;
import com.elphel.imagej.gpu.GpuQuad;
import com.elphel.imagej.gpu.JCuda_ImageJ_Example_Plugin;
import com.elphel.imagej.ims.EventLogger;
import com.elphel.imagej.ims.Imx5;
//import com.elphel.imagej.ims.Imx5;
import com.elphel.imagej.jp4.JP46_Reader_camera;
import com.elphel.imagej.lwir.LwirReader;
import com.elphel.imagej.orthomosaic.OrthoMap;
......@@ -131,7 +132,7 @@ import com.elphel.imagej.tileprocessor.SymmVector;
import com.elphel.imagej.tileprocessor.TwoQuadCLT;
import com.elphel.imagej.tileprocessor.lwoc.LwirWorld;
import com.elphel.imagej.vegetation.VegetationModel;
import com.elphel.imagej.vegetation.VegetationSegment;
//import com.elphel.imagej.vegetation.VegetationSegment;
import ij.CompositeImage;
import ij.IJ;
......@@ -867,7 +868,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
addButton("Vegetation LMA", panelOrange, color_process);
addButton("Combine LMA Segments", panelOrange, color_process);
addButton("Generate LWIR target", panelOrange, color_process);
addButton("Test LDLT Cholesky", panelOrange, color_process);
// addButton("Test LDLT Cholesky", panelOrange, color_process);
addButton("Test LLT Cholesky", panelOrange, color_process);
plugInFrame.add(panelOrange);
}
......@@ -5810,21 +5811,23 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
VegetationModel.processVegetation(
SYNC_COMMAND, // // SyncCommand SYNC_COMMAND,
CLT_PARAMETERS, //CLTParameters clt_parameters,
true); //boolean combine_segments);
true); //boolean combine_segments);
/*
} else if (label.equals("Test LDLT Cholesky")) {
ImagePlus imp_sel = WindowManager.getCurrentImage();
if (imp_sel == null) {
IJ.showMessage("Error", "No images selected");
return;
}
CholeskyLDLTMulti.testCholesky(imp_sel); // ImagePlus imp.);
CholeskyLDLTMulti.testCholesky(imp_sel); // ImagePlus imp.);
*/
} else if (label.equals("Test LLT Cholesky")) {
ImagePlus imp_sel = WindowManager.getCurrentImage();
if (imp_sel == null) {
IJ.showMessage("Error", "No images selected");
return;
}
CholeskyLLTMulti.testCholesky(imp_sel); // ImagePlus imp.);
CholeskyBlockTest.testCholesky(imp_sel); // ImagePlus imp.);
}
//
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment