package yi.ai.cluster;

import yi.iu.processor.filter.*;
import yi.ds.DataRange;
import java.util.*;
import yi.stat.hist.WeightedHistogram;

/**
 * Class OneDHistBasedCluster has the algorithms to cluster 1D data using
 * histogram shape based methods.
 */
public class OneDHistBasedCluster {

  /**
   * argument used in peak finding (in method localMaximaGuidedHistClustering),
   * if the peak is broader than WT, it is not considered to be prominent.  a
   * multiplicative weight WB is included in the salience measure.  WT_RATE is
   * normalized by the length of the histogram
   */
  public static double WT_RATE = 0.4;

  /**
   * see comments with WT
   * @see #WT_RATE
   */
  public static double WB = 0.1;

  /**
   * First the histogram will be smoothed using Gaussian.  Then local maximas
   * are found.  The boundary is found by searching from the local maximas until
   * reaching a local minima.
   * Moved from yi.dp.Histogram
   * @param data
   * @param sigma the sigma for Gaussian smooth
   * @param garma_thrsh
   * @param isPeriodic the histogram is periodic.  Eg. a line orientation
   *                   histogram is periodic since direction 0 is the same as
   *                   direction 180
   * @return
   */
  public static OneDClusters localMaximaGuidedHistClustering(double[] data, double sigma, double garma_thrsh, boolean isPeriodic)
  {
      int len = data.length;
      double WT = len * WT_RATE;

      // Gaussian smooth first
      double[] smoothed= (new GaussianSmooth(sigma)).smoothing(data, true);

      // 1-order difference diff[i] = data[i+1] - data[i]
      double[] diff= new double[len];
      for (int i=0;i<len-1;i++) diff[i] = smoothed[i+1] - smoothed[i];
      if (isPeriodic) {
        diff[len-1] = smoothed[0] - smoothed[len-1];
      } else {
        diff[len-1] = diff[0] = 0;
      }

      // find the average of the lower part of the data as a threshold of the
      // local maxima (the maxima must be no less than it to be qualified)
      double average = averageOfLowerPart(smoothed, 0.5);


      // To find the peaks, first find the point at which the difference changes
      // from positive to negative:
      //
      //   eg. data[1 4 5 2] -> diff[3 1 -3 -1]
      //
      //   diff[i] < 0 && diff[i-1] > 0 => data[i] is a local maxima
      //
      // first find the 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];
      //
      Vector data_ranges = new Vector(); // store the data ranges of clusters
      for (int i = isPeriodic? 0: 1; i < len; i++) {   // searching for local maxima begins

        // NOTE: (-1)%5=-1 !!!!
        // not local maxima
        if (diff[i] > 0) continue;
        if (isPeriodic ? diff[(i-1+len)%len] <= 0 : diff[i-1] <= 0) continue;

        // now index i satifies:
        //   diff[i] satisfy diff[i]<=0 && diff[i-1]>0

        // use average of lower part as the threshold
        if (smoothed[i] <= average) continue;

        // find the local minima as the left boundary
        // diff[max_l] is the local maxima in the diff histogram
        // data[max_l+1] and data[max_l] are around the inflection point
        int l = isPeriodic ? (i - 1 + len) % len : 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 = isPeriodic ? (l - 1 + len) % len : l - 1;  // backward
        } while ((isPeriodic || 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] is min
        l = isPeriodic ? (l+1)%len : l + 1;

        // find the local minima as the right boundary
        int r = i;
        double min = diff[r];
        int min_r = r;
        do {
          if (diff[r]<min) {
            min=diff[r]; min_r=r;
          }
          r = isPeriodic ? (r+1)%len : r+1;
        } while ((isPeriodic || r < diff.length) && 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

        // the boundaries of the peak
        if (!isPeriodic && r == diff.length ) r = diff.length - 1;
        int lb=l;
        int rb=r;

        // calculate the garma (a value to reflect whether the peak is "good")
        double hs = smoothed[i];
        // double dv= diff[max_l]-diff[min_r];
        double dv = max - min;
        double dh = (min_r-max_l+len)%len;
        double max_val = smoothed[i];    //maximum point
        double min_val = (smoothed[l]>smoothed[r]) ? smoothed[l] : smoothed[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;

        // compare garma with the threshold
        if (garma > garma_thrsh) {
          data_ranges.add(new DataRange(lb, rb));
          // [lb, rb] should be the range.  Since [rb, ?] could be the next
          // cluster too, set [lb, rb) as the range
        }

        // 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;
      }

      // convert the data format
      return new OneDClusters(data_ranges, 0, null);
  }

  /**
   * Clustering a histogram using suspension bridge algorithm.  This algorihtm
   * is good for clustering sparse histogram.  A non-zero value in the histogram
   * can be seen as a tower which can grab some nearby bins.  A cluster is a
   * continuous bins that be grabbed by towers.
   * @param towers
   * @param Rssp the suspension ability of a tower. range:(0,1]. the less the
   *             vaule, the wider the range.  0.18 was used and it looked OK
   * @param Tcluster the size threshold of a cluster.  ragne: [0,1].  1: must
   *                 including all the values; 0: not at all
   * @param max_suspension_width the maximum suspension width of a tower, 75
   *                 is a good choice for 640 image
   * @param min_suspension_width the minimum suspension width of a tower, 20
   *                 is a good choice for 640 image
   * @return
   */
  public static OneDClusters suspensionBridgeClustering(double[] towers,
    double Rssp, double Tcluster, int max_suspension_width, int min_suspension_width)
  {
      if (towers==null) return null;

      // average value of the bins
      double average = yi.stat.StatUtil.mean(towers);

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

      // calculate the suspension range of the towers
      int[] begin= new int[towers.length];
      int[] end= new int[towers.length];

      for (int i=0; i<towers.length; i++) {
          if (towers[i]>0) {
              int suspensionWidth = (int)Math.ceil(towers[i] / suspensionRatio);
              if (suspensionWidth < min_suspension_width) {
                suspensionWidth = min_suspension_width;
              }
              if (suspensionWidth > max_suspension_width) {
                suspensionWidth = max_suspension_width;
              }
              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] ++; // mark as a begining suspension point of a tower
              end[right] ++; // mark as a ending suspension point of a tower
          }
      }

      // calculate the suspended histogram (a value represents how many towers
      // is suspending the corresponding bin)
      int now = 0;
      int[] sh = new int[towers.length];
      for (int 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 histogram
      Vector clusterCans = new Vector();
      for (int l=0; l<sh.length; l++) {
          if (sh[l] == 0) continue; // sh[l] > 0 => l = left of a cluster
          int r = l; // searching the right
          while (r < sh.length-1 && sh[r+1] > 0) r ++; // sh[r+1] = 0 => r = right

          DataRange range = new DataRange(l, r+1); // range = [l, r] = [l, r+1)
          clusterCans.addElement(range);
          l = r + 1;  // sh[r+1] == 0, skip
      }

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

      Vector filtered_ranges = new Vector();

      // filter the clusters (sum must pass the threshold)
      double threshold = yi.stat.StatUtil.sum(towers) * Tcluster;
      for (int i=0; i<sumRange.length; i++) {
          if (sumRange[i] >= threshold) {
              filtered_ranges.addElement(clusterCans.elementAt(i));
          }
      }

      // return the cluster results
      return new OneDClusters(filtered_ranges, 0, null);
  }

  public static OneDClusters suspensionBridgeClustering(WeightedHistogram hist,
    double Rssp, double Tcluster, int max_suspension_width, int min_suspension_width) {
    OneDClusters clusters = suspensionBridgeClustering(
      hist.getWeightedHistogram(), Rssp, Tcluster, max_suspension_width, min_suspension_width);
    // change the data ranges since the first bin of histogram may not
    // correspond to value 0
    DataRange[] data_ranges = clusters.getClusters();
    for (int i=0; i<data_ranges.length; i++) data_ranges[i].shift(hist.getMinValue());
    return clusters;
  }

  /**
   * To cluster data using histogram breaking.  A break in a histogram is a
   * continuous bins in the histogram with value 0.  A histogram can be broken
   * to different part using those breakings.  However if the the length of a
   * breaking is not long enough, that breaking is not qualified.
   * @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
   */
  public static OneDClusters 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;
      }


      //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
      Vector clusters= new Vector();
      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.add(new DataRange(range[0], range[1]+1));
          }
      }

      return new OneDClusters(clusters, 0, null);
  }


  /**
   * Return the average of lower part of data.  Data will be first sorted and
   * then only the smaller data contribute to the average.  The value is good
   * be used as the threshold of local maxima.  If a local maxima is less than
   * this average, then it is not good enough
   * @param percentage how many data will be used
   * @return the average
   */
  static double averageOfLowerPart(double[] data, double percentage) {

    // copy the data
    double[] copy = (double[])data.clone();

    double sum=0;
    int length = (int)Math.round(data.length * percentage);

    // sort and calcuate the average
    for (int i=0; i<length; i++) {
      for (int j=i+1;j<data.length;j++) {
        if (copy[i]>copy[j]) {
          double temp= copy[i]; copy[i]= copy[j]; copy[j]=temp;
        }
      }
      sum+= copy[i];
    }
    double average= sum/length;
    return average;
  }

  /**
   *
   */
  public static void main(String[] args) {
    double[] data = new double[] {
      1,3,5,10,20,13,12,11,10,10,9,9,12,15,19,20,20,18,16,13,10
    };
    OneDClusters results = OneDHistBasedCluster.localMaximaGuidedHistClustering(data, 0.3, 0, true);
    yi.util.io.OutputResult.print(System.out, results.getClusterMap());
  }

}