[PATCH] S/Q relation: More Excel compat

[Wald Commits]

# HG changeset patch
# User Sascha L. Teichmann <teichmann at intevation.de>
# Date 1375896965 -7200
# Node ID 48f6780c372d4e73a56b3c2737b3a496a29f5a7f
# Parent  c146fd412f9d379fcc42bde727891e4581d3d048
S/Q relation: More Excel compat.

diff -r c146fd412f9d -r 48f6780c372d artifacts/src/main/java/org/dive4elements/river/artifacts/model/sq/Fitting.java
--- a/artifacts/src/main/java/org/dive4elements/river/artifacts/model/sq/Fitting.java	Wed Aug 07 18:24:07 2013 +0200
+++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/model/sq/Fitting.java	Wed Aug 07 19:36:05 2013 +0200
@@ -9,7 +9,6 @@
 package org.dive4elements.river.artifacts.model.sq;
 
 import org.dive4elements.river.artifacts.math.fitting.Function;
-import org.dive4elements.river.artifacts.math.fitting.Linear;
 
 import java.util.ArrayList;
 import java.util.List;
@@ -18,9 +17,8 @@
 
 import org.apache.commons.math.optimization.fitting.CurveFitter;
 
-import org.apache.commons.math.optimization.general.AbstractLeastSquaresOptimizer;
-import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
 import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
+import org.apache.commons.math.stat.regression.SimpleRegression;
 
 import org.apache.log4j.Logger;
 
@@ -82,53 +80,75 @@
     @Override
     public void initialize(List<SQ> sqs) throws MathException {
 
-        AbstractLeastSquaresOptimizer optimizer = getOptimizer();
+        if (USE_NON_LINEAR_FITTING
+        || function.getInitialGuess().length != 2) {
+            nonLinearFitting(sqs);
+        }
+        else {
+            linearFitting(sqs);
+        }
+    }
+
+    protected void linearFitting(List<SQ> sqs) {
+
+        coeffs = linearRegression(sqs);
+
+        instance = function.instantiate(coeffs);
+    }
+
+    protected double [] linearRegression(List<SQ> sqs) {
+
+        SimpleRegression reg = new SimpleRegression();
+
+        int invalidPoints = 0;
+        for (SQ sq: sqs) {
+            double s = sq.getS();
+            double q = sq.getQ();
+            if (s <= 0d || q <= 0d) {
+                ++invalidPoints;
+                continue;
+            }
+            reg.addData(Math.log(q), Math.log(s));
+        }
+
+        if (sqs.size() - invalidPoints < 2) {
+            log.debug("not enough points");
+            return new double [] { 0, 0 };
+        }
+
+        double a = Math.exp(reg.getIntercept());
+        double b = reg.getSlope();
+
+        if (log.isDebugEnabled()) {
+            log.debug("invalid points: " +
+                invalidPoints + " (" + sqs.size() + ")");
+            log.debug("a: " + a + " (" + Math.log(a) + ")");
+            log.debug("b: " + b);
+        }
+
+        return new double [] { a, b };
+    }
+
+
+    protected void nonLinearFitting(List<SQ> sqs) throws MathException {
+
+        LevenbergMarquardtOptimizer optimizer =
+            new LevenbergMarquardtOptimizer();
 
         CurveFitter cf = new CurveFitter(optimizer);
-        double [] values = new double[2];
+
         for (SQ sq: sqs) {
-            values[0] = sq.getQ();
-            values[1] = sq.getS();
-            transformInputValues(values);
-            cf.addObservedPoint(values[0], values[1]);
+            cf.addObservedPoint(sq.getS(), sq.getQ());
         }
 
         coeffs = cf.fit(
             function, function.getInitialGuess());
 
-        transformCoeffsBack(coeffs);
-
         instance = function.instantiate(coeffs);
 
         chiSqr = optimizer.getChiSquare();
     }
 
-    protected Function getFunction(Function function) {
-        return USE_NON_LINEAR_FITTING
-            ? function
-            : Linear.INSTANCE;
-    }
-
-    protected void transformInputValues(double [] values) {
-        if (!USE_NON_LINEAR_FITTING) {
-            for (int i = 0; i < values.length; ++i) {
-                values[i] = Math.log(values[i]);
-            }
-        }
-    }
-
-    protected AbstractLeastSquaresOptimizer getOptimizer() {
-        return USE_NON_LINEAR_FITTING
-            ? new LevenbergMarquardtOptimizer()
-            : new GaussNewtonOptimizer(false);
-    }
-
-    protected void transformCoeffsBack(double [] coeffs) {
-        if (!USE_NON_LINEAR_FITTING && coeffs.length > 0) {
-            coeffs[0] = Math.exp(coeffs[0]);
-        }
-    }
-
     @Override
     public double eval(SQ sq) {
         double s = instance.value(sq.q);