sascha@2729: package de.intevation.flys.artifacts.model.fixings; sascha@2729: sascha@2729: import de.intevation.flys.artifacts.FixationArtifactAccess; sascha@2729: sascha@2729: import de.intevation.flys.artifacts.math.fitting.Function; sascha@2729: import de.intevation.flys.artifacts.math.fitting.FunctionFactory; sascha@2729: sascha@2729: import de.intevation.flys.artifacts.model.Calculation; sascha@2729: import de.intevation.flys.artifacts.model.CalculationResult; sascha@2729: import de.intevation.flys.artifacts.model.FixingsColumn; sascha@2729: import de.intevation.flys.artifacts.model.FixingsColumnFactory; sascha@2729: import de.intevation.flys.artifacts.model.FixingsOverview; sascha@2729: import de.intevation.flys.artifacts.model.FixingsOverviewFactory; sascha@2729: import de.intevation.flys.artifacts.model.Parameters; sascha@2729: sascha@2744: import de.intevation.flys.artifacts.model.FixingsOverview.Range; sascha@2729: import de.intevation.flys.artifacts.model.FixingsOverview.Fixing; sascha@2729: import de.intevation.flys.artifacts.model.FixingsOverview.IdFilter; sascha@2744: import de.intevation.flys.artifacts.model.FixingsOverview.AndFilter; sascha@2744: import de.intevation.flys.artifacts.model.FixingsOverview.DateRangeFilter; sascha@2744: import de.intevation.flys.artifacts.model.FixingsOverview.SectorRangeFilter; sascha@2729: sascha@2729: import de.intevation.flys.utils.DoubleUtil; sascha@2729: sascha@2729: import java.util.ArrayList; sascha@2729: import java.util.List; sascha@2744: import java.util.Date; sascha@2729: sascha@2729: import org.apache.commons.math.MathException; sascha@2729: sascha@2729: import org.apache.commons.math.optimization.fitting.CurveFitter; sascha@2729: sascha@2729: import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer; sascha@2729: sascha@2729: import org.apache.log4j.Logger; sascha@2729: sascha@2729: public class FixCalculation sascha@2729: extends Calculation sascha@2729: { sascha@2786: private static Logger log = Logger.getLogger(FixCalculation.class); sascha@2729: sascha@2744: protected String river; sascha@2744: protected double from; sascha@2744: protected double to; sascha@2744: protected double step; sascha@2744: protected boolean preprocessing; sascha@2744: protected String function; sascha@2744: protected int [] events; sascha@2744: protected long [][] analysisPeriods; sascha@2744: protected int qSectorStart; sascha@2744: protected int qSectorEnd; sascha@2729: sascha@2729: public FixCalculation() { sascha@2729: } sascha@2729: sascha@2729: public FixCalculation(FixationArtifactAccess access) { sascha@2729: sascha@2744: String river = access.getRiver(); sascha@2744: Double from = access.getFrom(); sascha@2744: Double to = access.getTo(); sascha@2744: Double step = access.getStep(); sascha@2744: String function = access.getFunction(); sascha@2744: int [] events = access.getEvents(); sascha@2744: long [][] analysisPeriods = access.getAnalysisPeriods(); sascha@2744: Integer qSectorStart = access.getQSectorStart(); sascha@2744: Integer qSectorEnd = access.getQSectorEnd(); sascha@2729: sascha@2729: if (river == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.river"); sascha@2729: } sascha@2729: sascha@2729: if (from == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.from"); sascha@2729: } sascha@2729: sascha@2729: if (to == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.to"); sascha@2729: } sascha@2729: sascha@2729: if (step == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.step"); sascha@2729: } sascha@2729: sascha@2729: if (function == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.function"); sascha@2729: } sascha@2729: sascha@2729: if (events == null || events.length < 1) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.events"); sascha@2729: } sascha@2729: sascha@2744: if (analysisPeriods == null || analysisPeriods.length < 1) { sascha@2744: // TODO: i18n sascha@2744: addProblem("fix.missing.analysis.periods"); sascha@2744: } sascha@2744: sascha@2744: if (qSectorStart == null) { sascha@2744: // TODO: i18n sascha@2744: addProblem("fix.missing.qstart.sector"); sascha@2744: } sascha@2744: sascha@2744: if (qSectorEnd == null) { sascha@2744: // TODO: i18n sascha@2744: addProblem("fix.missing.qend.sector"); sascha@2744: } sascha@2744: sascha@2729: if (!hasProblems()) { sascha@2744: this.river = river; sascha@2744: this.from = from; sascha@2744: this.to = to; sascha@2744: this.step = step; sascha@2744: this.function = function; sascha@2744: this.events = events; sascha@2744: this.analysisPeriods = analysisPeriods; sascha@2744: this.qSectorStart = qSectorStart; sascha@2744: this.qSectorEnd = qSectorEnd; sascha@2729: } sascha@2729: } sascha@2729: sascha@2729: public CalculationResult calculate() { sascha@2729: sascha@2729: FixingsOverview overview = sascha@2729: FixingsOverviewFactory.getOverview(river); sascha@2729: sascha@2729: if (overview == null) { sascha@2729: addProblem("fix.no.overview.available"); sascha@2729: } sascha@2729: sascha@2729: Function func = FunctionFactory.getInstance() sascha@2729: .getFunction(function); sascha@2729: sascha@2729: if (func == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.invalid.function.name"); sascha@2729: } sascha@2729: sascha@2729: if (hasProblems()) { sascha@2729: return new CalculationResult(this); sascha@2729: } sascha@2729: sascha@2729: FixingsColumnFactory fcf = FixingsColumnFactory.getInstance(); sascha@2729: sascha@2729: List dataColumns = sascha@2729: new ArrayList(events.length); sascha@2729: sascha@2729: for (int eventId: events) { sascha@2729: IdFilter idFilter = new IdFilter(eventId); sascha@2729: sascha@2729: List columns = overview.filter(null, idFilter); sascha@2729: if (columns.isEmpty()) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.missing.column", eventId); sascha@2729: continue; sascha@2729: } sascha@2729: FixingsColumn dataColumn = fcf.getColumnData(columns.get(0)); sascha@2729: if (dataColumn == null) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.cannot.load.data", eventId); sascha@2729: continue; sascha@2729: } sascha@2729: dataColumns.add(dataColumn); sascha@2729: } sascha@2729: sascha@2729: if (dataColumns.size() < 2) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.too.less.data.columns"); sascha@2729: return new CalculationResult(this); sascha@2729: } sascha@2729: sascha@2729: double [] kms = DoubleUtil.explode(from, to, step); sascha@2729: sascha@2729: double [] ws = new double[dataColumns.size()]; sascha@2729: double [] qs = new double[ws.length]; sascha@2729: sascha@2729: String [] parameterNames = func.getParameterNames(); sascha@2729: sascha@2729: Parameters results = sascha@2729: new Parameters(createColumnNames(parameterNames)); sascha@2729: sascha@2729: boolean invalid = false; sascha@2729: sascha@2729: for (int i = 0; i < kms.length; ++i) { sascha@2729: double km = kms[i]; sascha@2729: sascha@2729: for (int j = 0; j < ws.length; ++j) { sascha@2729: FixingsColumn column = dataColumns.get(j); sascha@2729: qs[j] = column.getQ(km); sascha@2729: boolean interpolated = column.getW(km, ws, j); sascha@2729: // TODO: mark as interpolated. sascha@2729: } sascha@2729: sascha@2729: // TODO: Do preprocessing here! sascha@2729: double [] parameters = fit(func, km, ws, qs); sascha@2729: if (parameters == null) { // Problems are reported already. sascha@2729: continue; sascha@2729: } sascha@2729: sascha@2729: int row = results.newRow(); sascha@2729: sascha@2729: results.set(row, "km", km); sascha@2729: for (int j = 0; j < parameters.length; ++j) { sascha@2729: if (Double.isNaN(parameters[j])) { sascha@2729: invalid = true; sascha@2729: } sascha@2729: else { sascha@2729: results.set(row, parameterNames[j], parameters[j]); sascha@2729: } sascha@2729: } sascha@2729: // TODO: Calculate statistics, too! sascha@2729: } sascha@2729: sascha@2729: if (invalid) { sascha@2729: // TODO: i18n sascha@2729: addProblem("fix.invalid.values"); sascha@2729: results.removeNaNs(); sascha@2729: } sascha@2729: sascha@2744: // Calculate Delta W/t sascha@2783: DeltaWTsKM deltaWTsKM = calculateDeltaWTs( sascha@2744: func, sascha@2744: overview, sascha@2744: results); sascha@2744: sascha@2786: FixResult fr = new FixResult(results, deltaWTsKM); sascha@2786: sascha@2786: return new CalculationResult(fr, this); sascha@2744: } sascha@2744: sascha@2783: public DeltaWTsKM calculateDeltaWTs( sascha@2744: Function function, sascha@2744: FixingsOverview overview, sascha@2744: Parameters results sascha@2744: ) { sascha@2783: DeltaWTsKM deltaWTsKM = new DeltaWTsKM(results.size()); sascha@2744: sascha@2744: Column [][] analysisColumns = getAnalysisColumns(overview); sascha@2744: sascha@2744: int [] parameterIndices = sascha@2744: results.columnIndices(function.getParameterNames()); sascha@2744: sascha@2744: double [] parameterValues = sascha@2744: new double[parameterIndices.length]; sascha@2744: sascha@2744: double [] ow = new double[1]; sascha@2744: sascha@2744: for (int i = 0, N = results.size(); i < N; ++i) { sascha@2744: double km = results.get(i, "km"); sascha@2744: results.get(i, parameterIndices, parameterValues); sascha@2744: sascha@2783: DeltaWTsKM.KM dwtkm = new DeltaWTsKM.KM(km); sascha@2783: deltaWTsKM.add(dwtkm); sascha@2783: sascha@2744: // This is the paraterized function for a given km. sascha@2744: de.intevation.flys.artifacts.math.Function instance = sascha@2744: function.instantiate(parameterValues); sascha@2744: sascha@2744: // Evaluate all columns for all analysis periods. sascha@2744: for (int j = 0; j < analysisColumns.length; ++j) { sascha@2744: Column [] periodColumns = analysisColumns[j]; sascha@2744: sascha@2744: for (int k = 0; k < periodColumns.length; ++k) { sascha@2744: Column pc = periodColumns[k]; sascha@2744: sascha@2744: // Q from real data. sascha@2744: double q = pc.data.getQ(km); sascha@2744: if (Double.isNaN(q)) { sascha@2744: continue; sascha@2744: } sascha@2744: sascha@2744: // Calculate W from function. sascha@2744: double nw = instance.value(q); sascha@2744: if (Double.isNaN(nw)) { sascha@2744: continue; sascha@2744: } sascha@2744: sascha@2744: // W from real data. sascha@2744: pc.data.getW(km, ow); sascha@2744: sascha@2744: if (Double.isNaN(ow[0])) { sascha@2744: continue; sascha@2744: } sascha@2744: sascha@2744: double deltaW = ow[0] - nw; sascha@2744: sascha@2744: DeltaWT deltaWT = new DeltaWT( sascha@2744: deltaW, sascha@2744: pc.meta.getStartTime(), sascha@2744: pc.meta.getDescription()); sascha@2744: sascha@2783: dwtkm.add(deltaWT); sascha@2744: } sascha@2744: } sascha@2744: } sascha@2744: sascha@2783: return deltaWTsKM; sascha@2744: } sascha@2744: sascha@2744: /** Helper class to bundle the meta information of a column sascha@2744: * and the real data. sascha@2744: */ sascha@2744: protected static class Column { sascha@2744: sascha@2744: protected Fixing.Column meta; sascha@2744: protected FixingsColumn data; sascha@2744: sascha@2744: public Column() { sascha@2744: } sascha@2744: sascha@2744: public Column(Fixing.Column meta, FixingsColumn data) { sascha@2744: this.meta = meta; sascha@2744: this.data = data; sascha@2744: } sascha@2744: } // class Column sascha@2744: sascha@2744: /** Fetch meta and data columns for analysis periods. */ sascha@2744: protected Column [][] getAnalysisColumns(FixingsOverview overview) { sascha@2744: sascha@2744: Column columns [][] = new Column[analysisPeriods.length][]; sascha@2744: sascha@2744: Range range = new Range(from, to); sascha@2744: SectorRangeFilter sectorRangeFilter = sascha@2744: new SectorRangeFilter(qSectorStart, qSectorEnd); sascha@2744: sascha@2744: FixingsColumnFactory fcf = FixingsColumnFactory.getInstance(); sascha@2744: sascha@2744: for (int i = 0; i < columns.length; ++i) { sascha@2744: sascha@2744: // Construct filter for period. sascha@2744: long [] period = analysisPeriods[i]; sascha@2744: sascha@2744: AndFilter filter = new AndFilter(); sascha@2744: sascha@2744: DateRangeFilter dateRangeFilter = sascha@2744: new DateRangeFilter( sascha@2744: new Date(period[0]), sascha@2744: new Date(period[1])); sascha@2744: sascha@2744: filter.add(dateRangeFilter); sascha@2744: filter.add(sectorRangeFilter); sascha@2744: sascha@2744: List metaCols = overview.filter(range, filter); sascha@2744: sascha@2744: ArrayList cols = new ArrayList(metaCols.size()); sascha@2744: sascha@2744: // Only use columns which have data. sascha@2744: for (Fixing.Column meta: metaCols) { sascha@2744: FixingsColumn data = fcf.getColumnData(meta); sascha@2744: if (data != null) { sascha@2744: cols.add(new Column(meta, data)); sascha@2744: } sascha@2744: } sascha@2744: columns[i] = cols.toArray(new Column[cols.size()]); sascha@2744: } sascha@2744: sascha@2744: return columns; sascha@2729: } sascha@2729: sascha@2729: protected static String [] createColumnNames(String [] parameters) { sascha@2729: String [] result = new String[parameters.length + 1]; sascha@2729: result[0] = "km"; sascha@2729: // TODO: Add statistic columns, too. sascha@2729: System.arraycopy(parameters, 0, result, 1, parameters.length); sascha@2729: return result; sascha@2729: } sascha@2729: sascha@2729: protected double [] fit( sascha@2729: Function function, sascha@2729: double km, sascha@2729: double [] ws, sascha@2729: double [] qs sascha@2729: ) { sascha@2729: LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer(); sascha@2729: CurveFitter cf = new CurveFitter(lmo); sascha@2729: sascha@2729: boolean missingWs = false; sascha@2729: boolean missingQs = false; sascha@2729: sascha@2729: for (int i = 0; i < ws.length; ++i) { sascha@2729: boolean ignore = false; sascha@2729: if (Double.isNaN(ws[i])) { sascha@2729: ignore = true; sascha@2729: if (!missingWs) { sascha@2729: missingWs = true; sascha@2729: // TODO: i18n sascha@2729: addProblem(km, "fix.missing.w"); sascha@2729: } sascha@2729: } sascha@2729: if (Double.isNaN(qs[i])) { sascha@2729: ignore = true; sascha@2729: if (!missingQs) { sascha@2729: missingQs = true; sascha@2729: // TODO: i18n sascha@2729: addProblem(km, "fix.missing.q"); sascha@2729: } sascha@2729: } sascha@2729: if (!ignore) { sascha@2729: cf.addObservedPoint(ws[i], qs[i]); sascha@2729: } sascha@2729: } sascha@2729: sascha@2729: try { sascha@2729: return cf.fit(function, function.getInitialGuess()); sascha@2729: } sascha@2729: catch (MathException me) { sascha@2729: addProblem(km, "fix.fitting.failed"); sascha@2729: } sascha@2729: sascha@2729: return null; sascha@2729: } sascha@2729: } sascha@2729: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :