Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixAnalysisCalculation.java @ 3411:0ef83077c93f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
FixA: Renamed FixCalculation to FixAnalysisCalculation. flys-artifacts/trunk@5064 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 19 Jul 2012 14:11:11 +0000
parents
children e3c7a3228bc2
comparison
equal deleted inserted replaced
3410:f382127df48e 3411:0ef83077c93f
1 package de.intevation.flys.artifacts.model.fixings;
2
3 import de.intevation.flys.artifacts.access.FixAnalysisAccess;
4
5 import de.intevation.flys.artifacts.math.fitting.Function;
6 import de.intevation.flys.artifacts.math.fitting.FunctionFactory;
7
8 import de.intevation.flys.artifacts.model.Calculation;
9 import de.intevation.flys.artifacts.model.CalculationResult;
10 import de.intevation.flys.artifacts.model.DateRange;
11 import de.intevation.flys.artifacts.model.FixingsColumn;
12 import de.intevation.flys.artifacts.model.FixingsColumnFactory;
13
14 import de.intevation.flys.artifacts.model.FixingsOverview.AndFilter;
15 import de.intevation.flys.artifacts.model.FixingsOverview.DateRangeFilter;
16
17 import de.intevation.flys.artifacts.model.FixingsOverview.Fixing.Filter;
18
19 import de.intevation.flys.artifacts.model.FixingsOverview.Fixing;
20 import de.intevation.flys.artifacts.model.FixingsOverview.IdsFilter;
21 import de.intevation.flys.artifacts.model.FixingsOverview.KmFilter;
22 import de.intevation.flys.artifacts.model.FixingsOverview.SectorFilter;
23 import de.intevation.flys.artifacts.model.FixingsOverview.SectorRangeFilter;
24
25 import de.intevation.flys.artifacts.model.FixingsOverview;
26 import de.intevation.flys.artifacts.model.FixingsOverviewFactory;
27 import de.intevation.flys.artifacts.model.Parameters;
28 import de.intevation.flys.artifacts.model.Range;
29
30 import de.intevation.flys.utils.DateAverager;
31 import de.intevation.flys.utils.DoubleUtil;
32 import de.intevation.flys.utils.KMIndex;
33
34 import java.util.ArrayList;
35 import java.util.Date;
36 import java.util.HashMap;
37 import java.util.List;
38 import java.util.Map;
39
40 import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
41
42 import org.apache.log4j.Logger;
43
44 public class FixAnalysisCalculation
45 extends Calculation
46 {
47 private static Logger log = Logger.getLogger(FixAnalysisCalculation.class);
48
49 public static final double EPSILON = 1e-4;
50
51 protected String river;
52 protected double from;
53 protected double to;
54 protected double step;
55 protected boolean preprocessing;
56 protected String function;
57 protected int [] events;
58 protected DateRange referencePeriod;
59 protected DateRange [] analysisPeriods;
60 protected int qSectorStart;
61 protected int qSectorEnd;
62
63 public FixAnalysisCalculation() {
64 }
65
66 public FixAnalysisCalculation(FixAnalysisAccess access) {
67
68 String river = access.getRiver();
69 Double from = access.getFrom();
70 Double to = access.getTo();
71 Double step = access.getStep();
72 String function = access.getFunction();
73 int [] events = access.getEvents();
74 DateRange referencePeriod = access.getReferencePeriod();
75 DateRange [] analysisPeriods = access.getAnalysisPeriods();
76 Integer qSectorStart = access.getQSectorStart();
77 Integer qSectorEnd = access.getQSectorEnd();
78 Boolean preprocessing = access.getPreprocessing();
79
80 if (river == null) {
81 addProblem("fix.missing.river");
82 }
83
84 if (from == null) {
85 addProblem("fix.missing.from");
86 }
87
88 if (to == null) {
89 addProblem("fix.missing.to");
90 }
91
92 if (step == null) {
93 addProblem("fix.missing.step");
94 }
95
96 if (function == null) {
97 addProblem("fix.missing.function");
98 }
99
100 if (events == null || events.length < 1) {
101 addProblem("fix.missing.events");
102 }
103
104 if (referencePeriod == null) {
105 addProblem("fix.missing.reference.period");
106 }
107
108 if (analysisPeriods == null || analysisPeriods.length < 1) {
109 addProblem("fix.missing.analysis.periods");
110 }
111
112 if (qSectorStart == null) {
113 addProblem("fix.missing.qstart.sector");
114 }
115
116 if (qSectorEnd == null) {
117 addProblem("fix.missing.qend.sector");
118 }
119
120 if (preprocessing == null) {
121 addProblem("fix.missing.preprocessing");
122 }
123
124 if (!hasProblems()) {
125 this.river = river;
126 this.from = from;
127 this.to = to;
128 this.step = step;
129 this.function = function;
130 this.events = events;
131 this.referencePeriod = referencePeriod;
132 this.analysisPeriods = analysisPeriods;
133 this.qSectorStart = qSectorStart;
134 this.qSectorEnd = qSectorEnd;
135 this.preprocessing = preprocessing;
136 }
137 }
138
139 public CalculationResult calculate() {
140
141 boolean debug = log.isDebugEnabled();
142
143 FixingsOverview overview =
144 FixingsOverviewFactory.getOverview(river);
145
146 if (overview == null) {
147 addProblem("fix.no.overview.available");
148 }
149
150 Function func = FunctionFactory.getInstance()
151 .getFunction(function);
152
153 if (func == null) {
154 addProblem("fix.invalid.function.name");
155 }
156
157 if (hasProblems()) {
158 return new CalculationResult(this);
159 }
160
161 final List<Column> eventColumns = getEventColumns(overview);
162
163 if (eventColumns.size() < 2) {
164 addProblem("fix.too.less.data.columns");
165 return new CalculationResult(this);
166 }
167
168 double [] kms = DoubleUtil.explode(from, to, step / 1000.0);
169
170 final double [] qs = new double[eventColumns.size()];
171 final double [] ws = new double[qs.length];
172 final boolean [] interpolated = new boolean[ws.length];
173
174 Fitting.QWDFactory qwdFactory = new Fitting.QWDFactory() {
175 @Override
176 public QWD create(double q, double w) {
177 // Check all the event columns for close match
178 // and take the description and the date from meta.
179 for (int i = 0; i < qs.length; ++i) {
180 if (Math.abs(qs[i]-q) < EPSILON
181 && Math.abs(ws[i]-w) < EPSILON) {
182 Column column = eventColumns.get(i);
183 return new QWD(
184 qs[i], ws[i],
185 column.getDescription(),
186 column.getDate(),
187 interpolated[i],
188 0d);
189 }
190 }
191 log.warn("cannot find column for (" + q + ", " + w + ")");
192 return new QWD(q, w);
193 }
194 };
195
196 Fitting fitting = new Fitting(func, qwdFactory, preprocessing);
197
198 String [] parameterNames = func.getParameterNames();
199
200 Parameters results =
201 new Parameters(createColumnNames(parameterNames));
202
203 boolean invalid = false;
204
205 if (debug) {
206 log.debug("number of kms: " + kms.length);
207 }
208
209 KMIndex<QW []> outliers = new KMIndex<QW []>();
210 KMIndex<QWD []> referenced = new KMIndex<QWD []>(kms.length);
211
212 int kmIndex = results.columnIndex("km");
213 int chiSqrIndex = results.columnIndex("chi_sqr");
214 int maxQIndex = results.columnIndex("max_q");
215 int stdDevIndex = results.columnIndex("std-dev");
216 int [] parameterIndices = results.columnIndices(parameterNames);
217
218 int numFailed = 0;
219
220 for (int i = 0; i < kms.length; ++i) {
221 double km = kms[i];
222
223 // Fill Qs and Ws from event columns.
224 for (int j = 0; j < ws.length; ++j) {
225 interpolated[j] = !eventColumns.get(j).getQW(km, qs, ws, j);
226 }
227
228 fitting.reset();
229
230 if (!fitting.fit(qs, ws)) {
231 ++numFailed;
232 addProblem(km, "fix.fitting.failed");
233 continue;
234 }
235
236 referenced.add(km, fitting.referencedToArray());
237
238 if (fitting.hasOutliers()) {
239 outliers.add(km, fitting.outliersToArray());
240 }
241
242 int row = results.newRow();
243 double [] values = fitting.getParameters();
244
245 results.set(row, kmIndex, km);
246 results.set(row, chiSqrIndex, fitting.getChiSquare());
247 results.set(row, stdDevIndex, fitting.getStandardDeviation());
248 results.set(row, maxQIndex, fitting.getMaxQ());
249 invalid |= results.set(row, parameterIndices, values);
250
251 if (debug) {
252 log.debug("km: "+km+" " + toString(parameterNames, values));
253 }
254 }
255
256 if (debug) {
257 log.debug("success: " + (kms.length - numFailed));
258 log.debug("failed: " + numFailed);
259 }
260
261 if (invalid) {
262 addProblem("fix.invalid.values");
263 results.removeNaNs();
264 }
265
266 KMIndex<AnalysisPeriod []> analysisPeriods =
267 calculateAnalysisPeriods(func, results, overview);
268
269 outliers.sort();
270 referenced.sort();
271 analysisPeriods.sort();
272
273 FixResult fr = new FixResult(
274 results,
275 referenced, outliers,
276 analysisPeriods);
277
278 return new CalculationResult(fr, this);
279 }
280
281 protected String toString(String [] parameterNames, double [] values) {
282 StringBuilder sb = new StringBuilder();
283 for (int i = 0; i < parameterNames.length; ++i) {
284 if (i > 0) sb.append(", ");
285 sb.append(parameterNames[i]).append(": ").append(values[i]);
286 }
287 return sb.toString();
288 }
289
290 protected List<Column> getEventColumns(FixingsOverview overview) {
291
292 FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();
293
294 IdsFilter ids = new IdsFilter(events);
295 DateRangeFilter rdf = new DateRangeFilter(
296 referencePeriod.getFrom(),
297 referencePeriod.getTo());
298 Filter filter = new AndFilter().add(rdf).add(ids);
299
300 List<Fixing.Column> metas = overview.filter(null, filter);
301
302 List<Column> columns = new ArrayList<Column>(metas.size());
303
304 for (Fixing.Column meta: metas) {
305
306 FixingsColumn data = fcf.getColumnData(meta);
307 if (data == null) {
308 addProblem("fix.cannot.load.data");
309 }
310 else {
311 columns.add(new Column(meta, data));
312 }
313 }
314
315 return columns;
316 }
317
318
319 protected KMIndex<AnalysisPeriod []> calculateAnalysisPeriods(
320 Function function,
321 Parameters parameters,
322 FixingsOverview overview
323 ) {
324 Range range = new Range(from, to);
325
326 int kmIndex = parameters.columnIndex("km");
327 int maxQIndex = parameters.columnIndex("max_q");
328
329 ColumnCache cc = new ColumnCache();
330
331 double [] wq = new double[2];
332
333 int [] parameterIndices =
334 parameters.columnIndices(function.getParameterNames());
335
336 double [] parameterValues = new double[parameterIndices.length];
337
338 DateAverager dateAverager = new DateAverager();
339
340 KMIndex<AnalysisPeriod []> results =
341 new KMIndex<AnalysisPeriod []>(parameters.size());
342
343 IdsFilter idsFilter = new IdsFilter(events);
344
345 for (int row = 0, R = parameters.size(); row < R; ++row) {
346 double km = parameters.get(row, kmIndex);
347 parameters.get(row, parameterIndices, parameterValues);
348
349 // This is the paraterized function for a given km.
350 de.intevation.flys.artifacts.math.Function instance =
351 function.instantiate(parameterValues);
352
353 KmFilter kmFilter = new KmFilter(km);
354
355 ArrayList<AnalysisPeriod> periodResults =
356 new ArrayList<AnalysisPeriod>(analysisPeriods.length);
357
358 for (DateRange analysisPeriod: analysisPeriods) {
359 DateRangeFilter drf = new DateRangeFilter(
360 analysisPeriod.getFrom(),
361 analysisPeriod.getTo());
362
363 QWD [] qSectorAverages = new QWD[4];
364 double [] qSectorStdDevs = new double[4];
365
366 ArrayList<QWD> allQWDs = new ArrayList<QWD>();
367
368 // for all Q sectors.
369 for (int qSector = qSectorStart; qSector < qSectorEnd; ++qSector) {
370
371 Filter filter = new AndFilter()
372 .add(kmFilter)
373 .add(new SectorFilter(qSector))
374 .add(drf)
375 .add(idsFilter);
376
377 List<Fixing.Column> metas = overview.filter(range, filter);
378
379 if (metas.isEmpty()) {
380 // No fixings for km and analysis period
381 continue;
382 }
383
384 double sumQ = 0.0;
385 double sumW = 0.0;
386
387 StandardDeviation stdDev = new StandardDeviation();
388
389 List<QWD> qwds = new ArrayList<QWD>(metas.size());
390
391 dateAverager.clear();
392
393 for (Fixing.Column meta: metas) {
394 if (meta.findQSector(km) != qSector) {
395 // Ignore not matching sectors.
396 continue;
397 }
398
399 Column column = cc.getColumn(meta);
400 if (column == null || !column.getQW(km, wq)) {
401 continue;
402 }
403
404 double fw = instance.value(wq[1]);
405 if (Double.isNaN(fw)) {
406 continue;
407 }
408
409 double dw = (wq[0] - fw)*100.0;
410
411 stdDev.increment(dw);
412
413 Date date = column.getDate();
414 String description = column.getDescription();
415
416 QWD qwd = new QWD(
417 wq[1], wq[0], description, date, true, dw);
418
419 qwds.add(qwd);
420
421 sumW += wq[0];
422 sumQ += wq[1];
423
424 dateAverager.add(date);
425 }
426
427 // Calulate average per Q sector.
428 int N = qwds.size();
429 if (N > 0) {
430 allQWDs.addAll(qwds);
431 double avgW = sumW / N;
432 double avgQ = sumQ / N;
433
434 double avgFw = instance.value(avgQ);
435 if (!Double.isNaN(avgFw)) {
436 double avgDw = (avgW - avgFw)*100.0;
437 Date avgDate = dateAverager.getAverage();
438
439 String avgDescription = "avg.deltawt." + qSector;
440
441 QWD avgQWD = new QWD(
442 avgQ, avgW, avgDescription, avgDate, true, avgDw);
443
444 qSectorAverages[qSector] = avgQWD;
445 }
446 qSectorStdDevs[qSector] = stdDev.getResult();
447 }
448 else {
449 qSectorStdDevs[qSector] = Double.NaN;
450 }
451 } // for all Q sectors
452
453 QWD [] aqwds = allQWDs.toArray(new QWD[allQWDs.size()]);
454
455 AnalysisPeriod periodResult = new AnalysisPeriod(
456 analysisPeriod,
457 aqwds,
458 qSectorAverages,
459 qSectorStdDevs);
460 periodResults.add(periodResult);
461 }
462
463 double maxQ = -Double.MAX_VALUE;
464 for (AnalysisPeriod ap: periodResults) {
465 double q = ap.getMaxQ();
466 if (q > maxQ) {
467 maxQ = q;
468 }
469 }
470
471 double oldMaxQ = parameters.get(row, maxQIndex);
472 if (oldMaxQ < maxQ) {
473 parameters.set(row, maxQIndex, maxQ);
474 }
475
476 results.add(km, periodResults.toArray(
477 new AnalysisPeriod[periodResults.size()]));
478 }
479
480 return results;
481 }
482
483 /** Helper class to bundle the meta information of a column
484 * and the real data.
485 */
486 protected static class Column {
487
488 protected Fixing.Column meta;
489 protected FixingsColumn data;
490
491 public Column() {
492 }
493
494 public Column(Fixing.Column meta, FixingsColumn data) {
495 this.meta = meta;
496 this.data = data;
497 }
498
499 public Date getDate() {
500 return meta.getStartTime();
501 }
502
503 public String getDescription() {
504 return meta.getDescription();
505 }
506
507 public boolean getQW(
508 double km,
509 double [] qs,
510 double [] ws,
511 int index
512 ) {
513 qs[index] = data.getQ(km);
514 return data.getW(km, ws, index);
515 }
516
517 public boolean getQW(double km, double [] wq) {
518 data.getW(km, wq, 0);
519 if (Double.isNaN(wq[0])) return false;
520 wq[1] = data.getQ(km);
521 return !Double.isNaN(wq[1]);
522 }
523 } // class Column
524
525
526 /**
527 * Helper class to find the data belonging to meta info more quickly.
528 */
529 public static class ColumnCache {
530
531 protected Map<Integer, Column> columns;
532
533 public ColumnCache() {
534 columns = new HashMap<Integer, Column>();
535 }
536
537 public Column getColumn(Fixing.Column meta) {
538 Integer key = meta.getId();
539 Column column = columns.get(key);
540 if (column == null) {
541 FixingsColumn data = FixingsColumnFactory
542 .getInstance()
543 .getColumnData(meta);
544 if (data != null) {
545 column = new Column(meta, data);
546 columns.put(key, column);
547 }
548 }
549 return column;
550 }
551 } // class ColumnCache
552
553 /** Fetch meta and data columns for analysis periods. */
554 protected Column [][] getAnalysisColumns(FixingsOverview overview) {
555
556 boolean debug = log.isDebugEnabled();
557 if (debug) {
558 log.debug("number analysis periods: " + analysisPeriods.length);
559 }
560
561 Column columns [][] = new Column[analysisPeriods.length][];
562
563 Range range = new Range(from, to);
564 SectorRangeFilter sectorRangeFilter =
565 new SectorRangeFilter(qSectorStart, qSectorEnd);
566
567 FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();
568
569 for (int i = 0; i < columns.length; ++i) {
570
571 // Construct filter for period.
572 DateRange period = analysisPeriods[i];
573
574 AndFilter filter = new AndFilter();
575
576 DateRangeFilter dateRangeFilter =
577 new DateRangeFilter(period.getFrom(), period.getTo());
578
579 filter.add(dateRangeFilter);
580 filter.add(sectorRangeFilter);
581
582 List<Fixing.Column> metaCols = overview.filter(range, filter);
583
584 if (debug) {
585 log.debug("number of filtered columns: " + metaCols.size());
586 }
587
588 ArrayList<Column> cols = new ArrayList<Column>(metaCols.size());
589
590 // Only use columns which have data.
591 for (Fixing.Column meta: metaCols) {
592 FixingsColumn data = fcf.getColumnData(meta);
593 if (data != null) {
594 cols.add(new Column(meta, data));
595 }
596 }
597
598 if (debug) {
599 log.debug("failed loading: " + (metaCols.size()-cols.size()));
600 }
601 columns[i] = cols.toArray(new Column[cols.size()]);
602 }
603
604 return columns;
605 }
606
607 protected static String [] createColumnNames(String [] parameters) {
608 String [] result = new String[parameters.length + 4];
609 result[0] = "km";
610 result[1] = "chi_sqr";
611 result[2] = "max_q";
612 result[3] = "std-dev";
613 System.arraycopy(parameters, 0, result, 4, parameters.length);
614 return result;
615 }
616 }
617 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org