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@3008: sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.AndFilter; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.DateRangeFilter; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.Fixing; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.IdFilter; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.KmFilter; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.Range; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.SectorRangeFilter; sascha@3008: import de.intevation.flys.artifacts.model.FixingsOverview.SectorFilter; sascha@3008: 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@3008: import de.intevation.flys.utils.DateAverager; sascha@2729: import de.intevation.flys.utils.DoubleUtil; sascha@3010: import de.intevation.flys.utils.Pair; sascha@2729: sascha@2729: import java.util.ArrayList; sascha@3008: import java.util.Date; sascha@3008: import java.util.HashMap; sascha@2729: import java.util.List; sascha@3008: import java.util.Map; 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@3003: protected String river; sascha@3003: protected double from; sascha@3003: protected double to; sascha@3003: protected double step; sascha@3003: protected boolean preprocessing; sascha@3003: protected String function; sascha@3003: protected int [] events; sascha@3003: protected DateRange [] analysisPeriods; sascha@3003: protected int qSectorStart; sascha@3003: protected int qSectorEnd; sascha@2729: sascha@2729: public FixCalculation() { sascha@2729: } sascha@2729: sascha@2729: public FixCalculation(FixationArtifactAccess access) { sascha@2729: sascha@3003: String river = access.getRiver(); sascha@3003: Double from = access.getFrom(); sascha@3003: Double to = access.getTo(); sascha@3003: Double step = access.getStep(); sascha@3003: String function = access.getFunction(); sascha@3003: int [] events = access.getEvents(); sascha@3003: DateRange [] analysisPeriods = access.getAnalysisPeriods(); sascha@3003: Integer qSectorStart = access.getQSectorStart(); sascha@3003: 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@2992: boolean debug = log.isDebugEnabled(); sascha@2992: 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@2992: double [] kms = DoubleUtil.explode(from, to, step / 1000.0); 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@2992: if (debug) { sascha@2992: log.debug("number of kms: " + kms.length); sascha@2992: } sascha@2992: sascha@3010: int kmIndex = results.columnIndex("km"); sascha@3010: int chiSqrIndex = results.columnIndex("chi_sqr"); sascha@3010: sascha@2992: int numFailed = 0; sascha@2992: 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@3010: Pair fitResult = fit(func, km, ws, qs); sascha@3010: sascha@2729: // TODO: Do preprocessing here! sascha@3010: if (fitResult == null) { // Problems are reported already. sascha@2992: ++numFailed; sascha@2729: continue; sascha@2729: } sascha@2729: sascha@2729: int row = results.newRow(); sascha@2729: sascha@3010: results.set(row, kmIndex, km); sascha@3010: results.set(row, chiSqrIndex, fitResult.getB()); sascha@3010: sascha@3010: double [] parameters = fitResult.getA(); 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: } sascha@2729: sascha@2992: if (debug) { sascha@2992: log.debug("success: " + (kms.length - numFailed)); sascha@2992: log.debug("failed: " + numFailed); sascha@2992: } sascha@2992: 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@2993: boolean debug = log.isDebugEnabled(); sascha@2993: 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@2993: int kmIdx = results.columnIndex("km"); sascha@2993: sascha@2744: for (int i = 0, N = results.size(); i < N; ++i) { sascha@2993: double km = results.get(i, kmIdx); 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@2993: int failedQ = 0; sascha@2993: int failedW = 0; sascha@2993: int failedC = 0; sascha@2993: 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@2993: ++failedQ; 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@2993: ++failedC; 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@2993: ++failedW; sascha@2744: continue; sascha@2744: } sascha@2744: sascha@2994: double deltaW = (ow[0] - nw)*100.0; // in cm 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@2993: if (debug) { sascha@2993: log.debug("failed W: " + failedW); sascha@2993: log.debug("failed Q: " + failedQ); sascha@2993: log.debug("failed C: " + failedC); sascha@2993: log.debug("input size: " + periodColumns.length); sascha@2993: log.debug("outpt size: " + dwtkm.size()); sascha@2993: } sascha@2744: } sascha@2744: } sascha@2744: sascha@2783: return deltaWTsKM; sascha@2744: } sascha@2744: sascha@3008: protected List calculateAnalysisPeriods( sascha@3008: Function function, sascha@3008: Parameters parameters, sascha@3008: FixingsOverview overview sascha@3008: ) { sascha@3008: ArrayList results sascha@3008: = new ArrayList(parameters.size()); sascha@3008: sascha@3008: Range range = new Range(from, to); sascha@3008: sascha@3008: int kmIndex = parameters.columnIndex("km"); sascha@3008: sascha@3008: ColumnCache cc = new ColumnCache(); sascha@3008: sascha@3008: double [] wq = new double[2]; sascha@3008: sascha@3008: int [] parameterIndices = sascha@3008: parameters.columnIndices(function.getParameterNames()); sascha@3008: sascha@3008: double [] parameterValues = new double[parameterIndices.length]; sascha@3008: sascha@3008: DateAverager dateAverager = new DateAverager(); sascha@3008: sascha@3008: for (int row = 0, P = parameters.size(); row < P; ++row) { sascha@3008: double km = parameters.get(row, kmIndex); sascha@3008: parameters.get(row, parameterIndices, parameterValues); sascha@3008: sascha@3008: // This is the paraterized function for a given km. sascha@3008: de.intevation.flys.artifacts.math.Function instance = sascha@3008: function.instantiate(parameterValues); sascha@3008: sascha@3008: KmFilter kmFilter = new KmFilter(km); sascha@3008: sascha@3008: for (DateRange analysisPeriod: analysisPeriods) { sascha@3008: DateRangeFilter drf = new DateRangeFilter( sascha@3008: analysisPeriod.getFrom(), sascha@3008: analysisPeriod.getTo()); sascha@3008: sascha@3008: // for all Q sectors. sascha@3008: for (int qSector = qSectorStart; qSector < qSectorEnd; ++qSector) { sascha@3008: sascha@3008: AndFilter filter = new AndFilter(); sascha@3008: filter.add(kmFilter); sascha@3008: filter.add(new SectorFilter(qSector)); sascha@3008: filter.add(drf); sascha@3008: sascha@3008: List metas = overview.filter(range, filter); sascha@3008: sascha@3008: if (metas.isEmpty()) { sascha@3008: // No fixings for km and analysis period sascha@3008: continue; sascha@3008: } sascha@3008: sascha@3008: double sumQ = 0.0; sascha@3008: double sumW = 0.0; sascha@3008: sascha@3008: List qwds = new ArrayList(metas.size()); sascha@3008: sascha@3008: dateAverager.clear(); sascha@3008: sascha@3008: for (Fixing.Column meta: metas) { sascha@3008: if (meta.findQSector(km) != qSector) { sascha@3008: // Ignore not matching sectors. sascha@3008: continue; sascha@3008: } sascha@3008: sascha@3008: Column column = cc.getColumn(meta); sascha@3008: if (column == null || !column.getWQ(km, wq)) { sascha@3008: continue; sascha@3008: } sascha@3008: sascha@3008: double fw = instance.value(wq[1]); sascha@3008: if (Double.isNaN(fw)) { sascha@3008: continue; sascha@3008: } sascha@3008: sascha@3008: double dw = (wq[0] - fw)*100.0; sascha@3008: sascha@3008: Date date = column.getDate(); sascha@3008: String description = column.getDescription(); sascha@3008: sascha@3008: QWD qwd = new QWD(wq[1], wq[0], dw, date, description); sascha@3008: sascha@3008: qwds.add(qwd); sascha@3008: sascha@3008: sumW += wq[0]; sascha@3008: sumQ += wq[1]; sascha@3008: sascha@3008: dateAverager.add(date); sascha@3008: } sascha@3008: sascha@3008: // Calulate average per Q sector. sascha@3008: int N = qwds.size(); sascha@3008: if (N > 0) { sascha@3008: double avgW = sumW / N; sascha@3008: double avgQ = sumQ / N; sascha@3008: sascha@3008: double avgFw = instance.value(avgQ); sascha@3008: if (!Double.isNaN(avgFw)) { sascha@3008: double avgDw = (avgW - avgFw)*100.0; sascha@3008: Date avgDate = dateAverager.getAverage(); sascha@3008: sascha@3008: String avgDescription = "avg.deltawt." + qSector; sascha@3008: sascha@3008: QWD avgQWD = new QWD( sascha@3008: avgQ, avgW, avgDw, avgDate, avgDescription); sascha@3008: sascha@3008: // TODO: Store average value. sascha@3008: } sascha@3008: else { sascha@3008: // TODO: Store sascha@3008: } sascha@3008: } sascha@3008: } // for all Q sectors sascha@3008: } sascha@3008: } sascha@3008: sascha@3008: return results; sascha@3008: } sascha@3008: 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@3008: sascha@3008: public Date getDate() { sascha@3008: return meta.getStartTime(); sascha@3008: } sascha@3008: sascha@3008: public String getDescription() { sascha@3008: return meta.getDescription(); sascha@3008: } sascha@3008: sascha@3008: public boolean getWQ(double km, double [] wq) { sascha@3008: data.getW(km, wq, 0); sascha@3008: if (Double.isNaN(wq[0])) return false; sascha@3008: wq[1] = data.getQ(km); sascha@3008: return !Double.isNaN(wq[1]); sascha@3008: } sascha@2744: } // class Column sascha@2744: sascha@3008: sascha@3008: /** sascha@3008: * Helper class to find the data belonging to meta info more quickly. sascha@3008: */ sascha@3008: public static class ColumnCache { sascha@3008: sascha@3008: protected Map columns; sascha@3008: sascha@3008: public ColumnCache() { sascha@3008: columns = new HashMap(); sascha@3008: } sascha@3008: sascha@3008: public Column getColumn(Fixing.Column meta) { sascha@3008: Integer key = meta.getId(); sascha@3008: Column column = columns.get(key); sascha@3008: if (column == null) { sascha@3008: FixingsColumn data = FixingsColumnFactory sascha@3008: .getInstance() sascha@3008: .getColumnData(meta); sascha@3008: if (data != null) { sascha@3008: column = new Column(meta, data); sascha@3008: columns.put(key, column); sascha@3008: } sascha@3008: } sascha@3008: return column; sascha@3008: } sascha@3008: } // class ColumnCache sascha@3008: sascha@2744: /** Fetch meta and data columns for analysis periods. */ sascha@2744: protected Column [][] getAnalysisColumns(FixingsOverview overview) { sascha@2744: sascha@2993: boolean debug = log.isDebugEnabled(); sascha@2993: if (debug) { sascha@2993: log.debug("number analysis periods: " + analysisPeriods.length); sascha@2993: } sascha@2993: 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@3003: DateRange period = analysisPeriods[i]; sascha@2744: sascha@2744: AndFilter filter = new AndFilter(); sascha@2744: sascha@2744: DateRangeFilter dateRangeFilter = sascha@3003: new DateRangeFilter(period.getFrom(), period.getTo()); 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@2993: if (debug) { sascha@2993: log.debug("number of filtered columns: " + metaCols.size()); sascha@2993: } sascha@2993: 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@2993: sascha@2993: if (debug) { sascha@2993: log.debug("failed loading: " + (metaCols.size()-cols.size())); sascha@2993: } 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@3010: String [] result = new String[parameters.length + 2]; sascha@2729: result[0] = "km"; sascha@3010: result[0] = "chi_sqr"; sascha@3010: System.arraycopy(parameters, 0, result, 2, parameters.length); sascha@2729: return result; sascha@2729: } sascha@2729: sascha@3010: protected Pair 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: for (int i = 0; i < ws.length; ++i) { sascha@2994: if (!Double.isNaN(ws[i]) && !Double.isNaN(qs[i])) { sascha@2994: cf.addObservedPoint(qs[i], ws[i]); sascha@2729: } sascha@2729: } sascha@2729: sascha@2729: try { sascha@3010: double [] parameters = sascha@3010: cf.fit(function, function.getInitialGuess()); sascha@3010: double chiSqr = lmo.getChiSquare(); sascha@3010: sascha@3010: return new Pair(parameters, chiSqr); sascha@2729: } sascha@2729: catch (MathException me) { sascha@2992: log.warn(me, 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 :