Mercurial > dive4elements > river
view flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/fixings/FixCalculation.java @ 3215:750e98fc8b76
FixA: Tweaked the derivate diagram a bit and added chart info.
flys-artifacts/trunk@4838 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 29 Jun 2012 15:40:43 +0000 |
parents | ae14f412ba10 |
children | 50d61a2494cb |
line wrap: on
line source
package de.intevation.flys.artifacts.model.fixings; import de.intevation.flys.artifacts.access.FixationArtifactAccess; import de.intevation.flys.artifacts.math.fitting.Function; import de.intevation.flys.artifacts.math.fitting.FunctionFactory; import de.intevation.flys.artifacts.model.Calculation; import de.intevation.flys.artifacts.model.CalculationResult; import de.intevation.flys.artifacts.model.DateRange; import de.intevation.flys.artifacts.model.FixingsColumn; import de.intevation.flys.artifacts.model.FixingsColumnFactory; import de.intevation.flys.artifacts.model.FixingsOverview.AndFilter; import de.intevation.flys.artifacts.model.FixingsOverview.DateRangeFilter; import de.intevation.flys.artifacts.model.FixingsOverview.Fixing.Filter; import de.intevation.flys.artifacts.model.FixingsOverview.Fixing; import de.intevation.flys.artifacts.model.FixingsOverview.IdsFilter; import de.intevation.flys.artifacts.model.FixingsOverview.KmFilter; import de.intevation.flys.artifacts.model.FixingsOverview.SectorFilter; import de.intevation.flys.artifacts.model.FixingsOverview.SectorRangeFilter; import de.intevation.flys.artifacts.model.FixingsOverview; import de.intevation.flys.artifacts.model.FixingsOverviewFactory; import de.intevation.flys.artifacts.model.Parameters; import de.intevation.flys.artifacts.model.Range; import de.intevation.flys.utils.DateAverager; import de.intevation.flys.utils.DoubleUtil; import de.intevation.flys.utils.KMIndex; import java.util.ArrayList; import java.util.Date; import java.util.HashMap; import java.util.List; import java.util.Map; import org.apache.commons.math.stat.descriptive.moment.StandardDeviation; import org.apache.log4j.Logger; public class FixCalculation extends Calculation { private static Logger log = Logger.getLogger(FixCalculation.class); public static final double EPSILON = 1e-4; protected String river; protected double from; protected double to; protected double step; protected boolean preprocessing; protected String function; protected int [] events; protected DateRange referencePeriod; protected DateRange [] analysisPeriods; protected int qSectorStart; protected int qSectorEnd; public FixCalculation() { } public FixCalculation(FixationArtifactAccess access) { String river = access.getRiver(); Double from = access.getFrom(); Double to = access.getTo(); Double step = access.getStep(); String function = access.getFunction(); int [] events = access.getEvents(); DateRange referencePeriod = access.getReferencePeriod(); DateRange [] analysisPeriods = access.getAnalysisPeriods(); Integer qSectorStart = access.getQSectorStart(); Integer qSectorEnd = access.getQSectorEnd(); Boolean preprocessing = access.getPreprocessing(); if (river == null) { addProblem("fix.missing.river"); } if (from == null) { addProblem("fix.missing.from"); } if (to == null) { addProblem("fix.missing.to"); } if (step == null) { addProblem("fix.missing.step"); } if (function == null) { addProblem("fix.missing.function"); } if (events == null || events.length < 1) { addProblem("fix.missing.events"); } if (referencePeriod == null) { addProblem("fix.missing.reference.period"); } if (analysisPeriods == null || analysisPeriods.length < 1) { addProblem("fix.missing.analysis.periods"); } if (qSectorStart == null) { addProblem("fix.missing.qstart.sector"); } if (qSectorEnd == null) { addProblem("fix.missing.qend.sector"); } if (preprocessing == null) { addProblem("fix.missing.preprocessing"); } if (!hasProblems()) { this.river = river; this.from = from; this.to = to; this.step = step; this.function = function; this.events = events; this.referencePeriod = referencePeriod; this.analysisPeriods = analysisPeriods; this.qSectorStart = qSectorStart; this.qSectorEnd = qSectorEnd; this.preprocessing = preprocessing; } } public CalculationResult calculate() { boolean debug = log.isDebugEnabled(); FixingsOverview overview = FixingsOverviewFactory.getOverview(river); if (overview == null) { addProblem("fix.no.overview.available"); } Function func = FunctionFactory.getInstance() .getFunction(function); if (func == null) { addProblem("fix.invalid.function.name"); } if (hasProblems()) { return new CalculationResult(this); } final List<Column> eventColumns = getEventColumns(overview); if (eventColumns.size() < 2) { addProblem("fix.too.less.data.columns"); return new CalculationResult(this); } double [] kms = DoubleUtil.explode(from, to, step / 1000.0); 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(createColumnNames(parameterNames)); boolean invalid = false; 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(); } KMIndex<AnalysisPeriod []> analysisPeriods = calculateAnalysisPeriods(func, results, overview); outliers.sort(); referenced.sort(); analysisPeriods.sort(); FixResult fr = new FixResult( results, referenced, outliers, analysisPeriods); return new CalculationResult(fr, this); } protected String toString(String [] parameterNames, double [] values) { StringBuilder sb = new StringBuilder(); for (int i = 0; i < parameterNames.length; ++i) { if (i > 0) sb.append(", "); sb.append(parameterNames[i]).append(": ").append(values[i]); } return sb.toString(); } protected List<Column> getEventColumns(FixingsOverview overview) { FixingsColumnFactory fcf = FixingsColumnFactory.getInstance(); IdsFilter ids = new IdsFilter(events); DateRangeFilter rdf = new DateRangeFilter( referencePeriod.getFrom(), referencePeriod.getTo()); Filter filter = new AndFilter().add(rdf).add(ids); List<Fixing.Column> metas = overview.filter(null, filter); List<Column> columns = new ArrayList<Column>(metas.size()); for (Fixing.Column meta: metas) { FixingsColumn data = fcf.getColumnData(meta); if (data == null) { addProblem("fix.cannot.load.data"); } else { columns.add(new Column(meta, data)); } } return columns; } protected KMIndex<AnalysisPeriod []> calculateAnalysisPeriods( Function function, Parameters parameters, FixingsOverview overview ) { Range range = new Range(from, to); int kmIndex = parameters.columnIndex("km"); int maxQIndex = parameters.columnIndex("max_q"); ColumnCache cc = new ColumnCache(); double [] wq = new double[2]; int [] parameterIndices = parameters.columnIndices(function.getParameterNames()); double [] parameterValues = new double[parameterIndices.length]; DateAverager dateAverager = new DateAverager(); KMIndex<AnalysisPeriod []> results = new KMIndex<AnalysisPeriod []>(parameters.size()); IdsFilter idsFilter = new IdsFilter(events); for (int row = 0, R = parameters.size(); row < R; ++row) { double km = parameters.get(row, kmIndex); parameters.get(row, parameterIndices, parameterValues); // This is the paraterized function for a given km. de.intevation.flys.artifacts.math.Function instance = function.instantiate(parameterValues); KmFilter kmFilter = new KmFilter(km); ArrayList<AnalysisPeriod> periodResults = new ArrayList<AnalysisPeriod>(analysisPeriods.length); for (DateRange analysisPeriod: analysisPeriods) { DateRangeFilter drf = new DateRangeFilter( analysisPeriod.getFrom(), analysisPeriod.getTo()); QWD [] qSectorAverages = new QWD[4]; double [] qSectorStdDevs = new double[4]; ArrayList<QWD> allQWDs = new ArrayList<QWD>(); // for all Q sectors. for (int qSector = qSectorStart; qSector < qSectorEnd; ++qSector) { Filter filter = new AndFilter() .add(kmFilter) .add(new SectorFilter(qSector)) .add(drf) .add(idsFilter); List<Fixing.Column> metas = overview.filter(range, filter); if (metas.isEmpty()) { // No fixings for km and analysis period continue; } double sumQ = 0.0; double sumW = 0.0; StandardDeviation stdDev = new StandardDeviation(); List<QWD> qwds = new ArrayList<QWD>(metas.size()); dateAverager.clear(); for (Fixing.Column meta: metas) { if (meta.findQSector(km) != qSector) { // Ignore not matching sectors. continue; } Column column = cc.getColumn(meta); if (column == null || !column.getQW(km, wq)) { continue; } double fw = instance.value(wq[1]); if (Double.isNaN(fw)) { continue; } double dw = (wq[0] - fw)*100.0; stdDev.increment(dw); Date date = column.getDate(); String description = column.getDescription(); QWD qwd = new QWD( wq[1], wq[0], description, date, true, dw); qwds.add(qwd); sumW += wq[0]; sumQ += wq[1]; dateAverager.add(date); } // Calulate average per Q sector. int N = qwds.size(); if (N > 0) { allQWDs.addAll(qwds); double avgW = sumW / N; double avgQ = sumQ / N; double avgFw = instance.value(avgQ); if (!Double.isNaN(avgFw)) { double avgDw = (avgW - avgFw)*100.0; Date avgDate = dateAverager.getAverage(); String avgDescription = "avg.deltawt." + qSector; QWD avgQWD = new QWD( avgQ, avgW, avgDescription, avgDate, true, avgDw); qSectorAverages[qSector] = avgQWD; } qSectorStdDevs[qSector] = stdDev.getResult(); } else { qSectorStdDevs[qSector] = Double.NaN; } } // for all Q sectors QWD [] aqwds = allQWDs.toArray(new QWD[allQWDs.size()]); AnalysisPeriod periodResult = new AnalysisPeriod( analysisPeriod, aqwds, qSectorAverages, qSectorStdDevs); periodResults.add(periodResult); } double maxQ = -Double.MAX_VALUE; for (AnalysisPeriod ap: periodResults) { double q = ap.getMaxQ(); if (q > maxQ) { maxQ = q; } } double oldMaxQ = parameters.get(row, maxQIndex); if (oldMaxQ < maxQ) { parameters.set(row, maxQIndex, maxQ); } results.add(km, periodResults.toArray( new AnalysisPeriod[periodResults.size()])); } return results; } /** Helper class to bundle the meta information of a column * and the real data. */ protected static class Column { protected Fixing.Column meta; protected FixingsColumn data; public Column() { } public Column(Fixing.Column meta, FixingsColumn data) { this.meta = meta; this.data = data; } public Date getDate() { return meta.getStartTime(); } public String getDescription() { return meta.getDescription(); } public boolean getQW( double km, double [] qs, double [] ws, int index ) { qs[index] = data.getQ(km); return data.getW(km, ws, index); } public boolean getQW(double km, double [] wq) { data.getW(km, wq, 0); if (Double.isNaN(wq[0])) return false; wq[1] = data.getQ(km); return !Double.isNaN(wq[1]); } } // class Column /** * Helper class to find the data belonging to meta info more quickly. */ public static class ColumnCache { protected Map<Integer, Column> columns; public ColumnCache() { columns = new HashMap<Integer, Column>(); } public Column getColumn(Fixing.Column meta) { Integer key = meta.getId(); Column column = columns.get(key); if (column == null) { FixingsColumn data = FixingsColumnFactory .getInstance() .getColumnData(meta); if (data != null) { column = new Column(meta, data); columns.put(key, column); } } return column; } } // class ColumnCache /** Fetch meta and data columns for analysis periods. */ protected Column [][] getAnalysisColumns(FixingsOverview overview) { boolean debug = log.isDebugEnabled(); if (debug) { log.debug("number analysis periods: " + analysisPeriods.length); } Column columns [][] = new Column[analysisPeriods.length][]; Range range = new Range(from, to); SectorRangeFilter sectorRangeFilter = new SectorRangeFilter(qSectorStart, qSectorEnd); FixingsColumnFactory fcf = FixingsColumnFactory.getInstance(); for (int i = 0; i < columns.length; ++i) { // Construct filter for period. DateRange period = analysisPeriods[i]; AndFilter filter = new AndFilter(); DateRangeFilter dateRangeFilter = new DateRangeFilter(period.getFrom(), period.getTo()); filter.add(dateRangeFilter); filter.add(sectorRangeFilter); List<Fixing.Column> metaCols = overview.filter(range, filter); if (debug) { log.debug("number of filtered columns: " + metaCols.size()); } ArrayList<Column> cols = new ArrayList<Column>(metaCols.size()); // Only use columns which have data. for (Fixing.Column meta: metaCols) { FixingsColumn data = fcf.getColumnData(meta); if (data != null) { cols.add(new Column(meta, data)); } } if (debug) { log.debug("failed loading: " + (metaCols.size()-cols.size())); } columns[i] = cols.toArray(new Column[cols.size()]); } return columns; } protected static String [] createColumnNames(String [] parameters) { String [] result = new String[parameters.length + 4]; result[0] = "km"; result[1] = "chi_sqr"; result[2] = "max_q"; result[3] = "std-dev"; System.arraycopy(parameters, 0, result, 4, parameters.length); return result; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :