/*
 * LSpec.java
 *
 * Created on April 24, 2007, 11:06 PM
 */

package org.das2.qds.util;

import java.util.LinkedHashMap;
import java.util.Map;
import org.das2.datum.Datum;
import org.das2.datum.Units;
import org.das2.qds.ArrayDataSet;
import org.das2.qds.DDataSet;
import org.das2.qds.DataSetUtil;
import org.das2.qds.QDataSet;
import org.das2.qds.SemanticOps;
import org.das2.qds.ops.Ops;

/**
 * Form a rank 2 dataset with L and Time for tags by identifying monotonic sweeps
 * in two rank 1 datasets, and interpolating along the sweep.
 * 
 * @author jbf
 */
public class LSpec {

    /**
     * identifies the sweep of each record
     */
    public static final String USER_PROP_SWEEPS = "sweeps";
    
    /** Creates a new instance of LSpec */
    private LSpec() {
    }

    /**
     * just the inward sweeps
     */
    public static int DIR_INWARD= -1;
    
    /**
     * just the outward sweeps
     */
    public static int DIR_OUTWARD= 1;
    
    /**
     * include both sweeps
     */
    public static int DIR_BOTH= 0;
        
    /**
     * identify monotonically increasing or decreasing segments of the dataset.
     * @param lds the dataset which sweeps back and forth, such as LShell or MagLat (or a sine wave for testing).
     * @param dir 0=both 1=outward 2= inward
     * @return rank 2 data set of sweep indeces, dim0 is sweep number, dim1 is two-element [ start, end(inclusive) ].
     */
    public static QDataSet identifySweeps( QDataSet lds, int dir ) {

        DataSetBuilder builder= new DataSetBuilder( 2, 100, 2 );
        DataSetBuilder dep0builder= new DataSetBuilder( 2, 100, 2 );
        double slope0= lds.value(1) - lds.value(0);
        int start=0;
        int end=0; // index of the right point of slope0.

        QDataSet wds= SemanticOps.weightsDataSet(lds);
        ArrayDataSet tds= ArrayDataSet.copy( SemanticOps.xtagsDataSet(lds) );
        wds= Ops.multiply( wds, SemanticOps.weightsDataSet(tds) );
        
        Datum cadence= SemanticOps.guessXTagWidth( tds, lds );
        cadence= cadence.multiply(10.0); // kludge for data on 2012-10-25
        tds.putProperty( QDataSet.CADENCE, DataSetUtil.asDataSet(cadence) );
        QDataSet nonGap= SemanticOps.cadenceCheck(tds,lds);

        for ( int i=1; i<lds.length(); i++ ) {
            double slope1=  lds.value(i) - lds.value(i-1);
            if ( slope0 * slope1 <= 0. || ( nonGap.value(i-1)==1 && nonGap.value(i)==0 ) ) {
                if ( slope0!=0. && ( dir==0 || ( slope0*dir>0 && end-start>1 ) ) ) {  //TODO: detect jumps in the LShell, e.g. only downward: \\\\\
                    if ( start<end ) {
                        builder.putValue( -1, 0, start );
                        builder.putValue( -1, 1, end-1 );
                        builder.nextRecord();
                        dep0builder.putValue( -1, 0, tds.value(start) );
                        dep0builder.putValue( -1, 1, tds.value(end-1) );
                        dep0builder.nextRecord();
                    } else {
                        // all fill.
                    }
                }
                if ( slope1!=0. || nonGap.value(i)>0 ) {
                    start= i;
                }
            } else if ( wds.value(i)>0 ) {
                if ( nonGap.value(i-1)==0 && nonGap.value(i)==1 ) {
                    start= i;
                    end= i+1;
                } else {
                    end= i+1; // end is not inclusive
                }
            }
            slope0= slope1;
        }
        
        dep0builder.putProperty( QDataSet.BINS_1, "min,maxInclusive" );
        dep0builder.putProperty( QDataSet.UNITS, SemanticOps.getUnits(tds) );
        
        DDataSet result= builder.getDataSet();
        result.putProperty( QDataSet.DEPEND_0, dep0builder.getDataSet() );
        result.putProperty( QDataSet.RENDER_TYPE, "eventsBar" );
        
        return result;
        
    }
    
