package iDSS.math;

import java.util.*;
import java.sql.Timestamp;



public class Fourier {



    private Vector x;
    private Vector y;
    private double period;
    private double halfP;

    private Vector a, b;
    private double c;

    private Double aDouble;

    private int MAX_NUMBER_OF_ITERATIONS = 0;

    public Fourier (Vector xValues, Vector yValues,
                    double timePeriod, int numberOfFourierIterations) {

        x = xValues;
        y = yValues;
        period = timePeriod;
        halfP = period/2.0;

        if (numberOfFourierIterations >= 0)
            MAX_NUMBER_OF_ITERATIONS = numberOfFourierIterations;


        aDouble = (Double) x.elementAt(0);
        double xStart = aDouble.doubleValue();

        aDouble = (Double) x.elementAt(x.size()-1);
        double xEnd = aDouble.doubleValue();
        double dataLength = xEnd-xStart;

        if (Math.abs(period - dataLength) > (period/1000.0)) {
            System.out.println("period, dataLength = " + period +", " + dataLength);
            System.out.println("Length of observed data not sufficient");
            System.exit(0);
        }


        a = new Vector();
        b = new Vector();
        c = 0;

        calculateC();
        calculateA();
        calculateB();

    } //public Fourier (Vector xValues, Vector yValues, double timePeriod)







    public Fourier (Vector columnsOfData, int numberOfFourierIterations) {


        Vector xValues = (Vector) columnsOfData.elementAt(0);
        Vector yValues = (Vector) columnsOfData.elementAt(1);

        if (numberOfFourierIterations >= 0)
            MAX_NUMBER_OF_ITERATIONS = numberOfFourierIterations;


        x = new Vector();
        y = yValues;

        Timestamp firstDay = (Timestamp) xValues.elementAt(0);
        double firstDayValue = (double) firstDay.getTime();

        for(int i=0;i<xValues.size();i++) {

            Timestamp time = (Timestamp) xValues.elementAt(i);
            double aVal = (double) time.getTime();


            //aVal = aVal-firstDayValue;
            x.add(new Double(aVal));

        } //for(int i=0;i<xValues.size();i++)




        aDouble = (Double) x.elementAt(0);
        double xStart = aDouble.doubleValue();

        aDouble = (Double) x.elementAt(x.size()-1);
        double xEnd = aDouble.doubleValue();
        double dataLength = xEnd-xStart;
        period = dataLength;
        halfP = period/2.0;



        if (Math.abs(period - dataLength) > (period/1000.0)) {
            System.out.println("period, dataLength = " + period +", " + dataLength);
            System.out.println("Length of observed data not sufficient");
            System.exit(0);
        }


        a = new Vector();
        b = new Vector();
        c = 0;

        calculateC();
        calculateA();
        calculateB();

    } //public Fourier (Vector xValues, Vector yValues, double timePeriod)




    public Vector getFourierA() {

        return a;

    } //public Vector getFourierA()





    public Vector getFourierB() {

        return b;

    } //public Vector getFourierB()




    public double getFourierC() {

        return c;

    } //public double getFourierC()




    public int getFourierIterations() {

        return MAX_NUMBER_OF_ITERATIONS;

    } //public Vector getFourierIterations()





    public double getFourierPeriod() {

        return period;

    } //public double getFourierperiod()






    public void calculateC() {


        System.out.println("Calculating Fourier coefficient c");

        aDouble = (Double)x.elementAt(0);
        double x0 = aDouble.doubleValue();

        aDouble = (Double)y.elementAt(0);
        double y0 = aDouble.doubleValue();

        double x1;
        double y1;

        for (int i=1;i<x.size();i++) {

            aDouble = (Double)x.elementAt(i);
            x1 = aDouble.doubleValue();

            aDouble = (Double)y.elementAt(i);
            y1 = aDouble.doubleValue();
            c = c+(x1-x0)*(y0+y1)/2/halfP/2;

            x0 = x1;
            y0 = y1;

        } //for (int i=1;i<x.size();i++)

        System.out.println("c = " + c);


    } //public void calculateC()





