Mercurial > dive4elements > river
comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixAnalysisCalculation.java @ 3411:0ef83077c93f
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 : |