Cholesky done

src/main/java/com/elphel/imagej/common/CholeskyBlock.java

@@ -374,9 +374,9 @@ public class CholeskyBlock {
 			// Calculate A in one tile column of the remaining A2'
 			// start with diagonal (top) tile, and calculate its Cholesky
 			// In parallel, calculate A for all tiles in that column below diagonal
-			if (tile_diag == (n-1)) {
-				System.out.println("choleskyBlockMulti() last pass n= "+n+", tile_diag="+tile_diag);
-			}
+//			if (tile_diag == (n-1)) {
+//				System.out.println("choleskyBlockMulti() last pass n= "+n+", tile_diag="+tile_diag);
+//			}
 			ai.set(ftile_diag);
 			for (int ithread = 0; ithread < threads.length; ithread++) { // first sum for pairs
 				threads[ithread] = new Thread() {

src/main/java/com/elphel/imagej/common/CholeskyBlockTest.java

@@ -51,8 +51,9 @@ public class CholeskyBlockTest {
 		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());
+		Matrix wjtjlambda_copy2 = new Matrix(wjtjlambda.getArrayCopy());
 		int matrix_size = wjtjlambda.getRowDimension();
-		double [] starts = new double[4];
+		double [] starts = new double[6];
 		double start_time = (((double) System.nanoTime()) * 1E-9);
 		CholeskyDecomposition choleskyDecomposition = new CholeskyDecomposition(wjtjlambda_copy);
 		starts[0] = (((double) System.nanoTime()) * 1E-9) - start_time;
@@ -62,7 +63,10 @@ public class CholeskyBlockTest {
 		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;
-
+		CholeskyFloat choleskyFloat  = new CholeskyFloat(wjtjlambda_copy2);
+		starts[4] = (((double) System.nanoTime()) * 1E-9) - start_time;
+		Matrix mdelta_cholesky_float =  choleskyFloat.solve(jty);
+		starts[5] = (((double) System.nanoTime()) * 1E-9) - start_time;
 		
 		System.out.println("testCholesky(): matrix size=            "+matrix_size);
 		System.out.println("testCholesky(): block_size=             "+block_size);
@@ -70,9 +74,13 @@ public class CholeskyBlockTest {
 		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(): CholeskyFloat():        "+(starts[4]-starts[3])+" sec");
+		System.out.println("testCholesky(): float.solve():          "+(starts[5]-starts[4])+" 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(),
@@ -80,6 +88,15 @@ public class CholeskyBlockTest {
 				ch_diff.getRowPackedCopy()};
 		String[] dbg_titles = {"choleskyBlock","choleskyDecomposition",
 				"choleskyBlock-choleskyDecomposition"};
+				*/
+		Matrix ch_diff =      choleskyFloat.getL().minus(choleskyDecomposition.getL());
+		double [][] dbg_img = {
+				choleskyFloat.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,
@@ -87,7 +104,7 @@ public class CholeskyBlockTest {
 				true,
 				dbg_title, // String title)
 				dbg_titles);
-		return new Matrix[] {mdelta_cholesky, mdelta_cholesky_block};
+		return new Matrix[] {mdelta_cholesky, mdelta_cholesky_block,mdelta_cholesky_float};
 		/*
 		    First run, next MULTITHREADED ones are faster. Ran on an 8-core x86, 16 threads
 			choleskyBlockMulti() last pass n= 33, tile_diag=32

src/main/java/com/elphel/imagej/common/CholeskyFloat.java

+package com.elphel.imagej.common;
+
+import Jama.Matrix;
+
+// For testing for future GPU - is float enough
+public class CholeskyFloat {
+	   private float[][] L;
+	   private int n;
+	   
+	   public CholeskyFloat (Matrix mA) {
+		   double [][] M = mA.getArray();
+		   float [][] A = new float [M.length] [M[0].length];
+		   for (int i = 0; i < A.length; i++) {
+			   for (int j = 0; j < A[i].length; j++) {
+				   A[i][j] = (float) M[i][j];
+			   }
+		   }
+		   setupCholeskyFloat(A);
+	   }
+	  
+		public Matrix getL() {
+			   double [][] dL = new double [n] [n];
+			   for (int i = 0; i < n; i++) {
+				   for (int j = 0; j < n; j++) {
+					   dL[i][j] = L[i][j];
+				   }
+			   }
+			return new Matrix(dL,n,n);
+		}
+	   
+	   
+	   public CholeskyFloat (float[][] A) {
+		   setupCholeskyFloat(A);
+	   }
+	   
+	   public void setupCholeskyFloat (float[][] A) {
+		   // Initialize.
+		   n = A.length;
+		   L = new float[n][n];
+		   // Main loop.
+		   for (int j = 0; j < n; j++) {
+			   float[] Lrowj = L[j];
+			   float d = 0;
+			   for (int k = 0; k < j; k++) {
+				   float[] Lrowk = L[k];
+				   float s = 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;
+			   }
+			   d = A[j][j] - d;
+			   L[j][j] = (float) Math.sqrt(Math.max(d,0));
+			   for (int k = j+1; k < n; k++) {
+				   L[j][k] = 0;
+			   }
+		   }
+		   return;
+	   }
+
+
+	   public Matrix solve (Matrix B) {
+		   double [][] M = B.getArray();
+		   float [][] b = new float [M.length] [M[0].length];
+		   for (int i = 0; i < b.length; i++) {
+			   for (int j = 0; j < b[i].length; j++) {
+				   b[i][j] = (float) M[i][j];
+			   }
+		   }
+		   float  [][] x = solve (b);
+		   double [][] d = new double [M.length] [M[0].length];
+		   for (int i = 0; i < x.length; i++) {
+			   for (int j = 0; j < x[i].length; j++) {
+				   d[i][j] = x[i][j];
+			   }
+		   }
+		   return new Matrix(d,n,x[0].length);		   
+	   }
+	   
+	   public float[][] solve (float [][] b) {
+		   if (b.length != n) {
+			   throw new IllegalArgumentException("Matrix row dimensions must agree.");
+		   }
+		   // Copy right hand side.
+		   int nx = b[0].length;
+		   float [][] x = new float [n][];
+		   for (int i = 0; i < n; i++) {
+			   x[i] = b[i].clone();
+		   }
+
+		   // Solve L*Y = B;
+		   for (int k = 0; k < n; k++) {
+			   for (int j = 0; j < nx; j++) {
+				   for (int i = 0; i < k ; i++) {
+					   x[k][j] -= x[i][j]*L[k][i];
+				   }
+				   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++) {
+				   for (int i = k+1; i < n ; i++) {
+					   x[k][j] -= x[i][j]*L[i][k];
+				   }
+				   x[k][j] /= L[k][k];
+			   }
+		   }
+		   return x;
+	   }
+
+}

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -878,6 +878,7 @@ min_str_neib_fpn	0.35
 ///	public boolean   terr_only_pix =           true; // force per-pixel terrain elevation in terrain-only mode, overwrite fits_disable[TVAO_TERR_ELEV_PIX]
 	public int       terr_only_series =       -1;   // similar to terr_last_series but for terrain-only mode (<0 - length of terr_only_num_iters)
 	public int []    terr_only_num_iters =   {25};  // number of iterations
+	public int       terr_cholesky =           0;
 	
 // second run
 	public boolean [] terr_recalc_weights =    {false,false,true}; // recalculate weight depending on terrain visibility
@@ -2255,7 +2256,7 @@ min_str_neib_fpn	0.35
 ///		gd.addCheckbox    ("Adjust pixel elevation",terr_only_pix,               "Force per-pixel terrain elevation adjustment in terrain-only mode.");
 		gd.addNumericField("Terrain-only iterations",terr_only_series, 0,3,"",   "Last LMA series in terrain-only mode, -1 - to last available in the sequence below.");
 		gd.addStringField ("Maximal iterations", intsToString(terr_only_num_iters), 40, "Maximal number of LMA iterations per series, 1 or several values.");
-		
+		gd.addChoice      ("LMA matrix solution mode:", VegetationLMA.CHOLESKY_NAMES, VegetationLMA.CHOLESKY_NAMES[terr_cholesky]);
 		gd.addMessage     ("Second LMA run");
 		gd.addStringField ("Recalculate weighths", booleansToString(terr_recalc_weights,2), 40,
 				"Recalculate weights depending on terrain visibility. Use comma/space separated list of true/false, 1/0 or +/-." );
@@ -3053,8 +3054,8 @@ min_str_neib_fpn	0.35
 ///		terr_only_pix =            gd.getNextBoolean();// boolean
 		terr_only_series =   (int) gd.getNextNumber();// int    
 		terr_only_num_iters =      StringToInts(gd.getNextString());
-		
-		terr_recalc_weights = StringToBooleans(gd.getNextString());// booleans
+		terr_cholesky =            gd.getNextChoiceIndex();
+		terr_recalc_weights =      StringToBooleans(gd.getNextString());// booleans
 		terr_recalc_opaque =       gd.getNextNumber(); // double 
 		terr_recalc_pedestal =     gd.getNextNumber(); // double 
 		terr_recalc_frac =         gd.getNextNumber(); // double 
@@ -3827,7 +3828,7 @@ min_str_neib_fpn	0.35
 ///		properties.setProperty(prefix+"terr_only_pix",               terr_only_pix+"");     // boolean
 		properties.setProperty(prefix+"terr_only_series",            terr_only_series+"");       // int    
 		properties.setProperty(prefix+"terr_only_num_iters",         intsToString(terr_only_num_iters)+"");  // int []
-		
+		properties.setProperty(prefix+"terr_cholesky",               terr_cholesky+"");           // int
 		properties.setProperty(prefix+"terr_recalc_weights",         booleansToString(terr_recalc_weights,2));     // boolean
 		properties.setProperty(prefix+"terr_recalc_opaque",          terr_recalc_opaque+"");      // double 
 		properties.setProperty(prefix+"terr_recalc_pedestal",        terr_recalc_pedestal+"");    // double 
@@ -4615,6 +4616,7 @@ min_str_neib_fpn	0.35
 ///		if (properties.getProperty(prefix+"terr_only_pix")!=           null) terr_only_pix=Boolean.parseBoolean(properties.getProperty(prefix+"terr_only_pix"));
 		if (properties.getProperty(prefix+"terr_only_series")!=        null) terr_only_series=Integer.parseInt(properties.getProperty(prefix+"terr_only_series"));
 		if (properties.getProperty(prefix+"terr_only_num_iters")!=     null) terr_only_num_iters=StringToInts((String) properties.getProperty(prefix+"terr_only_num_iters"));
+		if (properties.getProperty(prefix+"terr_cholesky")!=           null) terr_cholesky=Integer.parseInt(properties.getProperty(prefix+"terr_cholesky"));
 		
 		if (properties.getProperty(prefix+"terr_recalc_weights")!=     null) {
 			terr_recalc_weights= StringToBooleans(properties.getProperty(prefix+"terr_recalc_weights"));// booleans
@@ -5367,7 +5369,7 @@ min_str_neib_fpn	0.35
 ///		imp.terr_only_pix =                      this.terr_only_pix;                           
 		imp.terr_only_series =                   this.terr_only_series;                           
 		imp.terr_only_num_iters =                this.terr_only_num_iters.clone();
-
+		imp.terr_cholesky =                      this.terr_cholesky;
 		imp.terr_recalc_weights =                this.terr_recalc_weights.clone();                           
 		imp.terr_recalc_opaque =                 this.terr_recalc_opaque;                           
 		imp.terr_recalc_pedestal =               this.terr_recalc_pedestal;

src/main/java/com/elphel/imagej/vegetation/VegetationLMA.java

@@ -12,6 +12,9 @@ import java.util.HashSet;
 import java.util.concurrent.atomic.AtomicBoolean;
 import java.util.concurrent.atomic.AtomicInteger;
 
+import com.elphel.imagej.common.CholeskyBlock;
+import com.elphel.imagej.common.CholeskyBlockTest;
+import com.elphel.imagej.common.CholeskyFloat;
 import com.elphel.imagej.common.CholeskyLDLTMulti;
 import com.elphel.imagej.common.DoubleGaussianBlur;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
@@ -80,6 +83,14 @@ public class VegetationLMA {
 	public static final int [] SAVE_TYPES =
 		{TVAO_TERRAIN, TVAO_VEGETATION, TVAO_ALPHA, TVAO_ELEVATION,	TVAO_TERR_ELEV_PIX, TVAO_SCENE_OFFSET};
 
+	public static final int CHOLESKY_BLOCK =         0; // multithreaded
+	public static final int CHOLESKY_JAMA =          1; // standard, from JAMA
+	public static final int CHOLESKY_FLOAT =         2; // floating point precision (for GPU eveluation
+	public static final int CHOLESKY_LU =            3; // LU instead of Cholesky (old way)
+	public final static String[] CHOLESKY_NAMES = {"multi-threaded", "jama", "float","LU"};
+	
+	
+	
 	private final Rectangle   full;
 	private final int         image_length;
 	// should be non-null everywhere where elevation is defined?
@@ -289,6 +300,7 @@ public class VegetationLMA {
 	public boolean            debug_source_vegetation = false;
 	public int [][][]         debug_source_veget = null;
 	
+	public int                cholesky_mode =     0; // 0 - cholesky_block; 1 - cholesky_standard; 2 - cholesky_float; 3  - old 
 	
 	/**
 	 * What WOI to apply to fX items
@@ -364,7 +376,9 @@ public class VegetationLMA {
 			double    alpha_min,
 			double    alpha_max,
 			double    alpha_sigma,
+			int       cholesky_mode,
 			int debugLevel) {
+		this.cholesky_mode = cholesky_mode;
 //		this.alpha_initial_contrast = alpha_initial_contrast;
 		full = vegetationModel.full;
 		image_length = full.width * full.height;
@@ -923,7 +937,8 @@ public class VegetationLMA {
 ///			final boolean            terr_only_special,//true; // special sequences for terrain-only tiles
 ///			final boolean            terr_only_pix,    //true; // force per-pixel terrain elevation in terrain-only mode, overwrite fits_disable[TVAO_TERR_ELEV_PIX]
 			final double             boost_parallax, // increase weight of scene with maximal parallax relative to the reference scene
-			final double             max_parallax,  // do not consider maximal parallax above this (consider it a glitch) 
+			final double             max_parallax,  // do not consider maximal parallax above this (consider it a glitch)
+			
 			final int                debugLevel,
 			final String             debug_path,
 			final boolean            debug_save_improved, // Save debug image after successful LMA step."); 
@@ -1856,8 +1871,70 @@ public class VegetationLMA {
 		}
 		Matrix jty = (new Matrix(last_jt_decimated)).times(y_minus_fx_weighted);
 		Matrix mdelta = null; //  jtjl_inv.times(jty);
-		boolean use_cholesky = true; // false;
+//		boolean use_cholesky = true; // false;
 		double matrix_start_time = ((double) System.nanoTime()) * 1E-9; 
+		switch (cholesky_mode) { // Cholesky does not copy array
+		case CHOLESKY_BLOCK:
+			try {
+				mdelta =  (new CholeskyBlock(wjtjlambda)).solve(jty); // does not modify arrays (both matrices)
+			} catch (RuntimeException e) {
+				rslt[1] = true;
+				if (debug_level > -2) {
+					System.out.println("CHOLESKY_BLOCK: Singular Matrix!");
+				}
+				return rslt;
+			}
+			break;
+		case CHOLESKY_JAMA:
+			try {
+				mdelta =  (new CholeskyDecomposition(new Matrix(wjtjlambda.getArrayCopy()))).solve(jty);
+			} catch (RuntimeException e) {
+				rslt[1] = true;
+				if (debug_level > -2) {
+					System.out.println("CHOLESKY_JAMA: Singular Matrix!");
+				}
+				return rslt;
+			}
+			break;
+		case CHOLESKY_FLOAT:
+			try {
+				mdelta =  (new CholeskyFloat(wjtjlambda)).solve(jty); // does not modify arrays (both matrices)
+			} catch (RuntimeException e) {
+				rslt[1] = true;
+				if (debug_level > -2) {
+					System.out.println("CHOLESKY_BLOCK: Singular Matrix!");
+				}
+				return rslt;
+			}
+			break;
+		case CHOLESKY_LU: // old way - inverse() using LU
+			Matrix jtjl_inv = null;
+			try {
+				jtjl_inv = wjtjlambda.inverse(); // check for errors
+			} catch (RuntimeException e) {
+				rslt[1] = true;
+				if (debug_level > -2) {
+					System.out.println("Singular Matrix!");
+				}
+
+				return rslt;
+			}
+			if (debug_level>4) {
+				System.out.println("(JtJ + lambda*diag(JtJ).inv()");
+				jtjl_inv.print(18, 6);
+			}
+			//last_jt has NaNs
+			//		Matrix jty = (new Matrix(last_jt_decimated)).times(y_minus_fx_weighted);
+			if (debug_level>2) {
+				System.out.println("Jt * (y-fx)");
+				jty.print(18, 6);
+			}
+			mdelta = jtjl_inv.times(jty);
+			break;
+		
+		}
+		
+		/*
 		if (use_cholesky) {
 			try {
 				mdelta =  (new CholeskyDecomposition(wjtjlambda)).solve(jty);
@@ -1892,12 +1969,12 @@ public class VegetationLMA {
 			}
 			mdelta = jtjl_inv.times(jty);
 		}
-		
+		*/
 		
 		
 		if (debug_level>-2) {
 			System.out.println("lmaStep(): Matrix inverted in "+((((double) System.nanoTime()) * 1E-9)-matrix_start_time)+
-					", used "+(use_cholesky? "Cholesky":"LU"));
+					", used matrix mode "+CHOLESKY_NAMES[cholesky_mode]+" ("+cholesky_mode+")");
 		}
 
 //		
@@ -1907,8 +1984,8 @@ public class VegetationLMA {
 		Matrix mdelta_cholesky = choleskyDecomposition.solve(jty);
 		Matrix mdelta_cholesky_multi = choleskyLDLTMulti.solve(jty);
 		 */
-		boolean save_cholesky = true;
-		boolean debug_cholesky = true;
+		boolean save_cholesky = false; // true;
+		boolean debug_cholesky = false; // true;
 		if (save_cholesky) {
 			int dbg_n = wjtjlambda.getRowDimension();
 			double [] dbg_a0 = wjtjlambda.getRowPackedCopy();
@@ -1929,9 +2006,15 @@ public class VegetationLMA {
 					"spd_a_"+debug_iter); // String title)
 		}
 		if (debug_cholesky) {
-			CholeskyLDLTMulti.testCholesky(
+//			CholeskyLDLTMulti.testCholesky(
+//					wjtjlambda, // Matrix wjtjlambda,
+//					jty);
+			CholeskyBlockTest.testCholesky(
 					wjtjlambda, // Matrix wjtjlambda,
-					jty);
+					jty,
+					"spd_a_"+debug_iter);
+
+			
 		}
 		if (debug_level>2) {
 			System.out.println("mdelta");

src/main/java/com/elphel/imagej/vegetation/VegetationModel.java

@@ -1607,6 +1607,8 @@ public class VegetationModel {
 		int       terr_only_series =  clt_parameters.imp.terr_only_series;         //  -1;   // similar to terr_last_series but for terrain-only mode (<0 - length of terr_only_num_iters)
 		int []    terr_only_num_iters=clt_parameters.imp.terr_only_num_iters;      // {25};  // number of iterations
 		if (terr_only_series < 0) terr_only_series = terr_only_num_iters.length - 1; 
+		int       terr_cholesky=      clt_parameters.imp.terr_cholesky;            // LMA matrix inversion mode 
+		
 		
 		// Combine  parameters
 		int       border_width =      clt_parameters.imp.terr_border_width;        // 6;
@@ -2037,8 +2039,9 @@ public class VegetationModel {
 				alpha_init_offs,  // 0.01; // double    alpha_offs,
 				alpha_init_min,   // 0.7;  // double    alpha_min,
 				alpha_init_max,   // 0.9;  // double    alpha_max,
-				alpha_sigma, // 8.0;  // double    alpha_sigma,
-				debugLevel); // int debugLevel) 
+				alpha_sigma,      // 8.0;  // double    alpha_sigma,
+				terr_cholesky,    //int       cholesky_mode,
+				debugLevel);      // int debugLevel) 
 		
 		
 		if (debugLevel > -2) {