    /**
     *
     * @param datax data series that is aggregation of monotonic series (up and down)
     * @param start start index limiting search
     * @param end end index inclusive limiting search
     * @param x the value to locate.
     * @param guess guess index, because we are calling this repeatedly
     * @param dir 1 if datax is increasing, -1 if decreasing
     * @return index
     */
    private static int findIndex( QDataSet datax, int start, int end, double x, int guess, int dir ) {
        int index= Math.max( Math.min( guess, end-1 ), start );
        if ( dir > 0 ) {
            while ( index<end && datax.value(index+1) < x ) index++;
            while ( index>start && datax.value(index) > x ) index--;
        } else {
            while ( index<end && datax.value(index+1) > x ) index++;
            while ( index>start && datax.value(index) < x ) index--;
        }
        return index;
    }
    
    /**
     * this is probably old and shouldn't be used.
     * @param lds
     * @param zds
     * @param start
     * @param end
     * @param col
     * @param lgrid
     * @param ds 
     */
    private static void interpolate( QDataSet lds, QDataSet zds, int start, int end, int col, QDataSet lgrid, DDataSet ds ) {
        
        Units u= SemanticOps.getUnits( ds );

        double fill= u.getFillDouble();
        ds.putProperty( QDataSet.FILL_VALUE, fill );
        
        if ( ! u.equals(  SemanticOps.getUnits( zds ) ) ) {
            throw new IllegalArgumentException("zds units must be the same as ds units!");
        }
        
        QDataSet wds= Ops.valid(zds);
        
        int index;
        
        int dir= (int) Math.signum( lds.value(end) - lds.value( start ) );
        if ( dir > 0 ) index= start; else index= end;
        
        for ( int i=0; i<lgrid.length(); i++ ) {
            double ll= lgrid.value(i);
            index= findIndex( lds, start, end, ll, index, dir );
            double alpha= ( ll - lds.value(index) ) / ( lds.value(index+1) - lds.value(index) );
            if ( alpha < 0 ) {
                ds.putValue( col, i, fill );
            } else if ( alpha > 1 ) {
                ds.putValue( col, i, fill );
            } else if ( alpha == 0 ) {
                ds.putValue( col, i, zds.value(index) );
            } else {
                if ( ( wds.value(index)==0  ) || wds.value( index+1 )==0 ) {
                    ds.putValue( col, i, fill );
                } else {
                    ds.putValue( col, i, zds.value(index) * ( 1.0 - alpha ) + zds.value(index+1) * alpha );
                }
            }
        }
        
    }
    
    private static double guessCadence( QDataSet xds, final int skip ) {
        double cadence = 0;
        
        // calculate average cadence for consistent points.  Preload to avoid extra branch.
        double cadenceS= 0;
        int cadenceN= 1;
        
        for ( int i=skip; i < xds.length(); i++) {
            double cadenceAvg;
            cadenceAvg= cadenceS/cadenceN;
            cadence = xds.value(i) - xds.value(i-skip);
            if ( cadence > 1.5 * cadenceAvg) {
                cadenceS= cadence;
                cadenceN= 1;
            } else if ( cadence < 1.5 * cadenceAvg ) {
                cadenceS+= cadence;
                cadenceN+= 1;
            }
        }
        return  cadenceS/cadenceN;
    }
    
    /**
     * rebin the datasets to rank 2 dataset ( time, LShell ), by interpolating along sweeps.  This
     * dataset has the property "sweeps", which is a dataset that indexes the input datasets.
     * @param lds rank 1 dataset of length N
     * @param zds rank 1 dataset of length N, indexed along with <tt>lds</tt>
     * @param lgrid rank 1 dataset indicating the dim 1 tags for the result dataset.
     * @return a rank 2 dataset, with one column per sweep, interpolated to <tt>lgrid</tt>
     */
    public static QDataSet rebin( QDataSet lds, QDataSet zds, QDataSet lgrid ) {
        return rebin( lds, zds, lgrid, 0 );
    }
    
    /**
     * alternate, convenient interface which where tlz is a bundle of buckshot data (T,L,Z)
     * 
     * @param tlz x,y,z (T,L,Z) bundle of Z values collected along inward and outward sweeps of L in time.
     * @param lgrid desired uniform grid of L values.
     * @param dir =1 increasing (outward) only, =-1 decreasing (inward) only, 0 both.  
     * @return a rank 2 dataset, with one column per sweep, interpolated to <tt>lgrid</tt>
     */
    public static QDataSet rebin( QDataSet tlz, QDataSet lgrid, int dir ) {
        QDataSet lds= Ops.link( Ops.slice1(tlz,0), Ops.slice1(tlz,1) );
        QDataSet zds= Ops.slice1(tlz,2);
        return rebin( lds, zds, lgrid, dir );
    }
    
    /**
     * alternate algorithm following Brian Larson's algorithm that rebin the datasets to rank 2 dataset ( time, LShell ), by interpolating along sweeps.  This
     * dataset has the property "sweeps", which is a dataset that indexes the input datasets.
     * @param lds The L values corresponding to y axis position, which should be a function of time.
     * @param tt The Time values corresponding to x axis position.  If null, then use lds.property(QDataSet.DEPEND_0).
     * @param zds the Z values corresponding to the parameter we wish to organize
     * @param tspace rank 0 cadence, such as dataset('9 hr')
     * @param lgrid rank 1 data is the grid points, such as  linspace( 2.,8.,30 )
     * @param dir =1 increasing (outward) only, =-1 decreasing (inward) only, 0 both
     * @return a rank 2 dataset, with one column per sweep, interpolated to <tt>lgrid</tt>
     */
    public static QDataSet rebin( QDataSet tt, QDataSet lds, QDataSet zds, QDataSet tspace, QDataSet lgrid, int dir ) {
        final QDataSet sweeps= identifySweeps( lds, dir );

        if ( tt==null ) {
            tt= (QDataSet) lds.property(QDataSet.DEPEND_0);
        }
        
        Units tu= SemanticOps.getUnits(tt);
        
        Number fill= (Number) zds.property( QDataSet.FILL_VALUE );
        double dfill= fill==null ? -1e31 : fill.doubleValue();

        QDataSet wds= DataSetUtil.weightsDataSet(zds);

        int ny= lgrid.length();
        double[] ss= new double[ny];
        double[] nn= new double[ny];

        DataSetBuilder builder= new DataSetBuilder( 2, 100, ny );
        DataSetBuilder tbuilder= new DataSetBuilder( 1, 100 );
        tbuilder.putProperty( QDataSet.UNITS, tu );
        
        builder.putProperty( QDataSet.FILL_VALUE,dfill );
        
        int ix=-1;
        double nextx= Double.MAX_VALUE;
        double dt= DataSetUtil.asDatum( tspace ).doubleValue( tu.getOffsetUnits() );
        double t0= -1;
        double dg= lgrid.value(1)-lgrid.value(0);
        double g0= lgrid.value(0);
        int n= ny-1;
        for ( int i=2; i<n; i++ ) {
           //if ( (lgrid.value(i)-lgrid.value(i-1)/dg ) throw new IllegalArgumentException( "lgrid must be uniform linear" );
        }
        for ( int i=0; i<sweeps.length(); i++ ) {
            int ist= (int)sweeps.value(i,0);
            int ien= (int)sweeps.value(i,1);
            for ( int j= ist; j<ien; j++ ) {
                if ( ix==-1 || tt.value(j)-nextx >= 0  ) { //reset the accumulators
                    nextx= dt * Math.ceil( tt.value(j) / dt );  // next threshold
                    if ( ix>-1 ) {
                        for ( int k=0; k<ny; k++ ) {
                            if ( nn[k]==0 ) {
                                builder.putValue( -1, k, dfill );
                            } else {
                                builder.putValue( -1, k, ss[k]/nn[k] );
                            }
                        }
                        builder.nextRecord();
                        tbuilder.putValue( -1, t0-dt/2 );
                        tbuilder.nextRecord();
                    }
                    t0= nextx;
                    
                    for ( int k=0; k<ny; k++ ) {
                        ss[k]=0;
                        nn[k]=0;
                    }

                    ix++;
                }
                int iy= (int)Math.floor( ( lds.value(j) - g0 ) / dg );
                if ( iy>=0 && iy<ny ) {
                    double w= wds.value(j);
                    ss[iy]+= w*zds.value(j);
                    nn[iy]+= w;
                }
            }
        }

        DDataSet result= builder.getDataSet();
        result.putProperty(QDataSet.DEPEND_0,tbuilder.getDataSet());
        result.putProperty(QDataSet.DEPEND_1,lgrid);
        DataSetUtil.copyDimensionProperties( zds, result );

        String title= (String) result.property( QDataSet.TITLE );
        if ( title==null ) title="";
        if ( dir==-1 ) {
            title+= "!cinward";
        } else if ( dir==1 ) {
            title+= "!coutward";
        }
        result.putProperty( QDataSet.TITLE, title );

        result.putProperty( QDataSet.FILL_VALUE, dfill );
        return result;

    }
    /**
     * rebin the datasets to rank 2 dataset ( time, LShell ), by interpolating along sweeps.  This
     * dataset has the property "sweeps", which is a dataset that indexes the input datasets.
     * @param lds rank 1 dataset of length N
     * @param zds rank 1 dataset of length N, indexed along with <tt>lds</tt>
     * @param lgrid rank 1 dataset indicating the dim 1 tags for the result dataset.
     * @param dir =1 increasing (outward) only, =-1 decreasing (inward) only, 0 both
     * @return a rank 2 dataset, with one column per sweep, interpolated to <tt>lgrid</tt>
     */
    public static QDataSet rebin( QDataSet lds, QDataSet zds, QDataSet lgrid, int dir ) {

        final QDataSet sweeps= identifySweeps( lds, dir );
        
        DDataSet result= DDataSet.createRank2( sweeps.length(), lgrid.length() );
        result.putProperty( QDataSet.UNITS, zds.property( QDataSet.UNITS ) );
        
        for ( int i=0; i<sweeps.length(); i++ ) {
            interpolate( lds, zds, (int) sweeps.value( i,0 ), (int) sweeps.value( i,1 ), i, lgrid, result );
        }
        
        DDataSet xtags= DDataSet.createRank2( sweeps.length(), 2 );
        
        QDataSet dep0= (QDataSet) lds.property(QDataSet.DEPEND_0);
        if ( dep0!=null ) {
            
            for ( int i=0; i<sweeps.length(); i++ ) {
                xtags.putValue( i,0,dep0.value( (int)sweeps.value( i,0 ) ) );
                xtags.putValue( i,1,dep0.value( (int)sweeps.value( i,1 ) ) );
            }

            DataSetUtil.putProperties( DataSetUtil.getDimensionProperties(dep0,null), xtags );
            Units xunits= SemanticOps.getUnits(dep0);
            xtags.putProperty( QDataSet.UNITS, xunits );
            xtags.putProperty( QDataSet.BINS_1, QDataSet.VALUE_BINS_MIN_MAX );
            xtags.putProperty(QDataSet.MONOTONIC, org.das2.qds.DataSetUtil.isMonotonic(dep0) );

        } else {
            for ( int i=0; i<sweeps.length(); i++ ) {
                xtags.putValue( i,0, sweeps.value( i,0 ) );
                xtags.putValue( i,1, sweeps.value( i,1 ) );
            }
        }

        Map<String,Object> userProps= new LinkedHashMap();
        userProps.put(USER_PROP_SWEEPS, sweeps);

        result.putProperty( QDataSet.USER_PROPERTIES, userProps );

        if ( lgrid.property(QDataSet.UNITS)==null ) {
            ArrayDataSet lgridCopy= ArrayDataSet.copy(lgrid); // often linspace( 4, 10, 30 ) is used to specify locations.  Go ahead and copy over units.
            Units u= (Units) lds.property( QDataSet.UNITS );
            if ( u!=null ) lgridCopy.putProperty( QDataSet.UNITS, u );
            lgrid= lgridCopy;
        }
        
        result.putProperty( QDataSet.DEPEND_1, lgrid );
        result.putProperty( QDataSet.DEPEND_0, xtags );

        DataSetUtil.putProperties( DataSetUtil.getDimensionProperties( zds, null ), result );
        
        return result;
    }
}