Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java @ 3011:ab81ffd1343e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
FixA: Reactivated rewrite of the outlier checks. flys-artifacts/trunk@4576 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 04 Jun 2012 16:44:56 +0000
parents 05a3fe8800b3
children e341606faeb4
comparison
equal deleted inserted replaced
3010:05a3fe8800b3 3011:ab81ffd1343e
23 import de.intevation.flys.artifacts.model.FixingsOverviewFactory; 23 import de.intevation.flys.artifacts.model.FixingsOverviewFactory;
24 import de.intevation.flys.artifacts.model.Parameters; 24 import de.intevation.flys.artifacts.model.Parameters;
25 25
26 import de.intevation.flys.utils.DateAverager; 26 import de.intevation.flys.utils.DateAverager;
27 import de.intevation.flys.utils.DoubleUtil; 27 import de.intevation.flys.utils.DoubleUtil;
28 import de.intevation.flys.utils.Pair; 28 import de.intevation.flys.utils.EpsilonComparator;
29 29
30 import java.util.ArrayList; 30 import java.util.ArrayList;
31 import java.util.Date; 31 import java.util.Date;
32 import java.util.HashMap; 32 import java.util.HashMap;
33 import java.util.List; 33 import java.util.List;
34 import java.util.Map; 34 import java.util.Map;
35 35 import java.util.TreeMap;
36 import org.apache.commons.math.MathException;
37
38 import org.apache.commons.math.optimization.fitting.CurveFitter;
39
40 import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
41 36
42 import org.apache.log4j.Logger; 37 import org.apache.log4j.Logger;
43 38
44 public class FixCalculation 39 public class FixCalculation
45 extends Calculation 40 extends Calculation
46 { 41 {
47 private static Logger log = Logger.getLogger(FixCalculation.class); 42 private static Logger log = Logger.getLogger(FixCalculation.class);
43
44 public static final double EPSILON = 1e-4;
48 45
49 protected String river; 46 protected String river;
50 protected double from; 47 protected double from;
51 protected double to; 48 protected double to;
52 protected double step; 49 protected double step;
69 String function = access.getFunction(); 66 String function = access.getFunction();
70 int [] events = access.getEvents(); 67 int [] events = access.getEvents();
71 DateRange [] analysisPeriods = access.getAnalysisPeriods(); 68 DateRange [] analysisPeriods = access.getAnalysisPeriods();
72 Integer qSectorStart = access.getQSectorStart(); 69 Integer qSectorStart = access.getQSectorStart();
73 Integer qSectorEnd = access.getQSectorEnd(); 70 Integer qSectorEnd = access.getQSectorEnd();
71 Boolean preprocessing = access.getPreprocessing();
74 72
75 if (river == null) { 73 if (river == null) {
76 // TODO: i18n 74 // TODO: i18n
77 addProblem("fix.missing.river"); 75 addProblem("fix.missing.river");
78 } 76 }
113 } 111 }
114 112
115 if (qSectorEnd == null) { 113 if (qSectorEnd == null) {
116 // TODO: i18n 114 // TODO: i18n
117 addProblem("fix.missing.qend.sector"); 115 addProblem("fix.missing.qend.sector");
116 }
117
118 if (preprocessing == null) {
119 // TODO: i18n
120 addProblem("fix.missing.preprocessing");
118 } 121 }
119 122
120 if (!hasProblems()) { 123 if (!hasProblems()) {
121 this.river = river; 124 this.river = river;
122 this.from = from; 125 this.from = from;
125 this.function = function; 128 this.function = function;
126 this.events = events; 129 this.events = events;
127 this.analysisPeriods = analysisPeriods; 130 this.analysisPeriods = analysisPeriods;
128 this.qSectorStart = qSectorStart; 131 this.qSectorStart = qSectorStart;
129 this.qSectorEnd = qSectorEnd; 132 this.qSectorEnd = qSectorEnd;
133 this.preprocessing = preprocessing;
130 } 134 }
131 } 135 }
132 136
133 public CalculationResult calculate() { 137 public CalculationResult calculate() {
134 138
151 155
152 if (hasProblems()) { 156 if (hasProblems()) {
153 return new CalculationResult(this); 157 return new CalculationResult(this);
154 } 158 }
155 159
156 FixingsColumnFactory fcf = FixingsColumnFactory.getInstance(); 160 final List<Column> eventColumns = getEventColumns(overview);
157 161
158 List<FixingsColumn> dataColumns = 162 if (eventColumns.size() < 2) {
159 new ArrayList<FixingsColumn>(events.length);
160
161 for (int eventId: events) {
162 IdFilter idFilter = new IdFilter(eventId);
163
164 List<Fixing.Column> columns = overview.filter(null, idFilter);
165 if (columns.isEmpty()) {
166 // TODO: i18n
167 addProblem("fix.missing.column", eventId);
168 continue;
169 }
170 FixingsColumn dataColumn = fcf.getColumnData(columns.get(0));
171 if (dataColumn == null) {
172 // TODO: i18n
173 addProblem("fix.cannot.load.data", eventId);
174 continue;
175 }
176 dataColumns.add(dataColumn);
177 }
178
179 if (dataColumns.size() < 2) {
180 // TODO: i18n 163 // TODO: i18n
181 addProblem("fix.too.less.data.columns"); 164 addProblem("fix.too.less.data.columns");
182 return new CalculationResult(this); 165 return new CalculationResult(this);
183 } 166 }
184 167
185 double [] kms = DoubleUtil.explode(from, to, step / 1000.0); 168 double [] kms = DoubleUtil.explode(from, to, step / 1000.0);
186 169
187 double [] ws = new double[dataColumns.size()]; 170 final double [] qs = new double[eventColumns.size()];
188 double [] qs = new double[ws.length]; 171 final double [] ws = new double[qs.length];
172
173 // Depending on preprocessing we need to find the outliers.
174 Fitting.QWFactory qwFactory = !preprocessing
175 ? null // No outliers
176 : new Fitting.QWFactory() {
177 @Override
178 public QW create(double q, double w) {
179 // Check all the event columns for close match
180 // and take the description and the date from meta.
181 for (int i = 0; i < qs.length; ++i) {
182 if (Math.abs(qs[i]-q) < EPSILON
183 && Math.abs(ws[i]-w) < EPSILON) {
184 Column column = eventColumns.get(i);
185 return new QW(
186 q, w,
187 column.getDescription(),
188 column.getDate());
189 }
190 }
191 log.warn("cannot find column for (" + q + ", " + w + ")");
192 return new QW(q, w);
193 }
194 };
195
196 Fitting fitting = new Fitting(func, qwFactory);
189 197
190 String [] parameterNames = func.getParameterNames(); 198 String [] parameterNames = func.getParameterNames();
191 199
192 Parameters results = 200 Parameters results =
193 new Parameters(createColumnNames(parameterNames)); 201 new Parameters(createColumnNames(parameterNames));
196 204
197 if (debug) { 205 if (debug) {
198 log.debug("number of kms: " + kms.length); 206 log.debug("number of kms: " + kms.length);
199 } 207 }
200 208
201 int kmIndex = results.columnIndex("km"); 209 TreeMap<Double, QW []> outliers =
202 int chiSqrIndex = results.columnIndex("chi_sqr"); 210 new TreeMap<Double, QW []>(EpsilonComparator.INSTANCE);
211
212 int kmIndex = results.columnIndex("km");
213 int chiSqrIndex = results.columnIndex("chi_sqr");
214 int [] parameterIndices = results.columnIndices(parameterNames);
203 215
204 int numFailed = 0; 216 int numFailed = 0;
205 217
206 for (int i = 0; i < kms.length; ++i) { 218 for (int i = 0; i < kms.length; ++i) {
207 double km = kms[i]; 219 double km = kms[i];
208 220
221 // Fill Qs and Ws from event columns.
209 for (int j = 0; j < ws.length; ++j) { 222 for (int j = 0; j < ws.length; ++j) {
210 FixingsColumn column = dataColumns.get(j); 223 boolean interpolated =
211 qs[j] = column.getQ(km); 224 eventColumns.get(j).getQW(km, qs, ws, j);
212 boolean interpolated = column.getW(km, ws, j);
213 // TODO: mark as interpolated. 225 // TODO: mark as interpolated.
214 } 226 }
215 227
216 Pair<double [], Double> fitResult = fit(func, km, ws, qs); 228 fitting.reset();
217 229
218 // TODO: Do preprocessing here! 230 if (!fitting.fit(qs, ws)) {
219 if (fitResult == null) { // Problems are reported already.
220 ++numFailed; 231 ++numFailed;
232 // TODO: i18n
233 addProblem(km, "fix.fitting.failed");
221 continue; 234 continue;
222 } 235 }
223 236
237 if (fitting.hasOutliers()) {
238 outliers.put(Double.valueOf(km), fitting.outliersToArray());
239 }
240
224 int row = results.newRow(); 241 int row = results.newRow();
225 242
226 results.set(row, kmIndex, km); 243 results.set(row, kmIndex, km);
227 results.set(row, chiSqrIndex, fitResult.getB()); 244 results.set(row, chiSqrIndex, fitting.getChiSquare());
228 245 invalid |= results.set(
229 double [] parameters = fitResult.getA(); 246 row, parameterIndices, fitting.getParameters());
230 for (int j = 0; j < parameters.length; ++j) {
231 if (Double.isNaN(parameters[j])) {
232 invalid = true;
233 }
234 else {
235 results.set(row, parameterNames[j], parameters[j]);
236 }
237 }
238 } 247 }
239 248
240 if (debug) { 249 if (debug) {
241 log.debug("success: " + (kms.length - numFailed)); 250 log.debug("success: " + (kms.length - numFailed));
242 log.debug("failed: " + numFailed); 251 log.debug("failed: " + numFailed);
252 DeltaWTsKM deltaWTsKM = calculateDeltaWTs( 261 DeltaWTsKM deltaWTsKM = calculateDeltaWTs(
253 func, 262 func,
254 overview, 263 overview,
255 results); 264 results);
256 265
257 FixResult fr = new FixResult(results, deltaWTsKM); 266 FixResult fr = new FixResult(results, deltaWTsKM, outliers);
258 267
259 return new CalculationResult(fr, this); 268 return new CalculationResult(fr, this);
269 }
270
271
272 protected List<Column> getEventColumns(FixingsOverview overview) {
273
274 FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();
275
276 List<Column> columns = new ArrayList<Column>(events.length);
277
278 for (int eventId: events) {
279 IdFilter idFilter = new IdFilter(eventId);
280
281 List<Fixing.Column> metas = overview.filter(null, idFilter);
282
283 if (metas.isEmpty()) {
284 // TODO: i18n
285 addProblem("fix.missing.column", eventId);
286 continue;
287 }
288
289 FixingsColumn data = fcf.getColumnData(metas.get(0));
290 if (data == null) {
291 // TODO: i18n
292 addProblem("fix.cannot.load.data", eventId);
293 continue;
294 }
295
296 columns.add(new Column(metas.get(0), data));
297 }
298
299 return columns;
260 } 300 }
261 301
262 public DeltaWTsKM calculateDeltaWTs( 302 public DeltaWTsKM calculateDeltaWTs(
263 Function function, 303 Function function,
264 FixingsOverview overview, 304 FixingsOverview overview,
411 // Ignore not matching sectors. 451 // Ignore not matching sectors.
412 continue; 452 continue;
413 } 453 }
414 454
415 Column column = cc.getColumn(meta); 455 Column column = cc.getColumn(meta);
416 if (column == null || !column.getWQ(km, wq)) { 456 if (column == null || !column.getQW(km, wq)) {
417 continue; 457 continue;
418 } 458 }
419 459
420 double fw = instance.value(wq[1]); 460 double fw = instance.value(wq[1]);
421 if (Double.isNaN(fw)) { 461 if (Double.isNaN(fw)) {
425 double dw = (wq[0] - fw)*100.0; 465 double dw = (wq[0] - fw)*100.0;
426 466
427 Date date = column.getDate(); 467 Date date = column.getDate();
428 String description = column.getDescription(); 468 String description = column.getDescription();
429 469
430 QWD qwd = new QWD(wq[1], wq[0], dw, date, description); 470 QWD qwd = new QWD(wq[1], wq[0], description, date, dw);
431 471
432 qwds.add(qwd); 472 qwds.add(qwd);
433 473
434 sumW += wq[0]; 474 sumW += wq[0];
435 sumQ += wq[1]; 475 sumQ += wq[1];
449 Date avgDate = dateAverager.getAverage(); 489 Date avgDate = dateAverager.getAverage();
450 490
451 String avgDescription = "avg.deltawt." + qSector; 491 String avgDescription = "avg.deltawt." + qSector;
452 492
453 QWD avgQWD = new QWD( 493 QWD avgQWD = new QWD(
454 avgQ, avgW, avgDw, avgDate, avgDescription); 494 avgQ, avgW, avgDescription, avgDate, avgDw);
455 495
456 // TODO: Store average value. 496 // TODO: Store average value.
457 } 497 }
458 else { 498 else {
459 // TODO: Store 499 // TODO: Store
488 528
489 public String getDescription() { 529 public String getDescription() {
490 return meta.getDescription(); 530 return meta.getDescription();
491 } 531 }
492 532
493 public boolean getWQ(double km, double [] wq) { 533 public boolean getQW(
534 double km,
535 double [] qs,
536 double [] ws,
537 int index
538 ) {
539 qs[index] = data.getQ(km);
540 return data.getW(km, ws, index);
541 }
542
543 public boolean getQW(double km, double [] wq) {
494 data.getW(km, wq, 0); 544 data.getW(km, wq, 0);
495 if (Double.isNaN(wq[0])) return false; 545 if (Double.isNaN(wq[0])) return false;
496 wq[1] = data.getQ(km); 546 wq[1] = data.getQ(km);
497 return !Double.isNaN(wq[1]); 547 return !Double.isNaN(wq[1]);
498 } 548 }
585 result[0] = "km"; 635 result[0] = "km";
586 result[0] = "chi_sqr"; 636 result[0] = "chi_sqr";
587 System.arraycopy(parameters, 0, result, 2, parameters.length); 637 System.arraycopy(parameters, 0, result, 2, parameters.length);
588 return result; 638 return result;
589 } 639 }
590
591 protected Pair<double [], Double> fit(
592 Function function,
593 double km,
594 double [] ws,
595 double [] qs
596 ) {
597 LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();
598 CurveFitter cf = new CurveFitter(lmo);
599
600 for (int i = 0; i < ws.length; ++i) {
601 if (!Double.isNaN(ws[i]) && !Double.isNaN(qs[i])) {
602 cf.addObservedPoint(qs[i], ws[i]);
603 }
604 }
605
606 try {
607 double [] parameters =
608 cf.fit(function, function.getInitialGuess());
609 double chiSqr = lmo.getChiSquare();
610
611 return new Pair<double [], Double>(parameters, chiSqr);
612 }
613 catch (MathException me) {
614 log.warn(me, me);
615 addProblem(km, "fix.fitting.failed");
616 }
617
618 return null;
619 }
620 } 640 }
621 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : 641 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org