package de.intevation.flys.artifacts.math;

import java.util.ArrayList;
import java.util.List;

import java.io.Serializable;

import org.apache.commons.math.analysis.interpolation.SplineInterpolator;

import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;

import org.apache.commons.math.ArgumentOutsideDomainException;

import org.apache.log4j.Logger;

public class BackJumpCorrector
implements   Serializable
{
    private static Logger log = Logger.getLogger(BackJumpCorrector.class);

    protected ArrayList<Double> backjumps;

    protected double [] corrected;

    public BackJumpCorrector() {
        backjumps = new ArrayList<Double>();
    }

    public boolean hasBackJumps() {
        return !backjumps.isEmpty();
    }

    public List<Double> getBackJumps() {
        return backjumps;
    }

    public double [] getCorrected() {
        return corrected;
    }

    private static final boolean isIncreasing(double [] ws) {
        int inc = 0;
        int dec = 0;
        int sweet = (ws.length-1)/2;
        for (int i = 1; i < ws.length; ++i) {
            if (ws[i] > ws[i-1]) {
                if (++inc > sweet) {
                    return true;
                }
            }
            else if (++dec > sweet) {
                return false;
            }
        }
        return inc > sweet;
    }

    public boolean doCorrection(double [] km, double [] ws) {
        boolean isIncreasing = isIncreasing(ws);

        if (isIncreasing) {
            // mirror along x axis to simulate decreasing values
            ws = (double [])ws.clone();
            for (int i = 0; i < ws.length; ++i) {
                ws[i] = -ws[i];
            }
        }

        boolean hasBackJumps = doCorrectionClean(km, ws);

        if (hasBackJumps && isIncreasing) {
            // mirror back
            for (int i = 0; i < corrected.length; ++i) {
                corrected[i] = -corrected[i];
            }
        }

        return hasBackJumps;
    }

    protected boolean doCorrectionClean(double [] km, double [] ws) {
        int N = km.length;

        if (N != ws.length) {
            throw new IllegalArgumentException("km.length != ws.length");
        }

        if (N < 2) {
            return false;
        }

        SplineInterpolator interpolator = null;

        for (int i = 1; i < N; ++i) {
            if (ws[i] <= ws[i-1]) {
                // no back jump
                continue;
            }
            backjumps.add(km[i]);

            if (corrected == null) {
                // lazy cloning
                ws = corrected = (double [])ws.clone();
            }

            int    back = i-2;
            double distance = Math.abs(km[i] - km[i-1]);
            double rest = 0.0;
            double ikm  = 0.0;

            while (back > -1) {
                if (ws[back] < ws[i]) { // under water
                    // continue scanning backwards
                    distance += Math.abs(km[back+1] - km[back]);
                    --back;
                    continue;
                }
                if (ws[back] > ws[i]) { // over water
                    // need to calculate intersection
                    log.debug("intersection case");
                    double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]);
                    double b = km[back]-ws[back]*m;
                    ikm = m*ws[i] + b;
                    distance += Math.abs(ikm - km[back+1]);
                    rest = Math.abs(km[back] - ikm);
                }
                else {
                    // equals height at ws[back]
                    log.debug("exact height match");
                    distance += Math.abs(km[back+1] - km[back]);
                    ikm = km[back];
                }
                break;
            }

            if (back < 0) {
                log.debug("run over left border");
                // fill with same as ws[i]
                for (int j = 0; j < i; ++j) {
                    ws[j] = ws[i];
                }
                continue;
            }

            double quarterDistance = 0.25*distance;

            // Now further seek back for the max height

            double restDistance = Math.max(0.0, quarterDistance - rest);
            --back;

            double mkmw = ws[i];
            double mkm  = km[0];

            while (back > -1) {
                double d = Math.abs(km[back+1] - km[back]);
                restDistance -= d;
                if (restDistance > 0.0) {
                    --back;
                    continue;
                }
                if (restDistance < 0.0) {
                    // need to calculate intersection
                    if (km[back] == km[back+1]) { // should not happen
                        mkm = km[back];
                        mkmw = 0.5*(ws[back] + ws[back+1]);
                    }
                    else {
                        double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]);
                        double b = ws[back] - km[back]*m;
                        mkm = km[back] + restDistance;
                        mkmw = m*mkm + b;
                    }
                }
                else {
                    // exact match
                    mkm  = km[back];
                    mkmw = ws[back];
                }
                break;
            }

            double factor = back >= 0 && Math.abs(restDistance) < 1e-4
                ? 1.0
                : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance));

            double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i];

            double end = ikm + quarterDistance*factor;

            double [] x = { mkm,  ikm,  end   };
            double [] y = { mkmw, ikmw, ws[i] };

            if (interpolator == null) {
                interpolator = new SplineInterpolator();
            }

            PolynomialSplineFunction spline = interpolator.interpolate(x, y);

            try {
                if (log.isDebugEnabled()) {
                    log.debug("spline points:");
                    for (int j = 0; j < x.length; ++j) {
                        log.debug(x[j] + " " + y[j] + " " + spline.value(x[j]));
                    }
                }

                for (back = Math.max(back, 0);
                     back < i && km[back] < end; 
                     ++back
                ) {
                    // to 3/4 point fill with spline values
                    ws[back] = spline.value(km[back]);
                }
                while (back < i) {
                    // fill the rest with ws[i]
                    ws[back++] = ws[i];
                }
            }
            catch (ArgumentOutsideDomainException aode) {
                log.error("spline interpolation failed", aode);
            }
        } // for all km

        return !backjumps.isEmpty();
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :