teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde teichmann@5863: * Software engineering by Intevation GmbH teichmann@5863: * teichmann@5994: * This file is Free Software under the GNU AGPL (>=v3) teichmann@5863: * and comes with ABSOLUTELY NO WARRANTY! Check out the teichmann@5994: * documentation coming with Dive4Elements River for details. teichmann@5863: */ teichmann@5863: teichmann@5831: package org.dive4elements.river.artifacts.model.sq; ingo@3072: teichmann@6787: import org.dive4elements.artifacts.common.utils.StringUtils; teichmann@5831: import org.dive4elements.river.artifacts.access.SQRelationAccess; sascha@3304: teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.Function; teichmann@5831: import org.dive4elements.river.artifacts.math.fitting.FunctionFactory; ingo@3072: teichmann@5831: import org.dive4elements.river.artifacts.model.Calculation; teichmann@5831: import org.dive4elements.river.artifacts.model.CalculationResult; teichmann@5831: import org.dive4elements.river.artifacts.model.DateRange; teichmann@5831: import org.dive4elements.river.artifacts.model.Parameters; teichmann@5831: teichmann@5831: import org.dive4elements.river.backend.SedDBSessionHolder; sascha@3289: sascha@3552: import java.util.ArrayList; teichmann@6787: import java.util.Arrays; sascha@3304: import java.util.List; sascha@3304: sascha@3222: import org.apache.log4j.Logger; ingo@3072: ingo@3072: public class SQRelationCalculation extends Calculation { ingo@3072: sascha@3222: private static final Logger log = ingo@3072: Logger.getLogger(SQRelationCalculation.class); ingo@3072: teichmann@6787: public static final boolean NON_LINEAR_FITTING = teichmann@6787: Boolean.getBoolean("minfo.sq.calcution.non.linear.fitting"); teichmann@6787: teichmann@6787: public static final String SQ_POW_FUNCTION_NAME = "sq-pow"; teichmann@6787: public static final String SQ_LIN_FUNCTION_NAME = "linear"; sascha@3304: teichmann@6796: public static final String [] EXTRA_PARAMETERS = { teichmann@6796: "chi_sqr", teichmann@6796: "std_dev", teichmann@6796: "max_q", teichmann@6797: "c_ferguson", teichmann@6796: "c_duan", teichmann@6796: "r2" teichmann@6796: }; teichmann@6796: ingo@3426: protected String river; ingo@3426: protected double location; ingo@3426: protected DateRange period; ingo@3426: protected double outliers; rrenkert@5396: private String method; sascha@3305: sascha@3222: public SQRelationCalculation() { ingo@3079: } ingo@3079: sascha@3222: public SQRelationCalculation(SQRelationAccess access) { ingo@3079: ingo@3426: String river = access.getRiver(); ingo@3426: Double location = access.getLocation(); sascha@3552: DateRange period = access.getPeriod(); ingo@3426: Double outliers = access.getOutliers(); rrenkert@5396: String method = access.getOutlierMethod(); ingo@3101: sascha@3222: if (river == null) { sascha@3222: // TODO: i18n sascha@3222: addProblem("sq.missing.river"); sascha@3222: } ingo@3101: sascha@3222: if (location == null) { sascha@3222: // TODO: i18n sascha@3222: addProblem("sq.missing.location"); sascha@3222: } ingo@3101: ingo@3426: if (period == null) { sascha@3222: // TODO: i18n sascha@3222: addProblem("sq.missing.periods"); sascha@3222: } sascha@3222: sascha@3222: if (outliers == null) { sascha@3222: // TODO: i18n sascha@3222: addProblem("sq.missing.outliers"); sascha@3222: } sascha@3222: rrenkert@5396: if (method == null) { rrenkert@5396: //TODO: i18n rrenkert@5396: addProblem("sq.missing.method"); rrenkert@5396: } rrenkert@5396: sascha@3222: if (!hasProblems()) { sascha@3222: this.river = river; sascha@3222: this.location = location; ingo@3426: this.period = period; sascha@3222: this.outliers = outliers; rrenkert@5396: this.method = method; sascha@3222: } ingo@3079: } ingo@3079: ingo@3105: sascha@3222: public CalculationResult calculate() { sascha@3222: log.debug("SQRelationCalculation.calculate"); ingo@3079: sascha@3222: if (hasProblems()) { sascha@3222: return new CalculationResult(this); ingo@3118: } ingo@3079: sascha@3289: SedDBSessionHolder.acquire(); sascha@3289: try { sascha@3289: return internalCalculate(); sascha@3289: } sascha@3289: finally { sascha@3289: SedDBSessionHolder.release(); sascha@3289: } sascha@3289: } sascha@3289: teichmann@6787: public interface TransformCoeffs { teichmann@6787: double [] transform(double [] coeffs); teichmann@6787: } teichmann@6787: teichmann@6787: public static final TransformCoeffs IDENTITY_TRANS = teichmann@6787: new TransformCoeffs() { teichmann@6787: @Override teichmann@6787: public double [] transform(double [] coeffs) { teichmann@6787: return coeffs; teichmann@6787: } teichmann@6787: }; teichmann@6787: teichmann@6787: public static final TransformCoeffs LINEAR_TRANS = teichmann@6787: new TransformCoeffs() { teichmann@6787: @Override teichmann@6787: public double [] transform(double [] coeffs) { teichmann@6787: log.debug("before transform: " + Arrays.toString(coeffs)); teichmann@6787: if (coeffs.length == 2) { teichmann@6787: coeffs = new double [] { Math.exp(coeffs[1]), coeffs[0] }; teichmann@6787: } teichmann@6787: log.debug("after transform: " + Arrays.toString(coeffs)); teichmann@6787: return coeffs; teichmann@6787: } teichmann@6787: }; teichmann@6787: sascha@3289: protected CalculationResult internalCalculate() { sascha@3289: teichmann@6787: Function powFunction = FunctionFactory sascha@3304: .getInstance() teichmann@6787: .getFunction(SQ_POW_FUNCTION_NAME); sascha@3304: teichmann@6787: if (powFunction == null) { teichmann@6787: log.error("No '" + SQ_POW_FUNCTION_NAME + "' function found."); sascha@3304: // TODO: i18n sascha@3304: addProblem("sq.missing.sq.function"); teichmann@6787: return new CalculationResult(new SQResult[0], this); sascha@3304: } sascha@3304: teichmann@6787: Function function; teichmann@6787: SQ.View sqView; teichmann@6787: SQ.Factory sqFactory; teichmann@6787: ParameterCreator pc; teichmann@6787: teichmann@6787: if (NON_LINEAR_FITTING) { teichmann@6787: log.debug("Use non linear fitting."); teichmann@6787: sqView = SQ.SQ_VIEW; teichmann@6787: sqFactory = SQ.SQ_FACTORY; teichmann@6787: function = powFunction; teichmann@6787: pc = new ParameterCreator( teichmann@6787: powFunction.getParameterNames(), teichmann@6796: powFunction.getParameterNames(), teichmann@6796: powFunction, teichmann@6796: sqView); teichmann@6787: } teichmann@6787: else { teichmann@6787: log.debug("Use linear fitting."); teichmann@6787: sqView = LogSQ.LOG_SQ_VIEW; teichmann@6787: sqFactory = LogSQ.LOG_SQ_FACTORY; teichmann@6787: function = FunctionFactory teichmann@6787: .getInstance() teichmann@6787: .getFunction(SQ_LIN_FUNCTION_NAME); teichmann@6787: if (function == null) { teichmann@6787: log.error("No '" + SQ_LIN_FUNCTION_NAME + "' function found."); teichmann@6787: // TODO: i18n teichmann@6787: addProblem("sq.missing.sq.function"); teichmann@6787: return new CalculationResult(new SQResult[0], this); teichmann@6787: } teichmann@6787: pc = new LinearParameterCreator( teichmann@6787: powFunction.getParameterNames(), teichmann@6796: function.getParameterNames(), teichmann@6796: function, teichmann@6796: sqView); teichmann@6787: } teichmann@6780: sascha@3552: Measurements measurements = teichmann@6780: MeasurementFactory.getMeasurements( teichmann@6780: river, location, period, sqFactory); sascha@3310: ingo@3426: SQFractionResult [] fractionResults = ingo@3426: new SQFractionResult[SQResult.NUMBER_FRACTIONS]; sascha@3304: teichmann@6787: ingo@3426: for (int i = 0; i < fractionResults.length; ++i) { ingo@3426: List sqs = measurements.getSQs(i); sascha@3310: ingo@3426: SQFractionResult fractionResult; sascha@3310: sascha@3552: List iterations = teichmann@6787: doFitting(function, sqs, sqView, pc); sascha@3552: sascha@3552: if (iterations == null) { ingo@3426: // TODO: i18n ingo@3426: addProblem("sq.fitting.failed." + i); ingo@3426: fractionResult = new SQFractionResult(); ingo@3426: } ingo@3426: else { sascha@3552: fractionResult = new SQFractionResult( sascha@3552: sqs.toArray(new SQ[sqs.size()]), sascha@3552: iterations); sascha@3552: } sascha@3310: ingo@3426: fractionResults[i] = fractionResult; sascha@3297: } sascha@3289: sascha@3310: return new CalculationResult( ingo@3426: new SQResult[] { new SQResult(location, fractionResults) }, sascha@3310: this); ingo@3072: } sascha@3304: sascha@3552: protected List doFitting( teichmann@6787: final Function function, teichmann@6787: List sqs, teichmann@6787: SQ.View sqView, teichmann@6787: final ParameterCreator pc sascha@3552: ) { sascha@3552: final List iterations = sascha@3552: new ArrayList(); sascha@3304: teichmann@6780: boolean success = new Fitting(function, outliers, sqView).fit( sascha@3552: sqs, rrenkert@5396: method, sascha@3552: new Fitting.Callback() { sascha@3552: @Override sascha@3552: public void afterIteration( sascha@3552: double [] coeffs, sascha@3552: SQ [] measurements, sascha@3552: SQ [] outliers, sascha@3552: double standardDeviation, sascha@3552: double chiSqr sascha@3552: ) { teichmann@6787: Parameters parameters = pc.createParameters( sascha@3552: coeffs, sascha@3552: standardDeviation, teichmann@6796: chiSqr, teichmann@6796: measurements); sascha@3552: iterations.add(new SQFractionResult.Iteration( sascha@3552: parameters, sascha@3552: measurements, sascha@3552: outliers)); sascha@3552: } sascha@3552: }); sascha@3552: sascha@3552: return success ? iterations : null; sascha@3552: } sascha@3552: teichmann@6787: public static class ParameterCreator { teichmann@6787: teichmann@6787: protected String [] origNames; teichmann@6787: protected String [] proxyNames; teichmann@6787: teichmann@6796: protected Function function; teichmann@6796: protected SQ.View view; teichmann@6796: teichmann@6796: public ParameterCreator( teichmann@6796: String [] origNames, teichmann@6796: String [] proxyNames, teichmann@6796: Function function, teichmann@6796: SQ.View view teichmann@6796: ) { teichmann@6787: this.origNames = origNames; teichmann@6787: this.proxyNames = proxyNames; teichmann@6796: this.view = view; teichmann@6787: } teichmann@6787: teichmann@6787: protected double [] transformCoeffs(double [] coeffs) { teichmann@6787: return coeffs; teichmann@6787: } teichmann@6787: teichmann@6796: private static double maxQ(SQ [] sqs) { teichmann@6796: double max = -Double.MAX_VALUE; teichmann@6796: for (SQ sq: sqs) { teichmann@6796: double q = sq.getQ(); // Don't use view here! teichmann@6796: if (q > max) { teichmann@6796: max = q; teichmann@6796: } teichmann@6796: } teichmann@6796: return Math.max(0d, max); teichmann@6796: } teichmann@6796: teichmann@6797: private double cFerguson( teichmann@6796: org.dive4elements.river.artifacts.math.Function instance, teichmann@6796: SQ [] sqs teichmann@6796: ) { teichmann@6796: double sqrSum = 0d; teichmann@6796: teichmann@6796: for (SQ sq: sqs) { teichmann@6796: double s = view.getS(sq); teichmann@6796: double q = view.getQ(sq); teichmann@6796: double diffS = s - instance.value(q); teichmann@6796: sqrSum += diffS*diffS; teichmann@6796: } teichmann@6796: teichmann@6796: return Math.exp(0.5d * sqrSum/(sqs.length-2)); teichmann@6796: } teichmann@6796: teichmann@6796: private double cDuan( teichmann@6796: org.dive4elements.river.artifacts.math.Function instance, teichmann@6796: SQ [] sqs teichmann@6796: ) { teichmann@6796: double sum = 0d; teichmann@6796: teichmann@6796: for (SQ sq: sqs) { teichmann@6796: double s = view.getS(sq); teichmann@6796: double q = view.getQ(sq); teichmann@6796: double diffS = s - instance.value(q); teichmann@6796: sum += Math.exp(diffS); teichmann@6796: } teichmann@6796: return sum / sqs.length; teichmann@6796: } teichmann@6796: teichmann@6796: private double r2( teichmann@6796: org.dive4elements.river.artifacts.math.Function instance, teichmann@6796: SQ [] sqs teichmann@6796: ) { teichmann@6796: double xm = 0; teichmann@6796: double ym = 0; teichmann@6796: for (SQ sq: sqs) { teichmann@6796: double s = view.getS(sq); teichmann@6796: double q = view.getQ(sq); teichmann@6796: double fs = instance.value(q); teichmann@6796: xm += s; teichmann@6796: ym += fs; teichmann@6796: } teichmann@6796: xm /= sqs.length; teichmann@6796: ym /= sqs.length; teichmann@6796: teichmann@6796: double mixXY = 0d; teichmann@6796: double sumX = 0d; teichmann@6796: double sumY = 0d; teichmann@6796: teichmann@6796: for (SQ sq: sqs) { teichmann@6796: double s = view.getS(sq); teichmann@6796: double q = view.getQ(sq); teichmann@6796: double fs = instance.value(q); teichmann@6796: teichmann@6796: double xDiff = xm - s; teichmann@6796: double yDiff = ym - fs; teichmann@6796: teichmann@6796: mixXY += xDiff*yDiff; teichmann@6796: teichmann@6796: sumX += xDiff*xDiff; teichmann@6796: sumY *= yDiff*yDiff; teichmann@6796: } teichmann@6796: teichmann@6796: double r = mixXY/Math.sqrt(sumX*sumY); teichmann@6796: return r*r; teichmann@6796: } teichmann@6796: teichmann@6796: teichmann@6787: public Parameters createParameters( teichmann@6787: double [] coeffs, teichmann@6787: double standardDeviation, teichmann@6796: double chiSqr, teichmann@6796: SQ [] measurements teichmann@6787: ) { teichmann@6796: String [] columns = StringUtils.join(EXTRA_PARAMETERS, origNames); teichmann@6796: teichmann@6787: Parameters parameters = new Parameters(columns); teichmann@6787: int row = parameters.newRow(); teichmann@6787: parameters.set(row, origNames, transformCoeffs(coeffs)); teichmann@6787: parameters.set(row, "chi_sqr", chiSqr); teichmann@6787: parameters.set(row, "std_dev", standardDeviation); teichmann@6796: parameters.set(row, "max_q", maxQ(measurements)); teichmann@6796: teichmann@6796: // We need to instantiate the function to calculate teichmann@6796: // the remaining values. teichmann@6796: org.dive4elements.river.artifacts.math.Function f = teichmann@6796: function.instantiate(coeffs); teichmann@6796: teichmann@6797: parameters.set(row, "c_ferguson", cFerguson(f, measurements)); teichmann@6796: parameters.set(row, "c_duan", cDuan(f, measurements)); teichmann@6796: parameters.set(row, "r2", r2(f, measurements)); teichmann@6796: teichmann@6787: return parameters; teichmann@6787: } teichmann@6787: } teichmann@6787: teichmann@6787: /** We need to transform the coeffs back to the original function. */ teichmann@6787: public static class LinearParameterCreator extends ParameterCreator { teichmann@6787: teichmann@6787: public LinearParameterCreator( teichmann@6787: String [] origNames, teichmann@6796: String [] proxyNames, teichmann@6796: Function function, teichmann@6796: SQ.View view teichmann@6787: ) { teichmann@6796: super(origNames, proxyNames, function, view); teichmann@6787: } teichmann@6787: teichmann@6787: @Override teichmann@6787: protected double [] transformCoeffs(double [] coeffs) { teichmann@6787: teichmann@6787: int bP = StringUtils.indexOf("m", proxyNames); teichmann@6787: int mP = StringUtils.indexOf("b", proxyNames); teichmann@6787: teichmann@6787: int aO = StringUtils.indexOf("a", origNames); teichmann@6787: int bO = StringUtils.indexOf("b", origNames); teichmann@6787: teichmann@6787: if (bP == -1 || mP == -1 || aO == -1 || bO == -1) { teichmann@6787: log.error("index not found: " teichmann@6787: + bP + " " + mP + " " teichmann@6787: + aO + " " + bO); teichmann@6787: return coeffs; teichmann@6787: } teichmann@6787: teichmann@6787: double [] ncoeffs = (double [])coeffs.clone(); teichmann@6787: ncoeffs[aO] = Math.exp(coeffs[mP]); teichmann@6787: ncoeffs[bO] = coeffs[bP]; teichmann@6787: teichmann@6787: if (log.isDebugEnabled()) { teichmann@6787: log.debug("before transform: " + Arrays.toString(coeffs)); teichmann@6787: log.debug("after transform: " + Arrays.toString(ncoeffs)); teichmann@6787: } teichmann@6787: teichmann@6787: return ncoeffs; teichmann@6787: } sascha@3304: } ingo@3072: } ingo@3072: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :