Improved wood yield function fitting.

src/ac/ed/lurg/forestry/WoodYieldResponse.java

@@ -50,7 +50,7 @@ public class WoodYieldResponse {
 
 		// Use Gauss-Newton algorithm to improve fit
 		Double[] initialGuess = {maxYield, slope, 1.0};
-		Double[] fittedParams = gaussNewton(scaledYield, initialGuess, 1000, 1e-6);
+		Double[] fittedParams = gaussNewton(scaledYield, initialGuess, 1000, 0.001, 0.1);
 
 		// If fit improved update parameters. Sometimes the GN algorithm fails if fit is very poor so use initial fit instead.
 		double resFitted = calcResidualSqr(scaledYield, fittedParams[0], fittedParams[1], fittedParams[2]);
@@ -120,7 +120,7 @@ public class WoodYieldResponse {
 	}
 
 	// Gauss-Newton optimization algorithm
-	private Double[] gaussNewton(Map<Integer, Double> yields, Double[] initialGuess, int maxIterations, double tolerance) {
+	private Double[] gaussNewton(Map<Integer, Double> yields, Double[] initialGuess, int maxIterations, double tolerance, double stepSize) {
 		double a = initialGuess[0];
 		double b = initialGuess[1];
 		double c = initialGuess[2];
@@ -131,7 +131,10 @@ public class WoodYieldResponse {
 			double sumBB = 0, sumBC = 0, sumCC = 0;
 
 			for (Map.Entry<Integer, Double> entry : yields.entrySet()) {
-				double t = Math.max(entry.getKey(), 0.1);
+				double t = entry.getKey();
+				if (t == 0) {
+					continue;
+				}
 				double y = entry.getValue();
 				double f = yieldFunction(a, b, c, t);
 				double error = y - f;
@@ -163,16 +166,16 @@ public class WoodYieldResponse {
 				break;
 			}
 
-			a += delta[0];
-			b += delta[1];
-			c += delta[2];
+			a += delta[0] * stepSize;
+			b += delta[1] * stepSize;
+			c += delta[2] * stepSize;
 
-			a = Math.max(a, 1e-6);
+			a = Math.max(a, 0.01);
 			a = Math.min(a, 50);
-			b = Math.max(-1, b);
-			b = Math.min(-1e-6, b);
-			c = Math.max(1e-6, c);
-			c = Math.min(10, c);
+			b = Math.max(-0.5, b);
+			b = Math.min(-0.001, b);
+			c = Math.max(0.5, c);
+			c = Math.min(3, c);
 
 			if (Math.abs(delta[0]) < tolerance && Math.abs(delta[1]) < tolerance && Math.abs(delta[2]) < tolerance) {
 				break;