    public void calculateA() {

        System.out.println("Calculating Fourier coefficient a");

        int n = 0;
        boolean keepGoing = true;

        while(keepGoing) {

            n++;

            aDouble = (Double)x.elementAt(0);
            double x0 = aDouble.doubleValue();

            aDouble = (Double)y.elementAt(0);
            double y0 = aDouble.doubleValue();

            y0 = y0 * Math.cos(n*Math.PI*x0/halfP);

            double x1;
            double y1;

            double a1 = 0;

            for (int i=1;i<x.size();i++) {


                aDouble = (Double)x.elementAt(i);
                x1 = aDouble.doubleValue();

                aDouble = (Double)y.elementAt(i);
                y1 = aDouble.doubleValue();

                y1 = y1 * Math.cos(n*Math.PI*x1/halfP);

                a1 = a1+Math.abs(x1-x0)*(y0+y1)/2;

                x0 = x1;
                y0 = y1;

            } //for (int i=1;i<x.size();i++)

            a1 = a1/halfP;
            a.add(new Double(a1));


            //if((Math.abs(a1) < 0.001) || (n > 1000))
            if(n > MAX_NUMBER_OF_ITERATIONS)
                keepGoing = false;

        } //while(keepGoing)



    } //public calculateA()






    public void calculateB() {

        System.out.println("Calculating Fourier coefficient b");

        int n = 0;
        boolean keepGoing = true;

        while(keepGoing) {

            n++;

            aDouble = (Double)x.elementAt(0);
            double x0 = aDouble.doubleValue();

            aDouble = (Double)y.elementAt(0);
            double y0 = aDouble.doubleValue();

            y0 = y0 * Math.sin(n*Math.PI*x0/halfP);

            double x1;
            double y1;

            double b1 = 0;

            for (int i=1;i<x.size();i++) {


                aDouble = (Double)x.elementAt(i);
                x1 = aDouble.doubleValue();

                aDouble = (Double)y.elementAt(i);
                y1 = aDouble.doubleValue();

                y1 = y1 * Math.sin(n*Math.PI*x1/halfP);

                b1 = b1+Math.abs(x1-x0)*(y0+y1)/2.0;

                x0 = x1;
                y0 = y1;

            } //for (int i=1;i<x.size();i++)

            b1 = b1/halfP;
            b.add(new Double(b1));


            //if((Math.abs(b1) < 0.001) || (n > 1000))
            if(n > MAX_NUMBER_OF_ITERATIONS)
                keepGoing = false;

        } //while(keepGoing) 


    } //public calculateB()











    public double getComputedY(double x) {

        double xVal = x;

        double componentA = 0;

        for(int j=0;j<a.size();j++) {

            int n = j+1;

            aDouble = (Double)a.elementAt(j);
            double aVal = aDouble.doubleValue();

            componentA = componentA +
                         aVal*Math.cos(n*Math.PI*xVal/halfP);

        } //for(int j=0;j<a.size();j++) for componentA

        double componentB = 0;

        for(int j=0;j<b.size();j++) {

            int n = j+1;

            aDouble = (Double)b.elementAt(j);
            double bVal = aDouble.doubleValue();

            componentB = componentB +
                         bVal*Math.sin(n*Math.PI*xVal/halfP);

        } //for(int j=0;j<a.size();j++) for componentB

        double yComputed = c+componentA+componentB;

        return yComputed;

    } //public double getComputedY(double x)











    public static double computeY(double x,Vector aa,Vector bb,double cc,double timePeriod) {

        double xVal = x;
        double halfPeriod = timePeriod/2.0;
        double componentA = 0;
        Double bDouble;

        for(int j=0;j<aa.size();j++) {

            int n = j+1;

            bDouble = (Double)aa.elementAt(j);
            double aVal = bDouble.doubleValue();

            componentA = componentA +
                         aVal*Math.cos(n*Math.PI*xVal/halfPeriod);

        } //for(int j=0;j<aa.size();j++) for componentA

        double componentB = 0;

        for(int j=0;j<bb.size();j++) {

            int n = j+1;

            bDouble = (Double)bb.elementAt(j);
            double bVal = bDouble.doubleValue();

            componentB = componentB +
                         bVal*Math.sin(n*Math.PI*xVal/halfPeriod);

        } //for(int j=0;j<aa.size();j++) for componentB

        double yComputed = cc+componentA+componentB;

        return yComputed;

    } //public static double computeY(double x,Vector aa,Vector bb,double cc,double timePeriod)


















/****


    public DataSet getComputedDataSet() {

        long xVal;
        double yComputed;

        Vector rows = new Vector();
        Vector columns = new Vector();

        for(int i=0;i<x.size();i++) {

            aDouble = (Double)x.elementAt(i);
            xVal = aDouble.longValue();
            Timestamp time = new Timestamp(xVal);

            yComputed = getComputedY(aDouble.doubleValue());

            Vector currentRow = new Vector();
            currentRow.add(time);
            currentRow.add(new Double(yComputed));

            rows.add(currentRow);

        } //for(int i=0;i<x.size();i++)


        columns.add("Date");
        columns.add("Computed Value");



    } //public DataSet getComputedDataSet();


****/












