Mercurial > dive4elements > river
diff flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java @ 3435:262e7d7e58fe
FixA: Made curve fitting over the given calculation range reusable. Removed dead code.
flys-artifacts/trunk@5098 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Sun, 22 Jul 2012 10:38:30 +0000 |
parents | 1a636be7612b |
children | e111902834d3 |
line wrap: on
line diff
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java Sun Jul 22 09:49:56 2012 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java Sun Jul 22 10:38:30 2012 +0000 @@ -1,7 +1,11 @@ package de.intevation.flys.artifacts.model.fixings; +import de.intevation.artifacts.common.utils.StringUtils; + import de.intevation.flys.artifacts.access.FixAccess; +import de.intevation.flys.artifacts.math.fitting.Function; + import de.intevation.flys.artifacts.model.Calculation; import de.intevation.flys.artifacts.model.FixingsColumn; import de.intevation.flys.artifacts.model.FixingsColumnFactory; @@ -12,6 +16,10 @@ import de.intevation.flys.artifacts.model.FixingsOverview.IdsFilter; import de.intevation.flys.artifacts.model.FixingsOverview; +import de.intevation.flys.artifacts.model.Parameters; + +import de.intevation.flys.utils.DoubleUtil; +import de.intevation.flys.utils.KMIndex; import java.util.ArrayList; import java.util.Date; @@ -26,6 +34,44 @@ { private static Logger log = Logger.getLogger(FixCalculation.class); + public static final double EPSILON = 1e-4; + + public static final String [] STANDARD_COLUMNS = { + "km", "chi_sqr", "max_q", "std-dev" + }; + + protected static class FitResult { + + protected Parameters parameters; + protected KMIndex<QWD []> referenced; + protected KMIndex<QW []> outliers; + + public FitResult() { + } + + public FitResult( + Parameters parameters, + KMIndex<QWD []> referenced, + KMIndex<QW []> outliers + ) { + this.parameters = parameters; + this.referenced = referenced; + this.outliers = outliers; + } + + public Parameters getParameters() { + return parameters; + } + + public KMIndex<QWD []> getReferenced() { + return referenced; + } + + public KMIndex<QW []> getOutliers() { + return outliers; + } + } // class FitResult + /** Helper class to bundle the meta information of a column * and the real data. */ @@ -208,5 +254,124 @@ return columns; } + + protected FitResult doFitting(FixingsOverview overview, Function func) { + + boolean debug = log.isDebugEnabled(); + + final List<Column> eventColumns = getEventColumns(overview); + + if (eventColumns.size() < 2) { + addProblem("fix.too.less.data.columns"); + return null; + } + + final double [] qs = new double[eventColumns.size()]; + final double [] ws = new double[qs.length]; + final boolean [] interpolated = new boolean[ws.length]; + + Fitting.QWDFactory qwdFactory = new Fitting.QWDFactory() { + @Override + public QWD create(double q, double w) { + // Check all the event columns for close match + // and take the description and the date from meta. + for (int i = 0; i < qs.length; ++i) { + if (Math.abs(qs[i]-q) < EPSILON + && Math.abs(ws[i]-w) < EPSILON) { + Column column = eventColumns.get(i); + return new QWD( + qs[i], ws[i], + column.getDescription(), + column.getDate(), + interpolated[i], + 0d); + } + } + log.warn("cannot find column for (" + q + ", " + w + ")"); + return new QWD(q, w); + } + }; + + Fitting fitting = new Fitting(func, qwdFactory, preprocessing); + + String [] parameterNames = func.getParameterNames(); + + Parameters results = + new Parameters( + StringUtils.join(STANDARD_COLUMNS, parameterNames)); + + boolean invalid = false; + + double [] kms = DoubleUtil.explode(from, to, step / 1000.0); + + if (debug) { + log.debug("number of kms: " + kms.length); + } + + KMIndex<QW []> outliers = new KMIndex<QW []>(); + KMIndex<QWD []> referenced = new KMIndex<QWD []>(kms.length); + + int kmIndex = results.columnIndex("km"); + int chiSqrIndex = results.columnIndex("chi_sqr"); + int maxQIndex = results.columnIndex("max_q"); + int stdDevIndex = results.columnIndex("std-dev"); + int [] parameterIndices = results.columnIndices(parameterNames); + + int numFailed = 0; + + for (int i = 0; i < kms.length; ++i) { + double km = kms[i]; + + // Fill Qs and Ws from event columns. + for (int j = 0; j < ws.length; ++j) { + interpolated[j] = !eventColumns.get(j).getQW(km, qs, ws, j); + } + + fitting.reset(); + + if (!fitting.fit(qs, ws)) { + ++numFailed; + addProblem(km, "fix.fitting.failed"); + continue; + } + + referenced.add(km, fitting.referencedToArray()); + + if (fitting.hasOutliers()) { + outliers.add(km, fitting.outliersToArray()); + } + + int row = results.newRow(); + double [] values = fitting.getParameters(); + + results.set(row, kmIndex, km); + results.set(row, chiSqrIndex, fitting.getChiSquare()); + results.set(row, stdDevIndex, fitting.getStandardDeviation()); + results.set(row, maxQIndex, fitting.getMaxQ()); + invalid |= results.set(row, parameterIndices, values); + + if (debug) { + log.debug("km: "+km+" " + toString(parameterNames, values)); + } + } + + if (debug) { + log.debug("success: " + (kms.length - numFailed)); + log.debug("failed: " + numFailed); + } + + if (invalid) { + addProblem("fix.invalid.values"); + results.removeNaNs(); + } + + outliers.sort(); + referenced.sort(); + + return new FitResult( + results, + referenced, + outliers); + } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :