package yi.dp;

import yi.util.io.*;
import yi.iu.processor.filter.GaussianSmooth;
import java.util.*;

public class Histogram
{
    public static int[] getLocalMaxima(int[] s, double maximaSize, double maximaThreshold)
    {
        int size= s.length;
        int[] rst= new int[size];
        int len=(int)(size*maximaSize);
        if (len<1) len=1;
        // alternatively, we can let len=1
        // len=1;
        for (int i=len;i<size-len;i++) {
            if (s[i]==0) continue;

            rst[i]=1;
            for (int j=i-len;j<=i+len;j++) {
                if (s[j]>s[i]) {
                    rst[i]=0;
                    break;
                }
            }

            // if some neighbours of the maxima are equal, then only select the middle one
            if (rst[i]==0) continue;
            int l=0,r=0;
            for (int k=1;k<=len;k++)
            	if (s[i-k]==s[i]) l++;
            for (int k=1;k<=len;k++)
            	if (s[i+k]==s[i]) r++;
            if (Math.abs(l-r)>1) s[i]=0;
            if (l>r) s[i]=0;

        }
        // eliminate the too small maxima
        int max=0;
        for (int i=len;i<size-len;i++) {
        	if (s[i]>max) max=s[i];
        }
        for (int i=len;i<size-len;i++) {
        	if (s[i]<max*maximaThreshold) rst[i]=0;
        }

        return rst;
    }

    //clustering a histogram using suspension bridge algorithm
    //input:
    //  towers are the histogram data
    //  Rssp: the suspension ability of a tower. range:(0,1]. the less the vaule, the wider, the range
    //  Tclass: the size threshold of a cluster [0,1].  1: must including all the values; 0: not at all
    //output:
    //  the returned vector includes a set of cluster range.
    public static Vector suspensionBridgeClustering(double[] towers, double Rssp, double Tcluster)
    {
        int i;
        Vector clusters= new Vector();

        if (towers==null) return clusters;

        //average value of the bins
        double sum=0;
        for (i=0;i<towers.length;i++) sum+= towers[i];
        double average=sum/towers.length;

        //calculate the suspension ration
        double suspensionRatio= average*Rssp;

        //distribute the bins
        int[] begin= new int[towers.length];
        int[] end= new int[towers.length];

        for (i=0;i<towers.length;i++) {
            if (towers[i]>0) {
                int suspensionWidth=(int)Math.ceil(towers[i]/suspensionRatio);
                if (suspensionWidth<20) suspensionWidth=20;
                if (suspensionWidth>75) suspensionWidth=75;
                int left_arm= (suspensionWidth-1)/2;
                int right_arm= suspensionWidth-1-left_arm;
                int left=Math.max(0,i-left_arm);
                int right=Math.min(towers.length-1,i+right_arm);
                begin[left]++;
                end[right]++;
            }
        }

        //re-calculate the suspended histogram
        int now=0;
        int[] sh= new int[towers.length];
        for (i=0;i<towers.length;i++) {
            now+=begin[i];
            if (i>0) now-=end[i-1];
            sh[i]=now;
        }

        //System.err.println("average="+average);
        //OutputResult.outputResult("towers", towers);
        //OutputResult.outputResult("suspended", sh);

        //find the cluster candidates in the suspended
        Vector clusterCans= new Vector();
        for (int l=0;l<sh.length;l++) {
            if (sh[l]==0) continue;          //find the left
            int r=l;
            while (r<sh.length-1&&sh[r+1]!=0) r++;      //find the right

            int[] c= {l,r};
            clusterCans.addElement(c);
            l=r+1;
        }

        //calculate the sum of the cluster candidates
        double[] sumRange= new double[clusterCans.size()];
        for (i=0;i<clusterCans.size();i++) {
            int[] range= (int[])clusterCans.elementAt(i);
            for (int j=range[0];j<=range[1];j++) {
                sumRange[i]+=towers[j];
            }
        }

        //find the range with the max sum
        for (i=0;i<sumRange.length;i++) {
            if (sumRange[i]>=sum*Tcluster) {
                clusters.addElement(clusterCans.elementAt(i));
            }
        }
        return clusters;

    }

    /** This function breaks a histogram with tolerance
      * @author Yi Li
      * @version 1.0.0
      * @param h The histogram to be processed
      * @param tolerance The tolerancable width of a break (in other words,
      *     only when the break is greater than the tolerance, it is considered as a break
      * @param Tcluster The threshold of a cluster.  After breaking the histogram,
      *     only those clusters taking a greater portion than Tcluster are considered.
      *     The range of this argument is [0..1].  1 means the cluster must contains all the
      *     value of the histogram; 0 means don't care
      * @return A Vector of int[] which is the range of the clusters
      */
    public static Vector breakWithTolerance(int[] h, double tolerance, double Tcluster)
    {
        int i;

        //find the possible break
        //OutputResult.outputResult("h", h, 1, h.length);
        Vector breaks= new Vector();
        for (int l=0;l<h.length;l++) {
            if (h[l]!=0) continue;       //find a break - left
            int r=l;
            while (r<h.length-1&&h[r+1]==0) r++;    //find the right of the break
            //ingore those breaks tolerancable
            if (r-l+1<=tolerance) continue;
            int[] b= {l,r};
            breaks.addElement(b);
            l=r+1;
        }

        Vector clusters= new Vector();

        //get the cluster candidates
        Vector clusterCans= new Vector();
        int l=0;    //left boundary of the cluster
        for (i=0;i<breaks.size();i++) {
            int[] b= (int[]) breaks.elementAt(i);
            if (b[0]-1<=l)      //the range is [l..b[i]-1], so ..
            {
                l=b[1]+1;
                continue;
            }
            int[] range= {l,b[0]-1};
            clusterCans.addElement(range);
            l=b[1]+1;                   //the left boundary of next range
        }
        int[] range= {l,h.length-1};    //the last one
        clusterCans.addElement(range);

        //calculate the sum of the ranges
        double sum=0;
        for (i=0;i<h.length;i++) sum+=h[i];
        double threshold= sum*Tcluster;
        Vector sumRanges= new Vector();
        for (i=0;i<clusterCans.size();i++) {
            range= (int[])clusterCans.elementAt(i);
            double sumRange=0;
            for (int j=range[0];j<=range[1];j++) sumRange+=h[j];
            if (sumRange>=threshold) {
                sumRanges.addElement(new Double(sumRange));
                clusters.addElement(range);
            }
        }

        return clusters;
    }

    public static Vector breakWithTolerance(int[] h,
        double tolerance, double Tcluster, double Tbreak)
    {
        int i;

        //min height of the bin
        int minheight=0,count=0;
        for (i=0;i<h.length;i++) {
            if (h[i]>0) {
                minheight+=h[i];
                count++;
            }
        }
        minheight=(int)((double)minheight*Tbreak/count);
        if (minheight==0) {
            System.err.println("minheight=0-->1");
            minheight=1;
        }
        //find the possible break
        //OutputResult.outputResult("h", h, 1, h.length);
        Vector breaks= new Vector();
        for (int l=0;l<h.length;l++) {
            if (h[l]>=minheight) continue;       //find a break - left
            int r=l;
            while (r<h.length-1&&h[r+1]<minheight) r++;    //find the right of the break
            //ingore those breaks tolerancable
            if (r-l+1<=tolerance) continue;
            int[] b= {l,r};
            breaks.addElement(b);
            l=r+1;
        }

        Vector clusters= new Vector();

        //get the cluster candidates
        Vector clusterCans= new Vector();
        int l=0;    //left boundary of the cluster
        for (i=0;i<breaks.size();i++) {
            int[] b= (int[]) breaks.elementAt(i);
            if (b[0]-1<=l)      //the range is [l..b[i]-1], so ..
            {
                l=b[1]+1;
                continue;
            }
            int[] range= {l,b[0]-1};
            clusterCans.addElement(range);
            l=b[1]+1;                   //the left boundary of next range
        }
        int[] range= {l,h.length-1};    //the last one
        clusterCans.addElement(range);

        //calculate the sum of the ranges
        double sum=0;
        for (i=0;i<h.length;i++) sum+=h[i];
        double threshold= sum*Tcluster;
        Vector sumRanges= new Vector();
        for (i=0;i<clusterCans.size();i++) {
            range= (int[])clusterCans.elementAt(i);
            double sumRange=0;
            for (int j=range[0];j<=range[1];j++) sumRange+=h[j];
            if (sumRange>=threshold) {
                sumRanges.addElement(new Double(sumRange));
                clusters.addElement(range);
            }
        }

        return clusters;
    }
    /**
     * This function find the position of the peaks (the tip of the peak)
     *
     * @param data
     * @return Return peak[3][len_of_histogram].
     * If there is a peak on histogram[i], then
     * peak[0][i]=1, peak[1][i]=left_boundary, peak[2][i]=right_boundary
     */
    public static int WT=72;       //argument used in peak finding, if the peak is broader than WT, it is not considered to be prominent
                            //a multiplicative weight WB is included in the salience measure
    public static double WB=0.1;   //see comments with WT
    public static double GARMA_THRSH=3; //the threshold for garma (garma is the result of peak salience measure)


    public static int[][] findPeaks(double[] data)
    {
        return findPeaks(data,GARMA_THRSH);
    }

    public static int[][] findPeaks(double[] data, double garma_thrsh)
    {
        int len= data.length;
        double[] diff= new double[len];
        int i;
        int[][] peaks= new int[3][len];
        double[][] result=null;

        boolean TRAINING=false;
        if (TRAINING) {
            result= new double[4][];
            result[0]= data;
            result[1]= new double[len];
            result[2]= new double[len];
            result[3]= new double[len];
        }

        //first difference of orientation histogram
        for (i=0;i<len-1;i++) diff[i]= data[i+1]-data[i];
        diff[len-1]= data[0]-data[len-1];

        //find the "average" of the lower half part of the histogram
        //1. sort the histogram first
        double[] copy= new double[len]; //make a copy of the data
        for (i=0;i<len;i++) copy[i]= data[i];
        double sum=0;
        for (i=0;i<len/2;i++) {
            for (int j=i+1;j<len;j++) {
                if (copy[i]>copy[j]) {
                    double temp= copy[i]; copy[i]= copy[j]; copy[j]=temp;
                }
            }
            sum+= copy[i];
        }
        double average= sum/(len/2);


        //to find the peaks, at that point the difference changes from positive to negative
        //eg data[1 4 5 2] -> diff[3 1 -3 -1]
        //first find a index at which the diff changes from positive to negative
        //then search diff backward to find a index l, satisfying diff[l-1]<diff[l]>diff[l+1];
        //also search diff forward to find a index r, satisfying diff[r-1]>diff[r]<diff[r+1];
        for (i=0;i<len;i++) {
            //REMEMBER (-1)%5=-1 !!!!
            //System.err.println("-1%5="+(-1)%5);
            if (diff[i]>0||diff[(i-1+len)%len]<=0) continue;
            //find a zero-crossing point
            //diff[i] satisfy diff[i]<=0 && diff[i-1]>0

            if (data[i]<=average) continue;
            //now the peak must be higher than the average

            //find the local maxima in the diff histogram
            int l=(i-1+len)%len;
            double max= diff[l];    //remember the maxima
            int max_l=l;            //remember the index of the maxima
            do {
                if (diff[l]>max) {
                    max=diff[l]; max_l=l;
                }
                l=(l-1+len)%len;
            } while (diff[l]>0);
            //diff[l]<=0, data[l+1] is a minima
            //eg data{3,0,3} -> diff{-3,3,?} -> diff[l]=-3, data[l+1]:min
            l=(l+1)%len;
            //now data[l] is the minima
            //data[max_l+1] and data[max_l] are around the inflection point

            //find the local minima in the diff histogram
            int r=i;
            double min= diff[r];
            int min_r=r;
            do {
                if (diff[r]<min) {
                    min=diff[r]; min_r=r;
                }
                r=(r+1)%len;
            } while (diff[r]<0);
            //diff[r]>=0, data[r] is a minima
            //eg data[3,0,3] -> diff[-3,3,?] -> diff[r]=3, data[r]:min

            //find the left boundary of the peak
            //calculate the average of the diff to the left of the peak
            double avrg_d= ((double)(data[i]-data[l]))/((i-l+len)%len);
            int lb=l;
            while (diff[lb]<avrg_d) lb=(lb+1)%len;
            lb=(lb+1)%len;
            //data[lb] is in the range of the peak, and data[lb-1] is not

            //find the right boundary of the peak
            avrg_d= -((double)(data[i]-data[r]))/((r-i+len)%len);
            int rb=r;
            while (diff[rb]>avrg_d) rb=(rb-1+len)%len;
            //data[rb] is in the range of the peak, and data[rb+1] is not

            //calculate the garma
            double hs= data[i];
            //double dv= diff[max_l]-diff[min_r];
            double dv= max-min;
            double dh= (min_r-max_l+len)%len;
            double max_val= data[i];    //maximum point
            double min_val= data[l]>data[r]?data[l]:data[r];    //the largest minimum value in data
            double wm= 1-(min_val/max_val);
            //double wb= ((min_r-max_l+len)%len)>WT?WB:1;
            double board= (rb-lb+len)%len;
            double wb= board>WT?WB:1;
            double garma= hs*dv/dh*wm*wb;
            if (TRAINING) {
                result[1][lb]=max_val;
                result[2][rb]=max_val;
                result[3][i]=garma;
            }
            //if (garma>0.1) System.err.println("garma="+garma);
            if (garma>garma_thrsh) {
                peaks[0][i]=1;
                peaks[1][i]=lb;
                peaks[2][i]=rb;
            }
            //because the r is already a minimum, so the next maxima must be the right to r
            if (r>i) i=r+1;
            else break;
        }
        //if (TRAINING) yi.util.io.OutputResult.outputResult("peaks",peaks);
        if (TRAINING) yi.util.io.OutputResult.outputResult("h-peaks",result);

        return peaks;

    }

    //the boundaries are at the local minima
    public static int[][] findFatPeaks(double[] data, double garma_thrsh)
    {
        int len= data.length;
        double[] diff= new double[len];
        int i;
        int[][] peaks= new int[3][len];
        double[][] result=null;

        boolean TRAINING=false;
        if (TRAINING) {
            result= new double[4][];
            result[0]= data;
            result[1]= new double[len];
            result[2]= new double[len];
            result[3]= new double[len];
        }

        //first difference of orientation histogram
        for (i=0;i<len-1;i++) diff[i]= data[i+1]-data[i];
        diff[len-1]= data[0]-data[len-1];

        //find the "average" of the lower half part of the histogram
        //1. sort the histogram first
        double[] copy= new double[len]; //make a copy of the data
        for (i=0;i<len;i++) copy[i]= data[i];
        double sum=0;
        for (i=0;i<len/2;i++) {
            for (int j=i+1;j<len;j++) {
                if (copy[i]>copy[j]) {
                    double temp= copy[i]; copy[i]= copy[j]; copy[j]=temp;
                }
            }
            sum+= copy[i];
        }
        double average= sum/(len/2);


        //to find the peaks, at that point the difference changes from positive to negative
        //eg data[1 4 5 2] -> diff[3 1 -3 -1]
        //first find a index at which the diff changes from positive to negative
        //then search diff backward to find a index l, satisfying diff[l-1]<diff[l]>diff[l+1];
        //also search diff forward to find a index r, satisfying diff[r-1]>diff[r]<diff[r+1];
        for (i=0;i<len;i++) {
            //REMEMBER (-1)%5=-1 !!!!
            //System.err.println("-1%5="+(-1)%5);
            if (diff[i]>0||diff[(i-1+len)%len]<=0) continue;
            //find a zero-crossing point
            //diff[i] satisfy diff[i]<=0 && diff[i-1]>0

            if (data[i]<=average) continue;
            //now the peak must be higher than the average

            //find the local maxima in the diff histogram
            int l=(i-1+len)%len;
            double max= diff[l];    //remember the maxima
            int max_l=l;            //remember the index of the maxima
            do {
                if (diff[l]>max) {
                    max=diff[l]; max_l=l;
                }
                l=(l-1+len)%len;
            } while (diff[l]>0);
            //diff[l]=<0, data[l+1] is a minima
            //eg data{3,0,3} -> diff{-3,3,?} -> diff[l]=-3, data[l+1]:min
            l=(l+1)%len;
            //now data[l] is the minima
            //data[max_l+1] and data[max_l] are around the inflection point

            //find the local minima in the diff histogram
            int r=i;
            double min= diff[r];
            int min_r=r;
            do {
                if (diff[r]<min) {
                    min=diff[r]; min_r=r;
                }
                r=(r+1)%len;
            } while (diff[r]<0);
            //diff[r]>=0, data[r] is a minima
            //eg data[3,0,3] -> diff[-3,3,?] -> diff[r]=3, data[r]:min

            //find the left boundary of the peak
            //calculate the average of the diff to the left of the peak
            //double avrg_d= ((double)(data[i]-data[l]))/((i-l+len)%len);
            int lb=l;
            //while (diff[lb]<avrg_d) lb=(lb+1)%len;
            //lb=(lb+1)%len;
            //data[lb] is in the range of the peak, and data[lb-1] is not

            //find the right boundary of the peak
            //avrg_d= -((double)(data[i]-data[r]))/((r-i+len)%len);
            int rb=r;
            //while (diff[rb]>avrg_d) rb=(rb-1+len)%len;
            //data[rb] is in the range of the peak, and data[rb+1] is not

            //calculate the garma
            double hs= data[i];
            //double dv= diff[max_l]-diff[min_r];
            double dv= max-min;
            double dh= (min_r-max_l+len)%len;
            double max_val= data[i];    //maximum point
            double min_val= data[l]>data[r]?data[l]:data[r];    //the largest minimum value in data
            double wm= 1-(min_val/max_val);
            //double wb= ((min_r-max_l+len)%len)>WT?WB:1;
            double board= (rb-lb+len)%len;
            double wb= board>WT?WB:1;
            double garma= hs*dv/dh*wm*wb;
            if (TRAINING) {
                result[1][lb]=max_val;
                result[2][rb]=max_val;
                result[3][i]=garma;
            }
            //if (garma>0.1) System.err.println("garma="+garma);
            if (garma>garma_thrsh) {
                peaks[0][i]=1;
                peaks[1][i]=lb;
                peaks[2][i]=rb;
            }
            //because the r is already a minimum, so the next maxima must be the right to r
            if (r>i) i=r+1;
            else break;
        }
        if (TRAINING) yi.util.io.OutputResult.outputResult("peaks",peaks);
        //if (TRAINING) yi.util.io.OutputResult.outputResult("h-result",result);
        if (TRAINING) {
            double[][] rst_t= new double[result.length][result[0].length];
            for (i=0;i<result.length;i++)
                for (int j=0;j<result[0].length;j++) {
                    rst_t[i][j]=result[i][result[0].length-1-j];
                }
            yi.util.io.OutputResult.outputResult("h-result",rst_t);
        }

        return peaks;

    }
    //the boundaries are at the local minima
    public static int[][] findFatPeaks_no_period(double[] data, double garma_thrsh)
    {
        int len= data.length;
        double[] diff= new double[len];
        int i;
        int[][] peaks= new int[3][len];
        double[][] result=null;

        boolean TRAINING=false;
        if (TRAINING) {
            result= new double[4][];
            result[0]= data;
            result[1]= new double[len];
            result[2]= new double[len];
            result[3]= new double[len];
        }

        //first difference of orientation histogram
        for (i=0;i<len-1;i++) diff[i]= data[i+1]-data[i];
        diff[len-1]= 0;
        diff[0]=0;

        //find the "average" of the lower half part of the histogram
        //1. sort the histogram first
        double[] copy= new double[len]; //make a copy of the data
        for (i=0;i<len;i++) copy[i]= data[i];
        double sum=0;
        for (i=0;i<len/2;i++) {
            for (int j=i+1;j<len;j++) {
                if (copy[i]>copy[j]) {
                    double temp= copy[i]; copy[i]= copy[j]; copy[j]=temp;
                }
            }
            sum+= copy[i];
        }
        double average= sum/(len/2);


        //to find the peaks, at that point the difference changes from positive to negative
        //eg data[1 4 5 2] -> diff[3 1 -3 -1]
        //first find a index at which the diff changes from positive to negative
        //then search diff backward to find a index l, satisfying diff[l-1]<diff[l]>diff[l+1];
        //also search diff forward to find a index r, satisfying diff[r-1]>diff[r]<diff[r+1];
        for (i=1;i<len;i++) {
            //REMEMBER (-1)%5=-1 !!!!
            //System.err.println("-1%5="+(-1)%5);
            if (diff[i]>0||diff[i-1]<=0) continue;
            //find a zero-crossing point
            //diff[i] satisfy diff[i]<=0 && diff[i-1]>0

            if (data[i]<=average) continue;
            //now the peak must be higher than the average

            //find the local maxima in the diff histogram
            int l=i-1;
            double max= diff[l];    //remember the maxima
            int max_l=l;            //remember the index of the maxima
            do {
                if (diff[l]>max) {
                    max=diff[l]; max_l=l;
                }
                l=l-1;
            } while (l>0&&diff[l]>0);
            //diff[l]=<0, data[l+1] is a minima
            //eg data{3,0,3} -> diff{-3,3,?} -> diff[l]=-3, data[l+1]:min
            l=l+1;
            //now data[l] is the minima
            //data[max_l+1] and data[max_l] are around the inflection point

            //find the local minima in the diff histogram
            int r=i;
            double min= diff[r];
            int min_r=r;
            do {
                if (diff[r]<min) {
                    min=diff[r]; min_r=r;
                }
                r=r+1;
            } while (r<diff.length&&diff[r]<0);
            if (r==diff.length) r=diff.length-1;
            //diff[r]>=0, data[r] is a minima
            //eg data[3,0,3] -> diff[-3,3,?] -> diff[r]=3, data[r]:min

            //find the left boundary of the peak
            //calculate the average of the diff to the left of the peak
            //double avrg_d= ((double)(data[i]-data[l]))/((i-l+len)%len);
            int lb=l;
            //while (diff[lb]<avrg_d) lb=(lb+1)%len;
            //lb=(lb+1)%len;
            //data[lb] is in the range of the peak, and data[lb-1] is not

            //find the right boundary of the peak
            //avrg_d= -((double)(data[i]-data[r]))/((r-i+len)%len);
            int rb=r;
            //while (diff[rb]>avrg_d) rb=(rb-1+len)%len;
            //data[rb] is in the range of the peak, and data[rb+1] is not

            //calculate the garma
            double hs= data[i];
            //double dv= diff[max_l]-diff[min_r];
            double dv= max-min;
            double dh= (min_r-max_l+len)%len;
            double max_val= data[i];    //maximum point
            double min_val= data[l]>data[r]?data[l]:data[r];    //the largest minimum value in data
            double wm= 1-(min_val/max_val);
            //double wb= ((min_r-max_l+len)%len)>WT?WB:1;
            double board= (rb-lb+len)%len;
            double wb= board>WT?WB:1;
            double garma= hs*dv/dh*wm*wb;
            if (TRAINING) {
                result[1][lb]=max_val;
                result[2][rb]=max_val;
                result[3][i]=garma;
            }
            //if (garma>0.1) System.err.println("garma="+garma);
            if (garma>garma_thrsh) {
                peaks[0][i]=1;
                peaks[1][i]=lb;
                peaks[2][i]=rb;
            }
            //because the r is already a minimum, so the next maxima must be the right to r
            if (r>i) i=r+1;
            else break;
        }
        if (TRAINING) yi.util.io.OutputResult.outputResult("peaks",peaks);
        //if (TRAINING) yi.util.io.OutputResult.outputResult("h-result",result);
        if (TRAINING) {
            double[][] rst_t= new double[result.length][result[0].length];
            for (i=0;i<result.length;i++)
                for (int j=0;j<result[0].length;j++) {
                    rst_t[i][j]=result[i][result[0].length-1-j];
                }
            yi.util.io.OutputResult.outputResult("h-result",rst_t);
        }

        return peaks;

    }

    public static int[] findHighestPeak_discrete(double[] data, double sigma, double garma_thrsh)
    {
        double[] smooth= (new GaussianSmooth(sigma)).smoothing(data);
        return findHighestPeak(smooth, garma_thrsh);
    }


    public static int[] findHighestPeak(double[] data, double garma_thrsh)
    {
        int len= data.length;
        double[] diff= new double[len];
        int i;
        int[][] peaks= new int[3][len];
        double[][] result=null;

        boolean TRAINING=true;
        if (TRAINING) {
            result= new double[4][];
            result[0]= data;
            result[1]= new double[len];
            result[2]= new double[len];
            result[3]= new double[len];
        }

        //first difference of orientation histogram
        for (i=0;i<len-1;i++) diff[i]= data[i+1]-data[i];
        diff[len-1]= data[0]-data[len-1];

        //find the "average" of the lower half part of the histogram
        //1. sort the histogram first
        double[] copy= new double[len]; //make a copy of the data
        for (i=0;i<len;i++) copy[i]= data[i];
        double sum=0;
        for (i=0;i<len/2;i++) {
            for (int j=i+1;j<len;j++) {
                if (copy[i]>copy[j]) {
                    double temp= copy[i]; copy[i]= copy[j]; copy[j]=temp;
                }
            }
            sum+= copy[i];
        }
        double average= sum/(len/2);


        //to find the peaks, at that point the difference changes from positive to negative
        //eg data[1 4 5 2] -> diff[3 1 -3 -1]
        //first find a index at which the diff changes from positive to negative
        //then search diff backward to find a index l, satisfying diff[l-1]<diff[l]>diff[l+1];
        //also search diff forward to find a index r, satisfying diff[r-1]>diff[r]<diff[r+1];
        for (i=0;i<len;i++) {
            //REMEMBER (-1)%5=-1 !!!!
            //System.err.println("-1%5="+(-1)%5);
            if (diff[i]>0||diff[(i-1+len)%len]<=0) continue;
            //find a zero-crossing point
            //diff[i] satisfy diff[i]<=0 && diff[i-1]>0

            if (data[i]<=average) continue;
            //now the peak must be higher than the average

            //find the local maxima in the diff histogram
            int l=(i-1+len)%len;
            double max= diff[l];    //remember the maxima
            int max_l=l;            //remember the index of the maxima
            do {
                if (diff[l]>max) {
                    max=diff[l]; max_l=l;
                }
                l=(l-1+len)%len;
            } while (diff[l]>=0);
            //diff[l]<0, data[l+1] is a minima
            //eg data{3,0,3} -> diff{-3,3,?} -> diff[l]=-3, data[l+1]:min
            l=(l+1)%len;
            //now data[l] is the minima
            //data[max_l+1] and data[max_l] are around the inflection point

            //find the local minima in the diff histogram
            int r=i;
            double min= diff[r];
            int min_r=r;
            do {
                if (diff[r]<min) {
                    min=diff[r]; min_r=r;
                }
                r=(r+1)%len;
            } while (diff[r]<=0);
            //diff[r]>0, data[r] is a minima
            //eg data[3,0,3] -> diff[-3,3,?] -> diff[r]=3, data[r]:min

            //find the left boundary of the peak
            //calculate the average of the diff to the left of the peak
            double avrg_d= ((double)(data[i]-data[l]))/((i-l+len)%len);
            int lb=l;
            while (diff[lb]<avrg_d) lb=(lb+1)%len;
            lb=(lb+1)%len;
            //data[lb] is in the range of the peak, and data[lb-1] is not

            //find the right boundary of the peak
            avrg_d= -((double)(data[i]-data[r]))/((r-i+len)%len);
            int rb=r;
            while (diff[rb]>avrg_d) rb=(rb-1+len)%len;
            //data[rb] is in the range of the peak, and data[rb+1] is not

            //calculate the garma
            double hs= data[i];
            //double dv= diff[max_l]-diff[min_r];
            double dv= max-min;
            double dh= (min_r-max_l+len)%len;
            double max_val= data[i];    //maximum point
            double min_val= data[l]>data[r]?data[l]:data[r];    //the largest minimum value in data
            double wm= 1-(min_val/max_val);
            //double wb= ((min_r-max_l+len)%len)>WT?WB:1;
            double board= (rb-lb+len)%len;
            double wb= board>WT?WB:1;
            double garma= hs*dv/dh*wm*wb;
            if (TRAINING) {
                result[1][lb]=max_val;
                result[2][rb]=max_val;
                result[3][i]=garma;
            }
            //if (garma>0.1) System.err.println("garma="+garma);
            if (garma>garma_thrsh) {
                peaks[0][i]=1;
                peaks[1][i]=lb;
                peaks[2][i]=rb;
            }
            //because the r is already a minimum, so the next maxima must be the right to r
            if (r>i) i=r+1;
            else break;
        }
        //if (TRAINING) yi.util.io.OutputResult.outputResult("peaks",peaks);
        if (TRAINING) yi.util.io.OutputResult.outputResult("h-peaks",result);

        //find the highest peak
        int[][] ranges=peaks2Ranges(peaks);
        if (ranges==null||ranges.length==0) return null;
        double[] sums= new double[ranges.length];
        for (i=0;i<data.length;i++) {
            for (int j=0;j<ranges.length;j++) {
                if (i>=ranges[j][0]&&i<=ranges[j][1]) {
                    sums[j]+=data[i];
                    break;
                }
            }
        }
        double max=0;
        int maxI=0;
        for (i=0;i<sums.length;i++) {
            if (sums[i]>max) {
                max= sums[i];
                maxI=i;
            }
        }
        return ranges[maxI];

    }

    /**
     * Convert the peaks array to ranges array
     *
     * @param peaks peaks[0][i]=1  there is a peak at degree i;
     *            =0  else
     * peaks[1][i]    the left boundary of the range, inclusively
     *     //  peaks[2][i]    the right boundary of the range, inclusively
     * @return range[i][0]    the left boundary of the i-th peak
     * range[i][1]    the right boundary of the i-th peak
     */
    public static int[][] peaks2Ranges(int[][] peaks)
    {
        int count=0;
        for (int i=0;i<peaks[0].length;i++) if (peaks[0][i]>0) count++;
        int[][] ranges= new int[count][2];
        //System.err.println("found "+count+" ranges");
        count=0;
        for (int i=0;i<peaks[0].length;i++)
            if (peaks[0][i]>0) {
                ranges[count][0]=peaks[1][i];
                ranges[count][1]=peaks[2][i];
                //System.err.println("  "+count+": "+ranges[count][0]+"-"+ranges[count][1]);
                count++;
            }
        return ranges;
    }


    //get the peaks in the histogram
    //only the break more than tolerance are considered as break
    //the peak must have more than TSize value
    public static Vector getPeaksSimple(int[] h, int tolerance, int TSize)
    {
        int i;

        //find the possible break
        Vector breaks= new Vector();
        for (int l=0;l<h.length;l++) {
            if (h[l]!=0) continue;       //find a break - left
            int r=l;
            while (r<h.length-1&&h[r+1]==0) r++;    //find the right of the break
            //ingore those breaks tolerancable
            if (r-l+1<=tolerance) continue;
            int[] b= {l,r};
            breaks.addElement(b);
            l=r+1;
        }

        Vector clusters= new Vector();

        //get the cluster candidates
        Vector clusterCans= new Vector();
        int l=0;    //left boundary of the cluster
        for (i=0;i<breaks.size();i++) {
            int[] b= (int[]) breaks.elementAt(i);
            if (b[0]-1<l)  //the range is [l..b[i]-1], so ..
            {
                l=b[1]+1;
                continue;
            }
            int[] range= {l,b[0]-1};
            clusterCans.addElement(range);
            l=b[1]+1;                   //the left boundary of next range
        }
        if (l<=h.length-1) {
            int[] range= {l,h.length-1};    //the last one
            clusterCans.addElement(range);
        }

        //calculate the sum of the ranges
        for (i=0;i<clusterCans.size();i++) {
            int[] range= (int[])clusterCans.elementAt(i);
            int sumRange=0;
            for (int j=range[0];j<=range[1];j++) sumRange+=h[j];
            if (sumRange>=TSize) {
                //sumRanges.addElement(new Integer(sumRange));
                clusters.addElement(range);
            }
        }

        return clusters;
    }