    public void compareValues() {

        double xVal;
        double yVal;
        double yComputed;

        for(int i=0;i<x.size();i++) {

            aDouble = (Double)x.elementAt(i);
            xVal = aDouble.doubleValue();

            aDouble = (Double)y.elementAt(i);
            yVal = aDouble.doubleValue();

            yComputed = getComputedY(xVal);


            System.out.println(xVal + ", " + yVal + ", " + yComputed);

        } //for(int i=0;i<x.size();i++)



    } //public void compareValues()






    public static void main(String args[]) {



        Vector xVector = new Vector();
        Vector yVector = new Vector();

        double period = 48;


        xVector.add(Timestamp.valueOf("1999-01-01 00:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 01:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 02:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 03:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 04:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 05:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 06:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 07:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 08:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 09:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 10:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 11:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 12:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 13:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 14:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 15:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 16:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 17:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 18:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 19:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 20:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 21:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 22:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 23:00:00"));
        xVector.add(Timestamp.valueOf("1999-01-01 24:00:00"));


        yVector.add(new Double( 2.595 ));
        yVector.add(new Double( 1.899 ));
        yVector.add(new Double( 1.341 ));
        yVector.add(new Double( 0.850 ));
        yVector.add(new Double( 0.587 ));
        yVector.add(new Double( 0.793 ));
        yVector.add(new Double( 1.360 ));
        yVector.add(new Double( 1.935 ));
        yVector.add(new Double( 2.562 ));
        yVector.add(new Double( 3.218 ));
        yVector.add(new Double( 3.795 ));
        yVector.add(new Double( 3.758 ));
        yVector.add(new Double( 3.019 ));
        yVector.add(new Double( 2.230 ));
        yVector.add(new Double( 1.619 ));
        yVector.add(new Double( 1.048 ));
        yVector.add(new Double( 0.628 ));
        yVector.add(new Double( 0.574 ));
        yVector.add(new Double( 0.972 ));
        yVector.add(new Double( 1.561 ));
        yVector.add(new Double( 2.165 ));
        yVector.add(new Double( 2.801 ));
        yVector.add(new Double( 3.454 ));
        yVector.add(new Double( 3.848 ));
        yVector.add(new Double( 3.379 ));



        Vector colsOfData = new Vector();
        colsOfData.add(xVector);
        colsOfData.add(yVector);

        Fourier f = new Fourier(colsOfData,100);

        f.compareValues();

        Timestamp anyTime = Timestamp.valueOf("1999-01-01 00:00:00");

        double xx = (double)anyTime.getTime();
        double yy = f.computeY(xx,f.getFourierA(),f.getFourierB(),f.getFourierC(),f.getFourierPeriod());

        System.out.println("xx = " + anyTime.toString());
        System.out.println("yy = " + yy);


        System.out.println("Extrapolation starts....\n");

        for (int i=0;i<xVector.size();i++) {

            anyTime = (Timestamp)xVector.elementAt(i);
            Calendar cal = Calendar.getInstance();
            cal.setTime(anyTime);
            cal.add(Calendar.DATE,1);
            anyTime.setTime(cal.getTime().getTime());
            xx = (double)cal.getTime().getTime();
            yy = f.computeY(xx,f.getFourierA(),f.getFourierB(),f.getFourierC(),f.getFourierPeriod());

            System.out.println(xx + "," + anyTime.toString() + "," + yy);

        } //for (int i=0;i<xVector.size();i++)



    } //public static void main(String args[])




} //public class Fourier




