Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/Fitting.java @ 9360:ddcd52d239cd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Outliers in fixation calculation are now shown within the other 'B' event themes and get a separate symbol (triangle). Removed old outliers theme. Also consider showpoints property. Also consider pointsize property.
author gernotbelger
date Wed, 01 Aug 2018 17:13:52 +0200
parents a3f318347707
children 9744ce3c3853
comparison
equal deleted inserted replaced
9359:df6f1b5806f6 9360:ddcd52d239cd
25 public class Fitting { 25 public class Fitting {
26 private static Logger log = Logger.getLogger(Fitting.class); 26 private static Logger log = Logger.getLogger(Fitting.class);
27 27
28 /** Use instance of this factory to find meta infos for outliers. */ 28 /** Use instance of this factory to find meta infos for outliers. */
29 public interface QWDFactory { 29 public interface QWDFactory {
30 QWD create(double q, double w, boolean isOutlier); 30 QWD create(double q, double w, double deltaW, boolean isOutlier);
31 } 31 }
32 32
33 private final boolean checkOutliers; 33 private final double chiSqr;
34
35 private final Function function;
36
37 private final QWDFactory qwdFactory;
38
39 private double chiSqr;
40
41 private double[] parameters;
42
43 private final List<QWD> removed = new ArrayList<>();
44
45 private QWD[] referenced;
46
47 private double standardDeviation;
48 34
49 public Fitting(Function function, QWDFactory qwdFactory, boolean checkOutliers) { 35 private final double[] parameters;
50 this.function = function;
51 this.qwdFactory = qwdFactory;
52 this.checkOutliers = checkOutliers;
53 }
54 36
55 public Function getFunction() { 37 private final double standardDeviation;
56 return function;
57 }
58 38
59 public boolean getCheckOutliers() { 39 private final List<QWD> qwds;
60 return checkOutliers; 40
41 public Fitting(final double[] parameters, final double standardDeviation, final double chiSqr, final List<QWD> qwds) {
42 this.parameters = parameters;
43 this.standardDeviation = standardDeviation;
44 this.chiSqr = chiSqr;
45 this.qwds = qwds;
61 } 46 }
62 47
63 public double getChiSquare() { 48 public double getChiSquare() {
64 return chiSqr; 49 return chiSqr;
65 } 50 }
66 51
67 public void reset() { 52 /**
68 chiSqr = 0.0; 53 * Returns all referenced and outliers as one array.
69 parameters = null; 54 */
70 removed.clear(); 55 public QWD[] getFixingsArray() {
71 referenced = null; 56 return qwds.toArray(new QWD[qwds.size()]);
72 standardDeviation = 0.0;
73 }
74
75 public boolean hasOutliers() {
76 return !removed.isEmpty();
77 }
78
79 public QWD[] outliersToArray() {
80 return removed.toArray(new QWD[removed.size()]);
81 }
82
83 public QWD[] referencedToArray() {
84 return referenced != null ? (QWD[]) referenced.clone() : null;
85 } 57 }
86 58
87 public double getMaxQ() { 59 public double getMaxQ() {
88 double maxQ = -Double.MAX_VALUE; 60 double maxQ = -Double.MAX_VALUE;
89 if (referenced != null) { 61
90 for (QWI qw : referenced) { 62 for (QWD qw : qwds) {
91 double q = qw.getQ(); 63 final double q = qw.getQ();
92 if (q > maxQ) { 64 if (!qw.isOutlier() && q > maxQ)
93 maxQ = q; 65 maxQ = q;
94 }
95 }
96 } 66 }
67
97 return maxQ; 68 return maxQ;
98 } 69 }
99 70
100 public double[] getParameters() { 71 public double[] getParameters() {
101 return parameters; 72 return parameters;
103 74
104 public double getStandardDeviation() { 75 public double getStandardDeviation() {
105 return standardDeviation; 76 return standardDeviation;
106 } 77 }
107 78
108 public boolean fit(double[] qs, double[] ws) { 79 public static Fitting fit(final Function function, final QWDFactory qwdFactory, final boolean checkOutliers, final double[] qs, final double[] ws) {
109 80
110 TDoubleArrayList xs = new TDoubleArrayList(qs.length); 81 final TDoubleArrayList xs = new TDoubleArrayList(qs.length);
111 TDoubleArrayList ys = new TDoubleArrayList(ws.length); 82 final TDoubleArrayList ys = new TDoubleArrayList(ws.length);
112 83
113 for (int i = 0; i < qs.length; ++i) { 84 for (int i = 0; i < qs.length; ++i) {
114 if (!Double.isNaN(qs[i]) && !Double.isNaN(ws[i])) { 85 if (!Double.isNaN(qs[i]) && !Double.isNaN(ws[i])) {
115 xs.add(qs[i]); 86 xs.add(qs[i]);
116 ys.add(ws[i]); 87 ys.add(ws[i]);
117 } 88 }
118 } 89 }
119 90
120 if (xs.size() < 2) { 91 if (xs.size() < 2) {
121 log.warn("Too less points."); 92 log.warn("Too less points.");
122 return false; 93 return null;
123 } 94 }
124 95
125 List<Double> inputs = new ArrayList<Double>(xs.size()); 96 final List<Double> inputs = new ArrayList<>(xs.size());
97 final List<QWD> qwds = new ArrayList<>(xs.size());
98 final List<QWD> outliers = new ArrayList<>(xs.size());
126 99
127 org.dive4elements.river.artifacts.math.Function instance = null; 100 org.dive4elements.river.artifacts.math.Function instance = null;
128
129 LevenbergMarquardtOptimizer lmo = null; 101 LevenbergMarquardtOptimizer lmo = null;
102 double[] parameters = null;
130 103
131 for (;;) { 104 for (;;) {
132 parameters = null; 105 parameters = null;
133 for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) { 106 for (double tolerance = 1e-10; tolerance < 1e-1; tolerance *= 10d) {
134 107
147 parameters = cf.fit(function, function.getInitialGuess()); 120 parameters = cf.fit(function, function.getInitialGuess());
148 break; 121 break;
149 } 122 }
150 catch (MathException me) { 123 catch (MathException me) {
151 if (log.isDebugEnabled()) { 124 if (log.isDebugEnabled()) {
152 log.debug("tolerance " + tolerance + " + failed."); 125 log.debug("tolerance " + tolerance + " + failed.", me);
153 } 126 }
154 } 127 }
155 } 128 }
129
156 if (parameters == null) { 130 if (parameters == null) {
157 /* 131 /*
158 * log.debug("Parameters is null"); 132 * log.debug("Parameters is null");
159 * for (int i = 0, N = xs.size(); i < N; ++i) { 133 * for (int i = 0, N = xs.size(); i < N; ++i) {
160 * log.debug("DATA: " + xs.getQuick(i) + " " + ys.getQuick(i)); 134 * log.debug("DATA: " + xs.getQuick(i) + " " + ys.getQuick(i));
161 * } 135 * }
162 */ 136 */
163 return false; 137 return null;
164 } 138 }
165 139
166 // This is the paraterized function for a given km. 140 // This is the paraterized function for a given km.
167 instance = function.instantiate(parameters); 141 instance = function.instantiate(parameters);
168 142
169 if (!checkOutliers) { 143 if (!checkOutliers)
170 break; 144 break;
171 }
172 145
173 inputs.clear(); 146 inputs.clear();
174 147
175 for (int i = 0, N = xs.size(); i < N; ++i) { 148 for (int i = 0, N = xs.size(); i < N; ++i) {
176 double y = instance.value(xs.getQuick(i)); 149 double y = instance.value(xs.getQuick(i));
178 y = Double.MAX_VALUE; 151 y = Double.MAX_VALUE;
179 } 152 }
180 inputs.add(Double.valueOf(ys.getQuick(i) - y)); 153 inputs.add(Double.valueOf(ys.getQuick(i) - y));
181 } 154 }
182 155
183 Integer outlier = GrubbsOutlier.findOutlier(inputs); 156 final Integer outlier = GrubbsOutlier.findOutlier(inputs);
157 if (outlier == null)
158 break;
184 159
185 if (outlier == null) { 160 final int idx = outlier.intValue();
186 break; 161 outliers.add(qwdFactory.create(xs.getQuick(idx), ys.getQuick(idx), Double.NaN, true));
187 }
188
189 int idx = outlier.intValue();
190 removed.add(qwdFactory.create(xs.getQuick(idx), ys.getQuick(idx), true));
191 xs.remove(idx); 162 xs.remove(idx);
192 ys.remove(idx); 163 ys.remove(idx);
193 } 164 }
165
166 for (QWD outlier : outliers) {
167
168 final double w = outlier.getW();
169 final double q = outlier.getQ();
170
171 final double dw = (w - instance.value(q)) * 100.0;
194 172
195 StandardDeviation stdDev = new StandardDeviation(); 173 outlier.setDeltaW(dw);
196 174
197 referenced = new QWD[xs.size()]; 175 qwds.add(outlier);
198 for (int i = 0; i < referenced.length; ++i) {
199 QWD qwd = qwdFactory.create(xs.getQuick(i), ys.getQuick(i), false);
200
201 if (qwd == null) {
202 log.warn("QW creation failed!");
203 } else {
204 referenced[i] = qwd;
205 double dw = (qwd.getW() - instance.value(qwd.getQ())) * 100.0;
206 qwd.setDeltaW(dw);
207 stdDev.increment(dw);
208 }
209 } 176 }
210 177
211 standardDeviation = stdDev.getResult(); 178 final StandardDeviation stdDev = new StandardDeviation();
212 179
213 chiSqr = lmo.getChiSquare(); 180 for (int i = 0; i < xs.size(); ++i) {
214 181
215 return true; 182 final QWD qwd = qwdFactory.create(xs.getQuick(i), ys.getQuick(i), Double.NaN, false);
183
184 double dw = (qwd.getW() - instance.value(qwd.getQ())) * 100.0;
185 qwd.setDeltaW(dw);
186
187 qwds.add(qwd);
188
189 stdDev.increment(dw);
190 }
191
192 final double standardDeviation = stdDev.getResult();
193
194 final double chiSqr = lmo.getChiSquare();
195
196 return new Fitting(parameters, standardDeviation, chiSqr, qwds);
216 } 197 }
217 } 198 }
http://dive4elements.wald.intevation.org