Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/org/dive4elements/river/artifacts/model/fixings/FixCalculation.java @ 5831:bd047b71ab37

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Repaired internal references
author Sascha L. Teichmann <teichmann@intevation.de>
date Thu, 25 Apr 2013 12:06:39 +0200
parents flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java@6e6753b9757f
children
comparison
equal deleted inserted replaced
5830:160f53ee0870 5831:bd047b71ab37
1 package org.dive4elements.river.artifacts.model.fixings;
2
3 import org.dive4elements.artifacts.common.utils.StringUtils;
4
5 import org.dive4elements.river.artifacts.access.FixAccess;
6
7 import org.dive4elements.river.artifacts.math.fitting.Function;
8 import org.dive4elements.river.artifacts.math.fitting.FunctionFactory;
9
10 import org.dive4elements.river.artifacts.model.Calculation;
11 import org.dive4elements.river.artifacts.model.CalculationResult;
12 import org.dive4elements.river.artifacts.model.FixingsColumn;
13 import org.dive4elements.river.artifacts.model.FixingsColumnFactory;
14
15 import org.dive4elements.river.artifacts.model.FixingsOverview.Fixing.Filter;
16
17 import org.dive4elements.river.artifacts.model.FixingsOverview.Fixing;
18 import org.dive4elements.river.artifacts.model.FixingsOverview.IdsFilter;
19
20 import org.dive4elements.river.artifacts.model.FixingsOverview;
21 import org.dive4elements.river.artifacts.model.FixingsOverviewFactory;
22 import org.dive4elements.river.artifacts.model.Parameters;
23
24 import org.dive4elements.river.utils.DoubleUtil;
25 import org.dive4elements.river.utils.KMIndex;
26
27 import java.util.ArrayList;
28 import java.util.Date;
29 import java.util.HashMap;
30 import java.util.List;
31 import java.util.Map;
32
33 import org.apache.log4j.Logger;
34
35 /** Calculation base class for fix. */
36 public abstract class FixCalculation
37 extends Calculation
38 {
39 private static Logger log = Logger.getLogger(FixCalculation.class);
40
41 public static final double EPSILON = 1e-4;
42
43 public static final String [] STANDARD_COLUMNS = {
44 "km", "chi_sqr", "max_q", "std-dev"
45 };
46
47 protected static class FitResult {
48
49 protected Parameters parameters;
50 protected KMIndex<QWD []> referenced;
51 protected KMIndex<QWI []> outliers;
52
53 public FitResult() {
54 }
55
56 public FitResult(
57 Parameters parameters,
58 KMIndex<QWD []> referenced,
59 KMIndex<QWI []> outliers
60 ) {
61 this.parameters = parameters;
62 this.referenced = referenced;
63 this.outliers = outliers;
64 }
65
66 public Parameters getParameters() {
67 return parameters;
68 }
69
70 public KMIndex<QWD []> getReferenced() {
71 return referenced;
72 }
73
74 public KMIndex<QWI []> getOutliers() {
75 return outliers;
76 }
77 } // class FitResult
78
79 /** Helper class to bundle the meta information of a column
80 * and the real data.
81 */
82 protected static class Column {
83
84 protected Fixing.Column meta;
85 protected FixingsColumn data;
86 protected int index;
87
88 public Column() {
89 }
90
91 public Column(Fixing.Column meta, FixingsColumn data, int index) {
92 this.meta = meta;
93 this.data = data;
94 this.index = index;
95 }
96
97 public Date getDate() {
98 return meta.getStartTime();
99 }
100
101 public String getDescription() {
102 return meta.getDescription();
103 }
104
105 public int getIndex() {
106 return index;
107 }
108
109 public boolean getQW(
110 double km,
111 double [] qs,
112 double [] ws,
113 int index
114 ) {
115 qs[index] = data.getQ(km);
116 return data.getW(km, ws, index);
117 }
118
119 public boolean getQW(double km, double [] wq) {
120 data.getW(km, wq, 0);
121 if (Double.isNaN(wq[0])) return false;
122 wq[1] = data.getQ(km);
123 return !Double.isNaN(wq[1]);
124 }
125 } // class Column
126
127 /**
128 * Helper class to find the data belonging to meta info more quickly.
129 */
130 protected static class ColumnCache {
131
132 protected Map<Integer, Column> columns;
133
134 public ColumnCache() {
135 columns = new HashMap<Integer, Column>();
136 }
137
138 public Column getColumn(Fixing.Column meta) {
139 Integer key = meta.getId();
140 Column column = columns.get(key);
141 if (column == null) {
142 FixingsColumn data = FixingsColumnFactory
143 .getInstance()
144 .getColumnData(meta);
145 if (data != null) {
146 column = new Column(meta, data, columns.size());
147 columns.put(key, column);
148 }
149 }
150 return column;
151 }
152 } // class ColumnCache
153
154
155 protected String river;
156 protected double from;
157 protected double to;
158 protected double step;
159 protected boolean preprocessing;
160 protected String function;
161 protected int [] events;
162 protected int qSectorStart;
163 protected int qSectorEnd;
164
165 public FixCalculation() {
166 }
167
168 public FixCalculation(FixAccess access) {
169 String river = access.getRiver();
170 Double from = access.getFrom();
171 Double to = access.getTo();
172 Double step = access.getStep();
173 String function = access.getFunction();
174 int [] events = access.getEvents();
175 Integer qSectorStart = access.getQSectorStart();
176 Integer qSectorEnd = access.getQSectorEnd();
177 Boolean preprocessing = access.getPreprocessing();
178
179 if (river == null) {
180 addProblem("fix.missing.river");
181 }
182
183 if (from == null) {
184 addProblem("fix.missing.from");
185 }
186
187 if (to == null) {
188 addProblem("fix.missing.to");
189 }
190
191 if (step == null) {
192 addProblem("fix.missing.step");
193 }
194
195 if (function == null) {
196 addProblem("fix.missing.function");
197 }
198
199 if (events == null || events.length < 1) {
200 addProblem("fix.missing.events");
201 }
202
203 if (qSectorStart == null) {
204 addProblem("fix.missing.qstart.sector");
205 }
206
207 if (qSectorEnd == null) {
208 addProblem("fix.missing.qend.sector");
209 }
210
211 if (preprocessing == null) {
212 addProblem("fix.missing.preprocessing");
213 }
214
215 if (!hasProblems()) {
216 this.river = river;
217 this.from = from;
218 this.to = to;
219 this.step = step;
220 this.function = function;
221 this.events = events;
222 this.qSectorStart = qSectorStart;
223 this.qSectorEnd = qSectorEnd;
224 this.preprocessing = preprocessing;
225 }
226 }
227
228 protected static String toString(
229 String [] parameterNames,
230 double [] values
231 ) {
232 StringBuilder sb = new StringBuilder();
233 for (int i = 0; i < parameterNames.length; ++i) {
234 if (i > 0) sb.append(", ");
235 sb.append(parameterNames[i]).append(": ").append(values[i]);
236 }
237 return sb.toString();
238 }
239
240
241 /** Create filter to accept only the chosen events.
242 * This factored out out to be overwritten.
243 */
244 protected Filter createFilter() {
245 return new IdsFilter(events);
246 }
247
248 protected List<Column> getEventColumns(
249 FixingsOverview overview,
250 ColumnCache cc
251 ) {
252 FixingsColumnFactory fcf = FixingsColumnFactory.getInstance();
253
254 Filter filter = createFilter();
255
256 List<Fixing.Column> metas = overview.filter(null, filter);
257
258 List<Column> columns = new ArrayList<Column>(metas.size());
259
260 for (Fixing.Column meta: metas) {
261
262 Column data = cc.getColumn(meta);
263 if (data == null) {
264 addProblem("fix.cannot.load.data");
265 }
266 else {
267 columns.add(data);
268 }
269 }
270
271 return columns;
272 }
273
274 // Fit a function to the given points from fixation.
275 protected FitResult doFitting(
276 FixingsOverview overview,
277 ColumnCache cc,
278 Function func
279 ) {
280 boolean debug = log.isDebugEnabled();
281
282 final List<Column> eventColumns = getEventColumns(overview, cc);
283
284 if (eventColumns.size() < 2) {
285 addProblem("fix.too.less.data.columns");
286 return null;
287 }
288
289 final double [] qs = new double[eventColumns.size()];
290 final double [] ws = new double[qs.length];
291 final boolean [] interpolated = new boolean[ws.length];
292
293 Fitting.QWDFactory qwdFactory = new Fitting.QWDFactory() {
294 @Override
295 public QWD create(double q, double w) {
296 // Check all the event columns for close match
297 // and take the description and the date from meta.
298 for (int i = 0; i < qs.length; ++i) {
299 if (Math.abs(qs[i]-q) < EPSILON
300 && Math.abs(ws[i]-w) < EPSILON) {
301 Column column = eventColumns.get(i);
302 return new QWD(
303 qs[i], ws[i],
304 column.getDescription(),
305 column.getDate(),
306 interpolated[i],
307 0d,
308 column.getIndex());
309 }
310 }
311 log.warn("cannot find column for (" + q + ", " + w + ")");
312 return new QWD(q, w);
313 }
314 };
315
316 Fitting fitting = new Fitting(func, qwdFactory, preprocessing);
317
318 String [] parameterNames = func.getParameterNames();
319
320 Parameters results =
321 new Parameters(
322 StringUtils.join(STANDARD_COLUMNS, parameterNames));
323
324 boolean invalid = false;
325
326 double [] kms = DoubleUtil.explode(from, to, step / 1000.0);
327
328 if (debug) {
329 log.debug("number of kms: " + kms.length);
330 }
331
332 KMIndex<QWI []> outliers = new KMIndex<QWI []>();
333 KMIndex<QWD []> referenced = new KMIndex<QWD []>(kms.length);
334
335 int kmIndex = results.columnIndex("km");
336 int chiSqrIndex = results.columnIndex("chi_sqr");
337 int maxQIndex = results.columnIndex("max_q");
338 int stdDevIndex = results.columnIndex("std-dev");
339 int [] parameterIndices = results.columnIndices(parameterNames);
340
341 int numFailed = 0;
342
343 for (int i = 0; i < kms.length; ++i) {
344 double km = kms[i];
345
346 // Fill Qs and Ws from event columns.
347 for (int j = 0; j < ws.length; ++j) {
348 interpolated[j] = !eventColumns.get(j).getQW(km, qs, ws, j);
349 }
350
351 fitting.reset();
352
353 if (!fitting.fit(qs, ws)) {
354 ++numFailed;
355 addProblem(km, "fix.fitting.failed");
356 continue;
357 }
358
359 referenced.add(km, fitting.referencedToArray());
360
361 if (fitting.hasOutliers()) {
362 outliers.add(km, fitting.outliersToArray());
363 }
364
365 int row = results.newRow();
366 double [] values = fitting.getParameters();
367
368 results.set(row, kmIndex, km);
369 results.set(row, chiSqrIndex, fitting.getChiSquare());
370 results.set(row, stdDevIndex, fitting.getStandardDeviation());
371 results.set(row, maxQIndex, fitting.getMaxQ());
372 invalid |= results.set(row, parameterIndices, values);
373
374 if (debug) {
375 log.debug("km: "+km+" " + toString(parameterNames, values));
376 }
377 }
378
379 if (debug) {
380 log.debug("success: " + (kms.length - numFailed));
381 log.debug("failed: " + numFailed);
382 }
383
384 if (invalid) {
385 addProblem("fix.invalid.values");
386 results.removeNaNs();
387 }
388
389 outliers.sort();
390 referenced.sort();
391
392 return new FitResult(
393 results,
394 referenced,
395 outliers);
396 }
397
398 public CalculationResult calculate() {
399 FixingsOverview overview =
400 FixingsOverviewFactory.getOverview(river);
401
402 if (overview == null) {
403 addProblem("fix.no.overview.available");
404 }
405
406 Function func = FunctionFactory.getInstance()
407 .getFunction(function);
408
409 if (func == null) {
410 addProblem("fix.invalid.function.name");
411 }
412
413 if (hasProblems()) {
414 return new CalculationResult(this);
415 }
416
417 return innerCalculate(overview, func);
418 }
419
420 protected abstract CalculationResult innerCalculate(
421 FixingsOverview overview,
422 Function function
423 );
424 }
425 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org