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: 
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: 
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@6787:                 powFunction.getParameterNames());
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@6787:                 function.getParameterNames());
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<SQ> sqs = measurements.getSQs(i);
sascha@3310: 
ingo@3426:             SQFractionResult fractionResult;
sascha@3310: 
sascha@3552:             List<SQFractionResult.Iteration> 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<SQFractionResult.Iteration> doFitting(
teichmann@6787:         final Function         function,
teichmann@6787:         List<SQ>               sqs,
teichmann@6787:         SQ.View                sqView,
teichmann@6787:         final ParameterCreator pc
sascha@3552:     ) {
sascha@3552:         final List<SQFractionResult.Iteration> iterations =
sascha@3552:             new ArrayList<SQFractionResult.Iteration>();
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,
sascha@3552:                         chiSqr);
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@6787:         public ParameterCreator(String [] origNames, String [] proxyNames) {
teichmann@6787:             this.origNames  = origNames;
teichmann@6787:             this.proxyNames = proxyNames;
teichmann@6787:         }
teichmann@6787: 
teichmann@6787:         protected double [] transformCoeffs(double [] coeffs) {
teichmann@6787:             return coeffs;
teichmann@6787:         }
teichmann@6787: 
teichmann@6787:         public Parameters createParameters(
teichmann@6787:             double [] coeffs,
teichmann@6787:             double    standardDeviation,
teichmann@6787:             double    chiSqr
teichmann@6787:         ) {
teichmann@6787:             String [] columns = new String[origNames.length + 2];
teichmann@6787:             columns[0] = "chi_sqr";
teichmann@6787:             columns[1] = "std_dev";
teichmann@6787:             System.arraycopy(origNames, 0, columns, 2, origNames.length);
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@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@6787:             String [] proxyNames
teichmann@6787:         ) {
teichmann@6787:             super(origNames, proxyNames);
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 :