Fixed previous bug, working snapshot

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

@@ -956,11 +956,12 @@ public class VegetationLMA {
 		}
 		ImageDtt.startAndJoin(threads);
 		// regularization weights and derivatives
-		
+		// splitting alpha_lpf from alpha_loss+alpha_push
 		int ind_next = y_vector.length;
-		if ((alpha_lpf >= 0) || (alpha_loss > 0)) {
+		
+		if ((alpha_loss > 0) || (alpha_push > 0)) {
 			int dbg_nx = -76340;
-			final int ind_y_alpha = ind_next;
+			final int ind_y_alpha_loss = ind_next;
 			ind_next += num_pars_vegetation_alpha;
 			ai.set(0);
 			for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -968,7 +969,7 @@ public class VegetationLMA {
 					public void run() {
 						for (int n = ai.getAndIncrement(); n < num_pars_vegetation_alpha; n = ai.getAndIncrement()) {
 							int np = ind_pars_vegetation_alpha + n; // index of the alpha parameter
-							int nx = n + ind_y_alpha; //  y_vector.length; // x - index
+							int nx = n + ind_y_alpha_loss; //  y_vector.length; // x - index
 							if (nx == dbg_nx) {
 								System.out.println("getFxDerivs(): n="+n+", nx="+nx);
 							}
@@ -1027,22 +1028,66 @@ public class VegetationLMA {
 									}
 								}								
 							}
+						}
+					}
+				};
+			}
+			ImageDtt.startAndJoin(threads);
+		} // if ((alpha_loss > 0) || (alpha_push > 0)){
+		
 		
 		
+		
+		if (alpha_lpf >= 0) {
+			int dbg_nx = -76340;
+			final int ind_y_alpha_lpf = ind_next;
+			ind_next += num_pars_vegetation_alpha;
+			ai.set(0);
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int n = ai.getAndIncrement(); n < num_pars_vegetation_alpha; n = ai.getAndIncrement()) {
+							int np = ind_pars_vegetation_alpha + n; // index of the alpha parameter
+							int nx = n + ind_y_alpha_lpf; //  y_vector.length; // x - index
+							if (nx == dbg_nx) {
+								System.out.println("getFxDerivs(): n="+n+", nx="+nx);
+							}
+							double d = 0;
+							fX[nx] = 0.0;
+							double avg = 0;
+							int nn = 0;
+							double neib_min = Double.POSITIVE_INFINITY, neib_max = Double.NEGATIVE_INFINITY;
+							for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
+								int di =  alpha_neibs[n][i];
+								d=0;
+								if (di >= 0) {
+									d = vector[di];  // d - full parameter index
+									avg+=d;
+									if (d < neib_min) neib_min = d; 
+									if (d > neib_max) neib_max = d; 
+									nn++;
+								} else if (di < -1) {
+									d = tvao[TVAO_VEGETATION_ALPHA][-di - 2]; 
+									avg+=d;
+									if (d < neib_min) neib_min = d; 
+									if (d > neib_max) neib_max = d; 
+									nn++;
+								}
+							}
+							avg /= nn; // average
 							// add cost for difference between this alpha and average of 4 neighbors (when they exist
 							// applies to alpha before cosine, so it will pull borders even when alpha<0 or alpha > 1 (zero derivatives)
 							//alpha_scale_avg
-							if (alpha_lpf > 0) { // should always be > 0 to provide stability for out-of-range alpha
 							double mm = neib_min + (neib_max-neib_min) * alpha_mm_hole;
 							double effective_alpha_lpf = alpha_lpf; 
 							if (!Double.isNaN(alpha_mm_hole) && (vector[np] <= mm) && ((neib_max-neib_min) >= alpha_diff_hole)) {
 								effective_alpha_lpf = 0.0; // disable alpha_lpf
 							}
 							// disable pull to average of neighbors for small holes in vegetation (local "almost minimum")
-///								if (Double.isNaN(alpha_mm_hole) || (vector[np] > mm) || ((neib_max-neib_min) < alpha_diff_hole)) {
+							///								if (Double.isNaN(alpha_mm_hole) || (vector[np] > mm) || ((neib_max-neib_min) < alpha_diff_hole)) {
 							//		this.alpha_mm_hole =     alpha_mm_hole;
 							if (alpha_scale_avg == 1) {
-									fX[nx] += effective_alpha_lpf * (vector[np] - avg);
+								fX[nx] = effective_alpha_lpf * (vector[np] - avg);
 								if (jt != null) {
 									jt[np][nx] += effective_alpha_lpf;
 									for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
@@ -1055,7 +1100,7 @@ public class VegetationLMA {
 							} else {
 								if ((avg > 0) && (avg < 1)) { 
 									double savg = 0.5 + (avg - 0.5) * alpha_scale_avg;
-										fX[nx] += effective_alpha_lpf * (vector[np] - savg);
+									fX[nx] = effective_alpha_lpf * (vector[np] - savg);
 									if (jt != null) {
 										//											jt[np][nx] += effective_alpha_lpf;
 										for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
@@ -1067,9 +1112,9 @@ public class VegetationLMA {
 									}
 								} else {
 									if (avg <= 0) {
-											fX[nx] += effective_alpha_lpf * (vector[np] + 0.5*(alpha_scale_avg - 1.0));
+										fX[nx] = effective_alpha_lpf * (vector[np] + 0.5*(alpha_scale_avg - 1.0));
 									} else { // avg > 1
-											fX[nx] += effective_alpha_lpf * (vector[np] - 0.5*(alpha_scale_avg + 1.0));
+										fX[nx] = effective_alpha_lpf * (vector[np] - 0.5*(alpha_scale_avg + 1.0));
 									}
 								}
 								if (jt != null) {
@@ -1078,7 +1123,6 @@ public class VegetationLMA {
 							}	
 						}
 					}
-					}
 				};
 			}
 			ImageDtt.startAndJoin(threads);
@@ -1178,231 +1222,6 @@ public class VegetationLMA {
 		return fX;		
 	}
 	
-	private double [] getFxDerivs_precos(
-			final double []   vector,
-			final double [][] jt, // should be null or initialized with [vector.length][]
-			final int         debug_level)
-	{
-		double [] fX = new double [weights.length]; // num_pairs + vector.length];
-		if (jt != null) {
-			for (int i = 0; i < jt.length; i++) {
-				jt[i] = new double [weights.length]; // weights.length];
-			}
-		}
-		final Thread[] threads = ImageDtt.newThreadArray();
-		final AtomicInteger ai = new AtomicInteger(0);
-		for (int ithread = 0; ithread < threads.length; ithread++) {
-			threads[ithread] = new Thread() {
-				public void run() {
-					double [] vegetation = new double [4]; 
-					double [] alpha =      new double [4];
-					for (int n = ai.getAndIncrement(); n < y_vector.length; n = ai.getAndIncrement()) {
-//						int nscene =   data_source[n][0][0];
-//						int indx =     data_source[n][0][1];
-						double terrain = vector[data_source[n][0][2]];
-						double [] cw = corners_weights[n]; 
-						double d; 
-						int [] indx_vegetation = data_source[n][1]; 
-						int [] indx_alpha =      data_source[n][2];
-						double sum_v =0, sum_a =0;
-						if (cw != null) {
-							for (int i = 0; i < 4; i++ ) {
-								int iv = indx_vegetation[i], ia=indx_alpha[i];
-								vegetation[i] = cw[i] * ((iv >= 0) ? vector[iv]: tvao[TVAO_VEGETATION][-1-iv]);
-								alpha[i] =      cw[i] * ((ia >= 0) ? vector[ia]: tvao[TVAO_VEGETATION_ALPHA][-1-ia]);
-								sum_v += vegetation[i];
-								sum_a += alpha[i];
-							}
-							d = terrain * (1.0 - sum_a) + sum_v * sum_a;
-							if (jt != null) {
-								jt[data_source[n][0][2]][n] = 1 - sum_a; // d/dterrain
-								for (int i = 0; i < 4; i++ ) {
-									if (indx_vegetation[i] >= 0) {
-//										jt[data_source[n][1][indx_vegetation[i]]][n] = cw[i] * sum_a; // d/dvegetation[i]
-										jt[data_source[n][1][i]][n] = cw[i] * sum_a; // d/dvegetation[i]
-									}
-									if (indx_alpha[i] >= 0) {
-//										jt[data_source[n][2][indx_alpha[i]]][n] = cw[i] * (sum_v - terrain); // d/dalpha[i]
-										jt[data_source[n][2][i]][n] = cw[i] * (sum_v - terrain); // d/dalpha[i]
-									}
-								}								
-								
-							}
-
-						} else {
-							d = terrain; 
-							if (jt != null) {
-								jt[data_source[n][0][2]][n] = 1;  // d/dterrain
-							}
-						}
-						double scene_offs = vector[data_source[n][0][3]];
-						fX[n] = d + scene_offs;
-						if (jt != null) {
-							jt[data_source[n][0][3]][n] = 1;
-						}
-					}
-				}
-			};
-		}
-		ImageDtt.startAndJoin(threads);
-		// regularization weights and derivatives
-		
-		int ind_next = y_vector.length;
-		if (alpha_lpf >= 0) {
-			final int ind_y_alpha = ind_next;
-			ind_next += num_pars_vegetation_alpha;
-			ai.set(0);
-			for (int ithread = 0; ithread < threads.length; ithread++) {
-				threads[ithread] = new Thread() {
-					public void run() {
-						for (int n = ai.getAndIncrement(); n < num_pars_vegetation_alpha; n = ai.getAndIncrement()) {
-							int np = ind_pars_vegetation_alpha + n; // index of the alpha parameter
-							int nx = n + ind_y_alpha; //  y_vector.length; // x - index
-							double d = 0;
-							fX[nx] = 0.0;
-							double alpha = vector[np];
-							if (alpha < alpha_offset) {
-								d = alpha- alpha_offset;
-							} else if (alpha > (1 - alpha_offset)) {
-								d = alpha - (1.0 - alpha_offset);
-							}
-							if (d != 0) {
-								fX[nx] = d * d * alpha_loss;
-								if (jt != null) {
-									jt[np][nx] = 2 * alpha_loss * d; // d/dalpha[i]		
-								}							
-							}
-							// add cost for difference between this alpha and average of 4 neighbors (when they exist
-							if (alpha_lpf > 0) {
-								double avg = 0;
-								int nn = 0;
-								for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
-									int di =  alpha_neibs[n][i];
-									d=0;
-									if (di >= 0) {
-										d = vector[di];  // d - full parameter index
-										avg+=d;
-										nn++;
-									} else if (di < -1) {
-										d = tvao[TVAO_VEGETATION_ALPHA][-di - 2]; 
-										avg+=d;
-										nn++;
-									}
-								}
-								avg /= nn; // average
-								fX[nx] += alpha_lpf * (vector[np] - avg);
-								if (jt != null) {
-									jt[np][nx] += alpha_lpf;
-									for (int i = 0; i < alpha_neibs[n].length; i++) { // now 4, may be increased
-										int di =  alpha_neibs[n][i];
-										if (di > 0) {
-											jt[di][nx] -= alpha_lpf/nn;
-										}
-									}								
-								}							
-							}
-						}
-					}
-				};
-			}
-			ImageDtt.startAndJoin(threads);
-		} // if (alpha_lpf >= 0) {
-		
-		if (terr_lpf >= 0) {
-			final int ind_y_terr = ind_next;
-			ind_next += num_pars_terrain;
-			ai.set(0);
-			for (int ithread = 0; ithread < threads.length; ithread++) {
-				threads[ithread] = new Thread() {
-					public void run() {
-						for (int n = ai.getAndIncrement(); n < num_pars_terrain; n = ai.getAndIncrement()) {
-							int np = ind_pars_terrain + n; // index of the alpha parameter
-							int nx = n + ind_y_terr; //  y_vector.length; // x - index
-							double d = 0;
-							if (terr_lpf > 0) {
-								double avg = 0;
-								int nn = 0;
-								for (int i = 0; i < terr_neibs[n].length; i++) { // now 4, may be increased
-									int di =  terr_neibs[n][i];
-									d=0;
-									if (di >= 0) {
-										d = vector[di];  // d - full parameter index
-										avg+=d;
-										nn++;
-									} else if (di < -1) {
-										d = tvao[TVAO_TERRAIN][-di - 2]; 
-										avg+=d;
-										nn++;
-									}
-								}
-								avg /= nn; // average
-								fX[nx] += terr_lpf * (vector[np] - avg);
-								if (jt != null) {
-									jt[np][nx] += terr_lpf;
-									for (int i = 0; i < terr_neibs[n].length; i++) { // now 4, may be increased
-										int di =  terr_neibs[n][i];
-										if (di > 0) {
-											jt[di][nx] -= terr_lpf/nn;
-										}
-									}								
-								}							
-							}
-						}
-					}
-				};
-			}
-			ImageDtt.startAndJoin(threads);
-		}
-		if (veget_lpf >= 0) {
-			final int ind_y_veget = ind_next;
-			ind_next += num_pars_vegetation;
-			ai.set(0);
-			
-			for (int ithread = 0; ithread < threads.length; ithread++) {
-				threads[ithread] = new Thread() {
-					public void run() {
-						for (int n = ai.getAndIncrement(); n < num_pars_vegetation; n = ai.getAndIncrement()) {
-							int np = ind_pars_vegetation + n; // index of the alpha parameter
-							int nx = n + ind_y_veget; //  y_vector.length; // x - index
-							double d = 0;
-							if (veget_lpf > 0) {
-								double avg = 0;
-								int nn = 0;
-								for (int i = 0; i < veget_neibs[n].length; i++) { // now 4, may be increased
-									int di =  veget_neibs[n][i];
-									d=0;
-									if (di >= 0) {
-										d = vector[di];  // d - full parameter index
-										avg+=d;
-										nn++;
-									} else if (di < -1) {
-										d = tvao[TVAO_VEGETATION][-di - 2]; 
-										avg+=d;
-										nn++;
-									}
-								}
-								avg /= nn; // average
-								fX[nx] += veget_lpf * (vector[np] - avg);
-								if (jt != null) {
-									jt[np][nx] += veget_lpf;
-									for (int i = 0; i < veget_neibs[n].length; i++) { // now 4, may be increased
-										int di =  veget_neibs[n][i];
-										if (di > 0) {
-											jt[di][nx] -= veget_lpf/nn;
-										}
-									}								
-								}							
-							}
-						}
-					}
-				};
-			}
-			ImageDtt.startAndJoin(threads);
-		}
-		return fX;		
-	}
-
-	
 	
 	private double [][] getFxDerivsDelta(
 			double []         vector,
@@ -1629,6 +1448,9 @@ public class VegetationLMA {
 //		int extra_samples = num_pars_vegetation_alpha; // in the future may be more regularization
 		int extra_samples = 0;
 		// using >=0 no use 0 as NOP but reserve space, <0 - do not reserve space
+		
+		//(alpha_loss > 0)
+		if ((alpha_loss > 0) || (alpha_push > 0)) extra_samples+= num_pars_vegetation_alpha; // need to split loss (always positive) from alpha_lpf
 		if (alpha_lpf >= 0) extra_samples+=	num_pars_vegetation_alpha;
 		if (terr_lpf >= 0)  extra_samples+=	num_pars_terrain;
 		if (veget_lpf >= 0) extra_samples+=	num_pars_vegetation;

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

@@ -677,7 +677,8 @@ public class VegetationModel {
 			String []          titles,
 			int                debugLevel) {
 		boolean diff_mode = true;
-		Rectangle woi50 =       new Rectangle(143,317,35,35);
+//		Rectangle woi50 =       new Rectangle(143,317,35,35);
+		Rectangle woi50 =       new Rectangle(160,317,35,35); // X+17
 		/* Bad, pull terr lower, so it all shifts to terrain
 		int       min_scenes =   10;
 		double    default_alpha = 0.8;
@@ -707,7 +708,7 @@ public class VegetationModel {
 		double    reg_weights =   0.25;        // fraction of the total weight used for regularization
 		double    alpha_loss =   100.0; // 10.0; /// 100.0; // 10.0;  // 10000.0; // 1000.0; // 100.; // 10.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
 		double    alpha_offset =  0.0; // 0.02; // 0.03;    // if >0,  start losses above 0.0 and below 1.0;
-		double    alpha_lpf =    10; /// 15; // 10.0; // 5.0; // 10.0; // 3; // 10;   // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
+		double    alpha_lpf =    5; /// 10; /// 15; // 10.0; // 5.0; // 10.0; // 3; // 10;   // 20; // 6.0; // 3.0; // 2.0;    // 1.5; // 5.0; // 0.5;         // pull to average of 4 neighbors
 		boolean   alpha_piece_linear = true; // false; // true;
 		double    alpha_scale_avg = 1.0; // 1.1; // 0.9; // 2.0; // 1.5; // scale average alpha (around 0.5) when pulling to it
 		double    alpha_push =       12; // 10.0; // 15.0;   // push from alpha==0.5