regularization to alpha

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

@@ -37,6 +37,11 @@ public class VegetationLMA {
 //	private double            default_alpha;
 	private int               num_pars_terrain; 
 	private int               num_pars_vegetation; 
+	private int               num_pars_vegetation_alpha; //== num_pars_vegetation
+	private int               ind_pars_terrain =   0;    // index of the first terrain parameter
+	private int               ind_pars_vegetation;       // index of the first vegetation parameter
+	private int               ind_pars_vegetation_alpha; // index of the first alpha parameter
+	private int               ind_pars_scenes;           // index of the first scene parameter
 	private int               num_pars_scenes;
 	private int [][][]        data_source; // [index][windx][3]:
 	                                       // 0:{ scene, full_indx, terr_indx, scale_indx} - always present as parameters
@@ -48,10 +53,13 @@ public class VegetationLMA {
 	private double []         y_vector = null;
 	private double            weight_pure =     0;
 	private double []         weights; // normalized so sum is 1.0 for all - samples and extra regularization
+	public double             alpha_loss;
+	public double             alpha_offset =    0; // if >0,  start losses above 0.0 and below 1.0;
 	public double             delta=            1e-5; // 7;
 	
-	
-	
+	public int                debug_index;
+	public double [][]        debug_image;
+	public static double []   debug_alpha_scales = {100,150}; 
 	
 	public int []             indices;
 	public int [][]           cpairs =          null;
@@ -101,8 +109,13 @@ public class VegetationLMA {
 			final int                min_scenes, // minimal number of scenes (inside woi) vegetation pixel must influence
 			final double             default_alpha,
 			final double []          scene_weights, // null or per-scene weights (higher for larger offsets)?
+			final double             reg_weights,        // fraction of the total weight used for regularization
+			final double             alpha_loss,         // quadratic loss when alpha reaches -1.0 or 2.0 
+			final double             alpha_offset,         // quadratic loss when alpha reaches -1.0 or 2.0			
 			final int                debugLevel) {
 		this.woi =          woi;
+		this.alpha_loss =   alpha_loss;
+		this.alpha_offset = alpha_offset;
 //		this.default_alpha = default_alpha;
 		int     min_scenes_uses = min_scenes;
 		int     min_scenes_used = min_scenes * 4;
@@ -114,16 +127,14 @@ public class VegetationLMA {
 		setupParametersVector (default_alpha); // areas where both terrain and vegetation are available		
 		setupDataSource(); // condensed
 		setupYVector();
-		setupWeights(scene_weights);
+		setupWeights( // after setupParametersIndices
+				scene_weights,
+				reg_weights);   // final double reg_weights);
 		last_jt = new double [parameters_vector.length][];
 		return weights.length;
 	}
 	
 	
-	
-	
-	
-	
 	private double [] getYminusFxWeighted(
 			final double []   fx,
 			final double []   rms_fp // null or [2]
@@ -154,6 +165,9 @@ public class VegetationLMA {
 //		System.out.println("ai.get()="+ai.get());
 		// important to set - after first cycle ai is left 16(number of threads) larger than number of cycles!
 		// It is so, because it first increments, then tests if (n < num_pairs)
+		
+		// Regularization part where y is assumed==0, only fX is used
+		
 		ai.set(y_vector.length); // not yet used
 		ati.set(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -312,22 +326,25 @@ public class VegetationLMA {
 		}
 		if (debug_level > 3) { // 0) {
 			double    delta =  this.delta; // 1E-3;
+			if (debug_level > 4) {
 				double delta_err= debugDerivs (
 						parameters_vector, // double [] vector,
 						delta,             // double    delta,
 						debug_level, //int       debugLevel,
 						false); // boolean   show_img)
 				System.out.println("\nMaximal error = "+delta_err);
+			}
 			int [] wh = new int[2];
 			double [] debug_image =  parametersImage (
 					parameters_vector, // double [] vector,
 					1,                // int       gap)
-					wh);
+					wh,
+					null); // String    title)
 			ShowDoubleFloatArrays.showArrays(
 					debug_image,
 					wh[0],
 					wh[1],
-					"parameters.vector");
+					"parameters.vector-"+debug_index);
 		}
 		
 		
@@ -402,13 +419,27 @@ public class VegetationLMA {
 			rslt[0] = true;
 			rslt[1] = rms[0] >=(this.last_rms[0] * (1.0 - rms_diff));
 			this.last_rms = rms.clone();
-
 			this.parameters_vector = new_vector.clone();
-			if (debug_level > 3) {
-				double [] debug_image =  parametersImage (
+			if (debug_level > 2) {
+				String title = (debug_level > 3)? "parameters.vector":null;
+				int [] wh = new int [2];
+				double [] debug_img =  parametersImage (
 						parameters_vector, // double [] vector,
 						1,                // int       gap)
-						null);
+						wh,
+						title);
+				if ((debug_image != null) && (debug_index < debug_image.length)) {
+					debug_image[debug_index++] = debug_img;
+				}
+				if (debug_level > 4) {
+					ShowDoubleFloatArrays.showArrays(
+							debug_image,
+							wh[0],
+							wh[1],
+							true,
+							"parameters_vector_live.tiff");
+					
+				}
 
 				// print vectors in some format
 				/*
@@ -514,7 +545,7 @@ public class VegetationLMA {
 				public void run() {
 					double [] vegetation = new double [4]; 
 					double [] alpha =      new double [4];
-					for (int n = ai.getAndIncrement(); n < weights.length; n = ai.getAndIncrement()) {
+					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]];
@@ -563,6 +594,35 @@ public class VegetationLMA {
 			};
 		}
 		ImageDtt.startAndJoin(threads);
+		// regularization weights and derivatives
+		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 + 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]		
+							}							
+						}
+					}
+				}
+			};
+		}
+		ImageDtt.startAndJoin(threads);
+		
 		return fX;		
 	}
 	
@@ -624,7 +684,8 @@ public class VegetationLMA {
 	public double [] parametersImage (
 			double [] vector,
 			int       gap,
-			int []    wh) {
+			int []    wh,
+			String    title) {
 		int [] indices = getParamDebugIndices (gap);
 		int woi_width3=woi.width*3 + 2 * gap;
 		double [] dbg_img = new double [indices.length];
@@ -633,18 +694,24 @@ public class VegetationLMA {
 				dbg_img[i] = Double.NaN;
 			} else {
 				dbg_img[i]= vector[indices[i]];
+				// scale alpha for the same image
+				int w = i % woi_width3;
+				if (w > 2 *( woi.width+ gap)) {
+					dbg_img[i] = debug_alpha_scales[0] + (debug_alpha_scales[1] - debug_alpha_scales[0]) * dbg_img[i];  
+				}
 			}
 		}
 		if (wh != null) {
 			wh[0] = woi_width3;
 			wh[1] = indices.length/woi_width3;
 		}
+		if (title != null) {
 			ShowDoubleFloatArrays.showArrays(
 					dbg_img,
 					woi_width3,
 					indices.length/woi_width3,
-				"parameters.vector");
-		
+					title);
+		}
 		return dbg_img;
 	}
 
@@ -683,9 +750,11 @@ public class VegetationLMA {
 	}
 	
 	
-	private void setupWeights(final double [] scene_weights) {
-		int extra_samples = 0; // in the future may be regularization
-		double extra_weights = 0.0;
+	private void setupWeights( // after setupParametersIndices
+			final double [] scene_weights,
+			final double reg_weights) {
+		int extra_samples = num_pars_vegetation_alpha; // in the future may be more regularization
+		double reg_sample_weight = reg_weights/num_pars_vegetation_alpha; // weight of each regularization sample
 		final double [] sw =  (scene_weights != null) ? scene_weights : new double [num_scenes];
 		if (scene_weights == null) {
 			Arrays.fill (sw, 1.0);
@@ -693,15 +762,17 @@ public class VegetationLMA {
 		weights = new double [data_source.length + extra_samples];
 		double s = 0;
 		for (int ny = 0; ny < data_source.length; ny++) {
-			s += sw[data_source[ny][0][0]];
+			int nscene =   data_source[ny][0][0];
+//			int indx =     data_source[ny][0][1];
+			s += sw[nscene];
 		}
 		final double s0 = s;
-		for (int i = 0; i < extra_samples; i++) {
-			s+=extra_weights;
-		}
+		s *= (1+ reg_weights);
 		final double k = 1.0/s;
 		weight_pure = s0/s;
-		
+//		for (int i = 0; i < extra_samples; i++) {
+//			s+=extra_weights;
+//		}
 		final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
 		final AtomicInteger ai = new AtomicInteger(0);
 		for (int ithread = 0; ithread < threads.length; ithread++) {
@@ -711,7 +782,7 @@ public class VegetationLMA {
 						if (nw  < data_source.length) {
 							weights[nw] = sw[data_source[nw][0][0]] * k;
 						} else {
-							weights[nw] = extra_weights * k;
+							weights[nw] = reg_sample_weight;
 						}
 					}
 				}
@@ -872,6 +943,7 @@ public class VegetationLMA {
 		Arrays.fill(par_index[TVAO_VEGETATION],       -1);
 		Arrays.fill(par_index[TVAO_VEGETATION_ALPHA], -1);
 		int par_num = 0;
+		ind_pars_terrain = 0;
 		for (int drow = 0; drow < woi.height; drow++) {
 			int row = woi.y + drow;
 			for (int dcol = 0; dcol < woi.width; dcol++) {
@@ -884,6 +956,7 @@ public class VegetationLMA {
 			}
 		}
 		num_pars_terrain =    par_num;
+		ind_pars_vegetation = par_num;
 		for (int drow = 0; drow < woi.height; drow++) {
 			int row = woi.y + drow;
 			for (int dcol = 0; dcol < woi.width; dcol++) {
@@ -895,7 +968,8 @@ public class VegetationLMA {
 				}
 			}
 		}
-		num_pars_vegetation = par_num - num_pars_terrain;
+		ind_pars_vegetation_alpha = par_num;
+		num_pars_vegetation = par_num - ind_pars_vegetation;
 		for (int drow = 0; drow < woi.height; drow++) {
 			int row = woi.y + drow;
 			for (int dcol = 0; dcol < woi.width; dcol++) {
@@ -907,11 +981,13 @@ public class VegetationLMA {
 				}
 			}
 		}
-		int num_pars_pixel = par_num;
+		num_pars_vegetation_alpha = par_num - ind_pars_vegetation_alpha; // normally should be == num_pars_vegetation
+//		int num_pars_pixel = par_num;
+		ind_pars_scenes =    par_num;
 		for (int i = 0; i < par_index[TVAO_SCENE_OFFSET].length; i++) {
 			par_index[TVAO_SCENE_OFFSET][i] = par_num++;
 		}
-		num_pars_scenes = par_num - num_pars_pixel;
+		num_pars_scenes = par_num - ind_pars_scenes;
 		parameters_vector = new double [par_num];
 	}
 	

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

@@ -698,11 +698,18 @@ public class VegetationModel {
 		int       min_scenes =  10;
 		double    default_alpha = 0.8;
 		double [] scene_weights = null;
+		double    reg_weights =   0.25;        // fraction of the total weight used for regularization
+		double    alpha_loss =    1.0;         // quadratic loss when alpha reaches -1.0 or 2.0 
+		double    alpha_offset =  0.0;         // if >0,  start losses above 0.0 and below 1.0;
+		
 		int       num_samples = vegetationLMA.prepareLMA(
 				woi50,         // final Rectangle   woi,
 				min_scenes,    // final int         min_scenes, // minimal number of scenes (inside woi) vegetation pixel must influence
 				default_alpha, // final double      default_alpha,
 				scene_weights, // final double []   scene_weights, // null or per-scene weights (higher for larger offsets)?
+				reg_weights,   // final double             reg_weights,        // fraction of the total weight used for regularization
+				alpha_loss,    // final double             alpha_loss,         // quadratic loss when alpha reaches -1.0 or 2.0 
+				alpha_offset,  // final double             alpha_offset,         // quadratic loss when alpha reaches -1.0 or 2.0 
 				debugLevel);   // final int         debugLevel);
 		
 		double  lambda =            0.1;
@@ -712,6 +719,9 @@ public class VegetationModel {
 		boolean last_run =       false;
 		double  rms_diff =          0.0001;
 		int     num_iter = 100;
+		vegetationLMA.debug_index = 0;
+		vegetationLMA.debug_image = new double [num_iter][];
+		
 		int lma_rslt= vegetationLMA.runLma( // <0 - failed, >=0 iteration number (1 - immediately)
 				lambda,           // double lambda,           // 0.1
 				lambda_scale_good,// double lambda_scale_good,// 0.5