    //get the peaks in the histogram
    //only the break more than tolerance are considered as break
    //the peak must have more than TSize value
    public static Vector getPeaksSimple_period(int[] h, int tolerance, int TSize)
    {
        int i;

        //find the possible break
        Vector breaks= new Vector();
        int first=0;
        if (h[first]==0) {
            first=h.length-1;
            while (h[first]==0&&first>0) first--;
            if (first==0) return null;
            first++;
        }

        int modul= h.length;
        for (int l=0;l<h.length;l++) {
            if (h[(first+l)%modul]!=0) continue;       //find a break - left
            int r=l;
            while (r<h.length-1&&h[(first+r+1)%modul]==0) r++;    //find the right of the break
            //ingore those breaks tolerancable
            if (r-l+1<=tolerance) continue;
            int[] b= {l,r};
            breaks.addElement(b);
            l=r+1;
        }

        Vector clusters= new Vector();

        //get the cluster candidates
        Vector clusterCans= new Vector();
        int l=0;    //left boundary of the cluster
        for (i=0;i<breaks.size();i++) {
            int[] b= (int[]) breaks.elementAt(i);
            if (b[0]-1<l)  //the range is [l..b[i]-1], so ..
            {
                l=b[1]+1;
                continue;
            }
            int[] range= {l,b[0]-1};
            clusterCans.addElement(range);
            l=b[1]+1;                   //the left boundary of next range
        }
        if (l<=h.length-1) {
            int[] range= {l,h.length-1};    //the last one
            clusterCans.addElement(range);
        }

        //calculate the sum of the ranges
        for (i=0;i<clusterCans.size();i++) {
            int[] range= (int[])clusterCans.elementAt(i);
            int sumRange=0;
            for (int j=range[0];j<=range[1];j++) sumRange+=h[(first+j)%modul];
            if (sumRange>=TSize) {
                //sumRanges.addElement(new Integer(sumRange));
                clusters.addElement(range);
            }
        }

        //adjust the range
        for (i=0;i<clusters.size();i++) {
            int[] range=(int[])clusters.elementAt(i);
            range[0]=(first+range[0])%modul;
            range[1]=(first+range[1])%modul;
        }

        return clusters;
    }

    public static int[][] getPeaksSimple_period_decent(int[] h, int tolerance, int TSize)
    {
        Vector rst= getPeaksSimple_period(h,tolerance,TSize);
        int[][] drst= new int[rst.size()][];
        for (int i=0;i<drst.length;i++)
            drst[i]= (int[]) rst.elementAt(i);
        return drst;
    }

    // size: the size of histogram
    // quan: quan
    // data[] = {10,10,20,30}, range = 4, size = 10 => hist = {0, 2, 1, 1}
    public static int[] makeHistogram(int[] data, int quan, int size) {
        int[] hist = new int[size];
        return hist;
    }

}