Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java @ 9646:0380717105ba

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implemented alternative fitting strategy for Log-Linear function.
author Gernot Belger <g.belger@bjoernsen.de>
date Mon, 02 Dec 2019 17:56:15 +0100
parents bc50ecfc58c5
children
comparison
equal deleted inserted replaced
9645:eb1a29fe823f 9646:0380717105ba
10 10
11 import java.util.ArrayList; 11 import java.util.ArrayList;
12 import java.util.Date; 12 import java.util.Date;
13 import java.util.List; 13 import java.util.List;
14 14
15 import org.apache.commons.math.MathException;
16 import org.apache.commons.math.optimization.fitting.CurveFitter;
17 import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
18 import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
19 import org.apache.log4j.Logger; 15 import org.apache.log4j.Logger;
20 import org.dive4elements.river.artifacts.math.GrubbsOutlier; 16 import org.dive4elements.river.artifacts.math.GrubbsOutlier;
21 import org.dive4elements.river.artifacts.math.fitting.Function; 17 import org.dive4elements.river.artifacts.math.fitting.Function;
18 import org.dive4elements.river.artifacts.math.fitting.LogLinearAlternative;
22 19
23 public class Fitting { 20 public class Fitting {
24 private static Logger log = Logger.getLogger(Fitting.class); 21 private static Logger log = Logger.getLogger(Fitting.class);
25 22
26 /** Use instance of this factory to find meta infos for outliers. */ 23 /** Use instance of this factory to find meta infos for outliers. */
27 public interface QWDFactory { 24 public interface QWDFactory {
28 QWD create(double q, double w, double deltaW, boolean isOutlier); 25 QWD create(double q, double w, double deltaW, boolean isOutlier);
29 }
30
31 private static final class FittingData {
32 public final FixingColumnWithData event;
33 public final double w;
34 public final double q;
35 public final boolean isInterpolated;
36
37 public FittingData(final FixingColumnWithData event, final double w, final double q, final boolean isInterpolated) {
38 this.event = event;
39 this.w = w;
40 this.q = q;
41 this.isInterpolated = isInterpolated;
42 }
43 } 26 }
44 27
45 private final double chiSqr; 28 private final double chiSqr;
46 29
47 private final double[] parameters; 30 private final double[] parameters;
99 } 82 }
100 83
101 final List<FittingData> outliers = new ArrayList<>(data.size()); 84 final List<FittingData> outliers = new ArrayList<>(data.size());
102 85
103 org.dive4elements.river.artifacts.math.Function instance = null; 86 org.dive4elements.river.artifacts.math.Function instance = null;
104 LevenbergMarquardtOptimizer lmo = null; 87 Fitting fitting = null;
105 double[] parameters = null;
106 88
89 final IFittingOperation operation = createOperation(function);
107 for (;;) { 90 for (;;) {
108 91
109 parameters = null; 92 fitting = operation.execute(data);
110 93
111 for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) { 94 if (fitting == null) {
112
113 lmo = new LevenbergMarquardtOptimizer();
114 lmo.setCostRelativeTolerance(tolerance);
115 lmo.setOrthoTolerance(tolerance);
116 lmo.setParRelativeTolerance(tolerance);
117
118 try {
119 final CurveFitter cf = new CurveFitter(lmo);
120
121 for (final FittingData fittingData : data)
122 cf.addObservedPoint(fittingData.q, fittingData.w);
123
124 parameters = cf.fit(function, function.getInitialGuess());
125 break;
126 }
127 catch (final MathException me) {
128 if (log.isDebugEnabled()) {
129 log.debug("tolerance " + tolerance + " + failed.", me);
130 }
131 }
132 }
133
134 if (parameters == null) {
135 /* 95 /*
136 * log.debug("Parameters is null"); 96 * log.debug("Parameters is null");
137 * for (int i = 0, N = xs.size(); i < N; ++i) { 97 * for (int i = 0, N = xs.size(); i < N; ++i) {
138 * log.debug("DATA: " + xs.getQuick(i) + " " + ys.getQuick(i)); 98 * log.debug("DATA: " + xs.getQuick(i) + " " + ys.getQuick(i));
139 * } 99 * }
140 */ 100 */
141 return null; 101 return null;
142 } 102 }
143 103
144 // This is the parameterized function for a given km. 104 // This is the parameterized function for a given km.
105 final double[] parameters = fitting.getParameters();
145 instance = function.instantiate(parameters); 106 instance = function.instantiate(parameters);
146 107
147 if (!checkOutliers) 108 if (!checkOutliers)
148 break; 109 break;
149 110
178 qwds.add(qwd); 139 qwds.add(qwd);
179 resultColumns.addQWD(outlier.event, km, qwd); 140 resultColumns.addQWD(outlier.event, km, qwd);
180 } 141 }
181 142
182 /* 143 /*
183 * calculate dW of used values against the resulting function and add them to results , also calculate standard * 144 * calculate dW of used values against the resulting function and add them to results
184 * deviation
185 */ 145 */
186 final StandardDeviation stdDev = new StandardDeviation();
187 double maxQ = -Double.MAX_VALUE;
188
189 for (final FittingData fittingData : data) { 146 for (final FittingData fittingData : data) {
190 147
191 final QWD qwd = createQWD(fittingData, instance, false); 148 final QWD qwd = createQWD(fittingData, instance, false);
192 qwds.add(qwd); 149 qwds.add(qwd);
193 resultColumns.addQWD(fittingData.event, km, qwd); 150 resultColumns.addQWD(fittingData.event, km, qwd);
194
195 stdDev.increment(qwd.getDeltaW());
196
197 final double q = qwd.getQ();
198 if (q > maxQ)
199 maxQ = q;
200 } 151 }
201 152
202 final double standardDeviation = stdDev.getResult(); 153 return fitting;
154 }
203 155
204 final double chiSqr = lmo.getChiSquare(); 156 private static IFittingOperation createOperation(final Function function) {
205 157
206 return new Fitting(parameters, standardDeviation, chiSqr, maxQ); 158 if (function instanceof LogLinearAlternative)
159 return new LogLinearFittingOperation((LogLinearAlternative) function);
160
161 return new LevenbergMarquardtFittingOperation(function);
207 } 162 }
208 163
209 private static QWD createQWD(final FittingData data, final org.dive4elements.river.artifacts.math.Function function, final boolean isOutlier) { 164 private static QWD createQWD(final FittingData data, final org.dive4elements.river.artifacts.math.Function function, final boolean isOutlier) {
210 165
211 final FixingColumnWithData event = data.event; 166 final FixingColumnWithData event = data.event;
http://dive4elements.wald.intevation.org