package yi.stat.em;

import java.io.*;
import java.util.*;

import org.jmat.data.*;
import org.jmat.gui.*;
import yi.iu.imgdb.learn.*;
import yi.iu.processor.color.*;
import yi.iu.processor.resample.*;
import yi.stat.dist.*;
import yi.stat.hist.*;
import yi.stat.kmeans.*;
import yi.util.config.*;
import yi.util.data.*;
import yi.util.io.*;

/**
 * GMMEM implements the EM algorithm for GMM (Gaussian Mixture
 * Model)
 *
 * The algorithm contains two steps: E-step and M-step.
 *
 * See main() for a simple example.
 */

public class GaussianMixtureModel extends MixtureModel implements Serializable{

  /**
   * When using k-means to init the original means, how many steps the k-means
   * algorithm should run is defined by KMEANS_MAX_STEPS
   */
  public int KMEANS_MAX_STEPS = 5;

  protected SingularGaussianModelFixable m_singularFixer = null;


  static final long serialVersionUID = -2535105722436793009L;

  /**
   * This function can't be called directly.  It is here just for inheritance
   */
  public GaussianMixtureModel(){
    // This function can't be called directly.  It is here just for inheritance
  }

  /**
   * Init the Gaussian Mixture Model with the number of kernel (m), the means of
   * the Gaussians (mu), the covariance matrixes (sigma), the weight of each
   * Gaussian (w), and the observed data (data)
   *
   * @param m
   * @param mu
   * @param sigma
   * @param w
   * @param data
   */
  public GaussianMixtureModel(int m, double[][] mu, double[][][] sigma, double[] w, double[][] data) {
    for (int i=0; i<m; i++) {
      GaussianComponent g = (GaussianComponent)createEMComponent();
      g.setMeans(mu[i]);
      g.setCovariance(sigma[i]);
      m_initComponents.put(new Integer(i), g);
      m_weights.put(new Integer(i), new Double(w[i]));
    }
    m_numComponents = m;
    m_data = data;
  }

  public void setSigularFixer(SingularGaussianModelFixable fixer) {
    m_singularFixer = fixer;
    for (Iterator i= m_initComponents.values().iterator(); i.hasNext();) {
      ((GaussianComponent)i.next()).setSigularFixer(fixer);
    }
  }

  /**
   * Create the same model with the same mean and covariance, but ignore all
   * the data.  This is useful to dump the model to disk and to use later.
   * The required space will be much smaller
   * @return
   */
  public GaussianMixtureModel createLeanInstance() {

    // get means, covariances, weights and num of components
    HashMap means = getComponentMeans();
    HashMap covs = getComponentCovariance();
    HashMap weights = getWeights();
    int m = getNumComponents();

    // create the lean instance
    GaussianMixtureModel lean = new GaussianMixtureModel();

    for (Iterator i=means.keySet().iterator(); i.hasNext();) {
      Object key = i.next();
      double[] mean = (double[])means.get(key);
      double[][] cov = (double[][])covs.get(key);
      Double w = (Double)weights.get(key);
      GaussianComponent g = (GaussianComponent)createEMComponent();
      g.setMeans(mean);
      g.setCovariance(cov);
      lean.m_initComponents.put(key, g);
      lean.m_weights.put(key, w);
    }
    lean.m_numComponents = m;
    lean.m_data = null;
    lean.m_rstComponents = lean.m_initComponents;
    return lean;
  }

  /**
   * Init the Gaussian Mixture Model with the number of kernel (m), and the
   * observed data (data).  The means of the Gaussians (mu), the covariance
   * matrixes (sigma) will be estimated by k-means.  The weight of each
   * Gaussian will be set the same in this initial step.
   *
   * The maximums declares the maximum value of
   * each dimension of the data.  This information is useful for K-Means since
   * that will lead to a correct initial means.  If passed in as null, than
   * a maximum of [1, 1, ..., 1] will be assumed.
   * @param m
   * @param data
   * @param minimums
   * @param maximums
   */
  public GaussianMixtureModel(int m, double[][] data, double[] scale) {

    /** @todo k-means or partition? */
    // k-means
    int[] labels;

    for (int k=2; k<=128; k*=2) {
      KMeans k_means = new KMeans();
      //k, KMCluster.RANDOM, int_data, null);
      labels = k_means.clustering(k, data[0].length, data, scale, 10, KMEANS_MAX_STEPS, false);
    }

    /** @todo test whether this more clusters idea is right */
    // Use a bigger K than the desired number of clusters since some clusters
    // gotten from K-Means may be empty.  It is not hurt to generate more
    // clusters in the k-means steps, since small clusters will be ignored and
    // only the first "m" biggest clusters will be used.  (Maybe? it is better
    // than having small clusters.)
    double MORE_INIT_CLUSTER_RATIO = 1.5;
    int k=(int)(m * MORE_INIT_CLUSTER_RATIO); // set k bigger than m

    while (true) {
      KMeans k_means = new KMeans();
      //k, KMCluster.RANDOM, int_data, null);
      labels = k_means.clustering(k, data[0].length, data, scale, 1, KMEANS_MAX_STEPS, true);

      int ncluster = k_means.getClusterNum();
      if (MORE_INIT_CLUSTER_RATIO <= 1) {
        m = ncluster;
        break;
      }
      if (ncluster < m) { // if the cluster number generated is less than the
                          // desired "m", increase k and re-cluster
        System.out.println("k-means cluster num = " + ncluster + ".  will increase");
        k = (int)(k * m / ncluster * MORE_INIT_CLUSTER_RATIO);  // if the generated cluster is less than m, increase k
      } else if (ncluster > m){ // if the generated clusters are more than
                                // desired, use only the biggest ones.
        double[][] means = k_means.means;
        double[][] picked_means = new double[m][];

        // pick up the biggest clusters
        SimpleHistogram hist = new SimpleHistogram(0, ncluster-1);
        hist.putToBin(labels);
        int[] counts = hist.getHistogram();
        for (int i=0; i<m; i++) {
          // find the highest now
          int highest =counts[0];
          int highest_index = 0;
          for (int j=1; j<ncluster; j++) {
            if (counts[j] > highest) {
              highest = counts[j];
              highest_index = j;
            }
          }
          picked_means[i] = means[highest_index];
          // so that it wouldn't be picked again
          counts[highest_index] = 0;
        }

        // do k-means again using the picked means
        KMeans k_means2 = new KMeans();
        labels = k_means2.clustering(data, m, picked_means, 1, false);
        if (ncluster != m) {
          System.out.println("k-means cluster num = " + ncluster +
                             ".  they shold be the same");
        }
        break;
      } else {
        break;
      }
    }


//    // partition
//    int[] labels = Partition.partition(data.length, data[0].length, m, data);

    System.out.println("component num = " + m);

    // assign the component weights

//    // assign a same weight to each component
//    for (int i=0; i<m; i++) {
//      EMCumulator g = createEMCumulator();
//      m_initAccumulators.put(new Integer(i), g);
//      m_weights.put(new Integer(i), new Double(1.0/m));
//    }

    // assign the weigth according to how many data is clustered to it using
    // k-means
    SimpleHistogram hist = new SimpleHistogram(0, m-1);
    hist.putToBin(labels);
    int[] counts = hist.getHistogram();
    for (int i=0; i<m; i++) {
      GaussianComponent g = (GaussianComponent)createEMComponent();
      m_initComponents.put(new Integer(i), g);
      m_weights.put(new Integer(i), new Double(counts[i]/(double)data.length));
      System.out.println("component #" + i + ", weight = " + m_weights.get(new Integer(i)));
    }

    // assign data to its closest cluster (using results from k-means)
    for (int i=0; i<data.length; i++) {
      int closest = labels[i];
      GaussianComponent g = (GaussianComponent)m_initComponents.get(new Integer(closest));
      g.add(data[i], 1);
    }
    m_data = data;
    m_numComponents = m;
  }

  /**
   * Init the Gaussian Mixture Model with the number of kernel (m), and the
   * observed data (data).  The means of the Gaussians (mu), the covariance
   * matrixes (sigma) will be estimated by the initial data classification
   * (specified by init_labels).
   *
   * @param m    number of kernel
   * @param data data to clustering
   * @param init_labels initial cluster label for each data point
   */
  public GaussianMixtureModel(int m, double[][] data, int[] init_labels) {

    // data point count of each labels
    SimpleHistogram hist = new SimpleHistogram(0, m-1);
    hist.putToBin(init_labels);
    int[] counts = hist.getHistogram();

    // weight of each component is how many of data points have been assign to
    // the corresponding label normalized by the total data points.
    for (int i=0; i<m; i++) {
      GaussianComponent g = (GaussianComponent)createEMComponent();
      m_initComponents.put(new Integer(i), g);
      m_weights.put(new Integer(i), new Double(counts[i]/(double)data.length));
    }

    // assign data to the corresponding component (specified by init_labels)
    for (int i=0; i<data.length; i++) {
      int assign = init_labels[i];
      GaussianComponent g = (GaussianComponent)m_initComponents.get(new Integer(assign));
      g.add(data[i], 1);
    }
    m_data = data;
    m_numComponents = m;
  }

  public GaussianComponent getGaussianComponent(Object index) {

    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    return (GaussianComponent)m_rstComponents.get(index);
  }

  /**
   * Return the probabilities of the data from all the components.  (If the
   * data is the mean of a Gaussian component, then it gets prob of 1, else
   * if will get density(d)/density(mean)
   * @param data
   * @param component_weighted
   * @param normalize
   * @return
   */
  public HashMap getProbUnderComponents(double[] data) {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    Iterator iter = m_rstComponents.keySet().iterator();
    HashMap probs = new HashMap();
    double prob_sum = 0;
    while (iter.hasNext()) {
      Object key = iter.next();
      GaussianComponent component = (GaussianComponent) m_rstComponents.get(key);
      try {
        double density = component.getDensity(data);
        double highest_density = component.getMaxProb();
        double prob = density/highest_density;
        probs.put(key, new Double(prob));
      }
      catch (ModelSingularException e) {
        e.printStackTrace();
      }
    }
    return probs;
  }

  /**
   * Return the minimum allowed Prob on each result component.  (just a
   * guess).  The returned results are sorted by the index of the result
   * components
   * @param use_component_weight
   * @return
   * @throws GaussianModelSingularException
   */
  public double[] getRstComponentMinProb(boolean use_component_weight) throws GaussianModelSingularException{
    ArrayList keys = new ArrayList(getRstComponentKeys());
    java.util.Collections.sort(keys);

    double[] prob = new double[getRstNumComponents()];
    int j_index = 0;
    for (Iterator j = keys.iterator(); j.hasNext(); j_index++) {
      Object key = j.next();
      GaussianComponent gc = getGaussianComponent( key);
      double p = gc.getMinProb(getWeight(key));
      if (use_component_weight) {
        p *= getWeight(key);
      }
      prob[j_index] = p;
    }
    return prob;
  }

  /**
   * Return the minimum allowed Prob on each result component.  (just a
   * guess).  The returned results are sorted by the index of the result
   * components
   * The training data will be assigned to a class.  Then the threshold is
   * selected so that a portion of the data will be assigned to class anyway.
   * @param portion
   * @param use_component_weight
   * @return
   * @throws GaussianModelSingularException
   */
  public double[] getRstComponentMinProbToCoverPortionData(double portion, boolean use_component_weight) throws GaussianModelSingularException{
    ArrayList keys = new ArrayList(getRstComponentKeys());
    java.util.Collections.sort(keys);

    int n_component = getRstNumComponents();

    double[][] prob = new double[n_component][m_data.length];
    for (int i=0; i<m_data.length; i++) {
      HashMap scores = getDensitiesUnderComponents(m_data[i], true, false);
      // classify it
      int best_index = -1;
      double best_p = 0;
      int j_index = 0;
      for (Iterator j = keys.iterator(); j.hasNext(); j_index++) {
        Object key = j.next();
        double p = ((Double)scores.get(key)).doubleValue();
        if (p > best_p) {
          best_p = p;
          best_index = j_index;
        }
      }
      // remember the prob on the classified obj
      if (best_p != 0) {
        prob[best_index][i] = best_p;
      }
    }

    double[] min_prob = new double[n_component];
    int j_index = 0;
    for (Iterator j = keys.iterator(); j.hasNext(); j_index++) {
      Object key = j.next();
      int instance_count = 0;
      for (int i=0; i<m_data.length; i++) {
        if (prob[j_index][i] != 0) instance_count ++;
      }
      int covered_count = (int)(instance_count * portion);
      if (covered_count == 0) covered_count = 1; // at least cover 1

      // sort the probs and keep the largest covered_count
      Arrays.sort(prob[j_index]);
      min_prob[j_index] = prob[j_index][m_data.length - covered_count];
      min_prob[j_index] = min_prob[j_index];

      if (!use_component_weight) {
        min_prob[j_index] *= getWeight(key);
      }
      System.out.println("out of " + prob[j_index].length + " instance");
      System.out.println("  " + instance_count + " are classified as " + key);
      System.out.println("  to keep " + portion + ", " +  covered_count + " will be covered");
      System.out.println("  the cut point is " +  (m_data.length - covered_count) + ", p = " + prob[j_index][m_data.length - covered_count]);
      System.out.println("  after log " +  min_prob[j_index]);
      System.out.println("  after use_component_weight " +  min_prob[j_index]);
    }
    return min_prob;
  }


  /**
   * Create the data set for component.  This method can be overrided.  For
   * example, if return "AlignedGaussianComponent"s, aligned Gaussians will be
   * used.
   * @return
   */
  MixtureComponential createEMComponent() {
    GaussianComponent c = new GaussianComponent();
    c.setSigularFixer(m_singularFixer);
    return c;
  }

  /**
   * Return the covariance matrix of each component
   * @return
   */
  public HashMap getComponentCovariance() {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    HashMap covs = new HashMap();
    Iterator iter = m_rstComponents.keySet().iterator();
    while (iter.hasNext()) {
      Object key = iter.next();
      GaussianComponent accumulator =  (GaussianComponent)m_rstComponents.get(key);
      try {
        double[][] cov = accumulator.getCovariance();
        covs.put(key, cov);
      } catch (GaussianModelSingularException e) {
        e.printStackTrace();
      }
    }

    return covs;
  }



  /**
   * Return the frequency function on one dimension (suppose the Gaussian models
   * are aligned Gaussians, or the covariance matrixes only have number on their
   * diagnal).
   * @param min
   * @param max
   * @param slice_num
   * @return
   */
  public double[] getOneDFrequencyFunction(
      int dim_index, double min, double max, int slice_num) throws GaussianModelSingularException
  {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    double[] ff_total = new double[slice_num];
    Iterator iter = m_rstComponents.keySet().iterator();
    while (iter.hasNext()) {
      Object key = iter.next();
      GaussianComponent accumulator =  (GaussianComponent)m_rstComponents.get(key);
      double[] mean = accumulator.getMeans();
      double[][] covariance = accumulator.getCovariance();
      double mu = mean[dim_index];
      double sigma = covariance[dim_index][dim_index];
      double[] ff = NormalDistribution.getFrequencyFunction(new NormalDistribution(mu, sigma), min, max, slice_num);

      double w = ((Double)m_weights.get(key)).doubleValue();
      for (int i=0; i<ff.length; i++) {
        ff_total[i] += w * ff[i];
      }
    }
    return ff_total;
  }

  private void drawComponentHistogram(ComponentHashMap clusters, int dim) {
    ArrayList cluster_names = new ArrayList(clusters.keySet());
    Collections.sort(cluster_names);
    double[][] hists = new double[cluster_names.size()][];
    Plot2DPanel pp = new Plot2DPanel(); // plot panel
    int i_index = 0;
    for (Iterator i=cluster_names.iterator(); i.hasNext(); i_index++) {
      Object key = i.next();
      GaussianComponent accumulator =  (GaussianComponent)clusters.get(key);
      // get the histogram
      hists[i_index] = accumulator.getHistogram(dim, 0, 1, 100);
      // change to 2D
      double[][] line = HistUtil.histToXYPoints(hists[i_index], 0, 1);
      Matrix matrix = new Matrix(line);
      matrix.toPlot2DPanel(pp, key.toString(), "LINE");
    }
    new FrameView("hist", pp);
  }

  /**
   * Construct a ".m" file format string.  So that the data and the model can
   * be read from Matlab
   * @param data_var variable name for the data
   * @param mu_var variable name for means
   * @param cov_var variable name for covariance matrix
   * @return
   */
  public String toMatlab(String data_var, String mu_var, String cov_var) {
    StringWriter str = new StringWriter();
    PrintWriter out = new PrintWriter(str);

    // y=[1 1 1 2; 2 2 2 1; 1 2 1 2]; one data per column
    for (int i = 0; i < m_data.length; i++) {
      OutputResult.printMatlabMatrix (out, data_var, i+1, m_data[i]);
    }

    // k=2; number of kernel
    //estmu(:,1)=[1 ;1 ;1 ];
    //estmu(:,2)=[2 ;2 ;2 ];
    //estcov(:,:,1) = [0.1 0; 0 0.1];
    //estcov(:,:,2) = [0.1 0; 0 0.1];
    HashMap means = getComponentMeans();
    HashMap covs = getComponentCovariance();
    int i_index = 0; // matrix index starts from 1
    for (Iterator i=means.keySet().iterator(); i.hasNext(); i_index++) {
      Integer key = (Integer)i.next();
      double[] mean = (double[])means.get(key);
      double[][] cov = (double[][])covs.get(key);
      OutputResult.printMatlabMatrix (out, mu_var, i_index+1, mean);
      OutputResult.printMatlabMatrix (out, cov_var, i_index+1, cov);

    }
    return str.toString();
  }


  /** @todo consider where to put  calculateDistributions*/
  /**
   * add the distribution (each one is 1D) of current model to the passed in
   * linked list
   */
  void calculateDistributions(
      GaussianMixtureModel gmm, LinkedList[] ffs,
      int min, int max, int slice_num) throws GaussianModelSingularException{
    for (int i=0; i<ffs.length; i++) {
      double[] ff = gmm.getOneDFrequencyFunction(i, min, max, slice_num);
      if (ffs[i] == null) ffs[i] = new LinkedList();
      ffs[i].add(ff);
    }
  }

  /**
   * Show the
   * @param ffs
   * @return
   */
  public Plot2DPanel showDataModelOnPanel(LinkedList ffs, int i, double min, double max, int DISTRIBUTION_DRAW_SLICE_NUM) {
       // the draw panel
      Plot2DPanel pp = new Plot2DPanel();
      // the data
      double[] oned_data = new double[m_data.length];
      for (int j = 0; j < m_data.length; j++)
        oned_data[j] = m_data[j][i];
//      // output the result
//      OutputResult.outputResult("data", oned_data, 1, oned_data.length);

      /** @todo debug max */
      System.out.println("max = " + max);

      // construct the histogram for data
      int[] hist = HistUtil.histogram(oned_data, min, max,
                                      DISTRIBUTION_DRAW_SLICE_NUM);
      // convert the histogram to double type
      double[] hist_double = Conversion.toDouble(hist);
      // convert to (x,y) style
      double[][] hist_matrix = HistUtil.histToXYPoints(hist_double, min, max);
      Matrix oned_matrix = new Matrix(hist_matrix);
      oned_matrix.toPlot2DPanel(pp, "data", "STAIRCASE");

      // the models
      int m_index = 1;
      for (Iterator j = ffs.iterator(); j.hasNext(); m_index++) {
        double[] dist = (double[]) j.next();
        //        OutputResult.outputResult("model_D" + m_index, dist, 1, dist.length);

        // construct the histogram for model
        double[] dist_norm = new double[dist.length];
        for (int kk = 0; kk < dist.length; kk++) {
          dist_norm[kk] = dist[kk] * m_data.length * (max-min) / DISTRIBUTION_DRAW_SLICE_NUM;
        }
        double[][] dist_m = HistUtil.histToXYPoints(dist_norm, min, max);

        Matrix matrix = new Matrix(dist_m);
        matrix.toPlot2DPanel(pp, Integer.toString(m_index), "LINE");
      }

    return pp;

  }


  public static class ModelData {

    public final static String USAGE = "Usage: [-d img_data_center_file] [-i img_list_file] [-f feature_class_name] [-m M] [-r rst_file] [-c cut_ratio]";

    public static void main(String[] args) {

      try {
        // init configuration
        Configuration.init(args);

        // get the datacenter
        String data_center_name = ReadArgu.getOption(
            'd', args,
            "groundtruth.ser"
            );

        ImgDataCenter data_center = ImgDataCenter.createFromFile(
            data_center_name);

        // get the feature data class names
        String feature_name = ReadArgu.getOption(
            'f', args,
            "yi.iu.processor.color.reg.DS_PolygonRegionColorStat"
            );

        // get the image list file name
        String img_list_file = ReadArgu.getOption('i', args);
        // read image list from a file if it is specified
        LinkedList img_list = null;
        if (img_list_file != null) {
          String[] images = yi.util.io.ReadArgu.readStrings(img_list_file);
          img_list = yi.util.data.Conversion.asLinkedList(images);
        }

        double[][] data = (img_list_file != null) ?
            data_center.groupFeatureVectors(img_list, feature_name) :
            data_center.groupFeatureVectors(feature_name);

        System.out.println("total feature vectors #= " + data.length + "\n");

        String cut_ratio_str = ReadArgu.getOption('c', args, "5");  // 5%
        double cut_ratio = Double.parseDouble(cut_ratio_str);
        LinkedList cut_data = new LinkedList();
        Random r = new Random(0);
        for (int i=0;i<data.length; i++) {
          if (r.nextInt(100) < cut_ratio) cut_data.add(data[i]);
        }
        data = yi.util.data.Conversion.to2DDoubleArray(cut_data);
        // garbage collection
        cut_data.clear();
        cut_data = null;
        System.out.println("total feature vectors cut #= " + data.length + "\n");

        // normalize the data */
        double[] max_data = yi.stat.StatUtil.maximumOnCol(data);
        double[] min_data = yi.stat.StatUtil.minimumOnCol(data);

        // get the negative model first
        // get the datacenter
        String M_str = ReadArgu.getOption(
            'm', args,
            "64"
            );
        int M = Integer.parseInt(M_str);
        System.out.println("finding model... m = " + M);

        GaussianMixtureModel gmm = new AlignedGaussianMixtureModel(
            M, data, (double[])null); // suppose the data is within [0,1)
        gmm.learnTheModel();
        String rst_file = ReadArgu.getOption(
            'r', args,
            "model.ser"
            );
        gmm.writeToFile(rst_file);

        // for SHOW_DISTRIBUTION, save the frequency functions
        // each entry is for one dimension.  Each entry of the linked list
        // is the one for a specified m (as in m mixture models)
        int feature_len = data_center.getFeatureLength(feature_name);
        int slice_num = 100;
        double maximum = yi.stat.StatUtil.maximum(max_data);
        double minimum = yi.stat.StatUtil.minimum(min_data);

        for (int i = 0; i < feature_len; i++) {
          double[] ff = gmm.getOneDFrequencyFunction(i, minimum, maximum, slice_num);
          LinkedList ffs = new LinkedList();
          ffs.add(ff);
          Plot2DPanel pp = gmm.showDataModelOnPanel(ffs, i, minimum,
              maximum,
              slice_num);
          new FrameView("obj_D" + i, pp);
        }

      }
      catch (Exception e) {
        e.printStackTrace();
      }
    }

  }

  public static class Classify {

    public final static String USAGE = "Usage: [-d img_data_center_file] [-m model_file] [-i img_list_file] [-f feature]";

    public static void main(String[] args) {

      // init configuration
      Configuration.init(args);

      try {

        // read in an ObjGMMAdapter object
        String model_file = ReadArgu.getOption(
            'm', args,
            "model0.ser"
            );

        final GaussianMixtureModel gmm = (GaussianMixtureModel)GaussianMixtureModel.createFromFile(model_file);
        int m = gmm.getNumComponents();
        HashMap kernel_means = gmm.getComponentMeans();
        double[][] means = new double[m][];
        /** @todo this only works on HSI */
        int[] mean_color = new int[m + 1]; // make room for background
        mean_color[0] = 0xff000000;
        for (Iterator i=kernel_means.keySet().iterator(); i.hasNext(); ) {
          Integer key = (Integer)i.next();
          int key_i = key.intValue();
          means[key_i] = (double[])kernel_means.get(key);
          mean_color[key_i + 1] = java.awt.Color.HSBtoRGB((float)means[key_i][0],
              (float)means[key_i][1],(float)means[key_i][2]);
        }

        /** @todo color spectrum */
        int[] color_spectrum = ColorUtil.createSpectrum(mean_color, 40, 50);
        OutputResult.outputGIFImage("color_means.gif", color_spectrum, 40*mean_color.length,50);
        int[] color_label = new int[mean_color.length];
        for (int i=0; i<color_label.length; i++) color_label[i] = i;
        int[] label_spectrum = ColorUtil.createSpectrum(color_label, 40, 50);
        int[] label_spectrum_img = yi.util.ImageUtil.makeNewImageFromLabel(
              label_spectrum,
              m+1, // set the maximum cluster number (since we know).
                   //  +1 because the background takes label 0
              true // 0 is treated as background
              );
        OutputResult.outputGIFImage("color_labels.gif", label_spectrum_img, 40*mean_color.length,50);
        //        if (label_spectrum_img != null) System.exit(0);



        // init the datacenter
        String data_center_name = ReadArgu.getOption(
            'd', args,
            "groundtruth.rgn.ser"
            );

        ImgDataCenter data_center = ImgDataCenter.createFromFile(
            data_center_name);

        // get the image list file name
        String img_list_file = ReadArgu.getOption('i', args);

        // image list
        LinkedList img_list = null;
        if (img_list_file != null) {
          String[] images = yi.util.io.ReadArgu.readStrings(img_list_file);
          img_list = yi.util.data.Conversion.asLinkedList(images);
        }

        /** @todo custom the feature name */
        String feature = ReadArgu.getOption('f', args, "yi.iu.processor.texture.reg.DS_PolygonRegionTextureStat");

        HashMap rst_file_suffix_map = new HashMap();
        rst_file_suffix_map.put("yi.iu.processor.texture.reg.DS_PolygonRegionTextureStat", ".s256.text.reg.gif");
        rst_file_suffix_map.put("yi.iu.processor.color.reg.DS_PolygonRegionColorStat", ".s256.color.reg.gif");
        rst_file_suffix_map.put("yi.iu.processor.clc.LineClusterSets", ".CLC_LineStruct.bound.gif");
        rst_file_suffix_map.put("yi.iu.processor.color.reg.ColorLabeledRgnFts", ".s256.color.rgn.gif");

        /** @todo make scale setable */
        // a resampler to resize the image
        NearestResampler resampler = new NearestResampler(128, 128);
        int[][] cluster_dat = new int[img_list.size()][];

        // loop through the image and their features
        int i_index = 0;
        for (Iterator i=img_list.iterator(); i.hasNext(); i_index++) {
          String img = (String)i.next();
          RgnFeatureRepresentationList rfrl =
              (RgnFeatureRepresentationList)data_center.getImgFeatureList(img, feature);
          int[] score_img = rfrl.drawScoreDemo(
              null,
              new FeatureScoringCapable(){
                public double getFeatureScore(double[] feature, Object obj) {
                  return ((Integer)gmm.classify(feature)).intValue();
                }
              },
              1, // shift = 1 because the cluster #0 should be labeled as 1 to
                 // avoid using the same number as background
              1);

          // in the pixels_img, each pixel is represented by a cluster label
          // 0 -> background, others -> cluster_index
          int[] pixels_img = yi.util.ImageUtil.makeNewImageFromLabel(
              score_img,
              m+1, // set the maximum cluster number (since we know).
                   //  +1 because the background takes label 0
              true // 0 is treated as background
              );
          /** @todo make image suffix setable */
          OutputResult.outputGIFImage(img + ".cc.gif", pixels_img,
                                                 rfrl.imgWidth, rfrl.imgHeight); // color_cluster

          int[] color_mean_img = yi.util.ImageUtil.makeNewImageFromLabel(score_img, mean_color);
          yi.util.io.OutputResult.outputGIFImage(img + ".cc.mc.gif", color_mean_img,
                                                 rfrl.imgWidth, rfrl.imgHeight); // color_cluster mean_color


          System.out.println("<P>" +
                             "<IMG SRC=\"" + img + "\" width=256>" +
                             "<IMG SRC=\"" + img + (String)rst_file_suffix_map.get(feature) + "\" width=256>" +
                             "<IMG SRC=\"" + img + ".cc.gif\" width=256>" +
                             "<IMG SRC=\"" + img + ".cc.mc.gif\" width=256>" +
                             "</P>");

          cluster_dat[i_index] = resampler.strecth(score_img, rfrl.imgWidth, rfrl.imgHeight);

        }
        OutputResult.print("cc.dat", cluster_dat, 128*128, cluster_dat.length);

      } catch (Exception e) {
        e.printStackTrace();
      }
    }
  }

  public static class ModelToMatlab {

    public final static String USAGE = "Usage: [-m model_file]";

    public static void main(String[] args) {

      // init configuration
      Configuration.init(args);

      try {

        // read in an ObjGMMAdapter object
        String model_file = ReadArgu.getOption(
            'm', args,
            "model0.ser"
            );

        final GaussianMixtureModel gmm = (GaussianMixtureModel)GaussianMixtureModel.createFromFile(model_file);
        System.out.println(gmm.toMatlab("y", "mu", "cov"));

      } catch (Exception e) {
        e.printStackTrace();
      }
    }
  }

  public static class MatlabToModel {
    public final static String USAGE = "Usage: [-d directory_of_matlab_files] [-r rst_file]";

    public static void main(String[] args) {

      // init configuration
      Configuration.init(args);

      try {

        // read in an ObjGMMAdapter object
        String directory = ReadArgu.getOption(
            'd', args,
            "mixturecode"
            );

        // read mu
        double[][] mu = ReadArgu.loadMatlab(directory + "/" + "bestmu.dat");
        mu = yi.util.data.Conversion.swapXandY(mu, mu[0].length, mu.length);
        double[][] pp = ReadArgu.loadMatlab(directory + "/" + "bestpp.dat");

        int m = pp[0].length;
        double[][][] cov = new double[m][][];
        for (int i=0; i<m; i++) {
          String cov_name = directory + "/" + "bestcov" + (i+1) + ".dat";
          cov[i] = ReadArgu.loadMatlab(cov_name);
        }

        GaussianMixtureModel gmm = new GaussianMixtureModel(m, mu, cov, pp[0], null);
        gmm.m_rstComponents = gmm.m_initComponents;

        String rst_file = ReadArgu.getOption(
            'r', args,
            "model.ser"
            );
        gmm.writeToFile(rst_file);

      } catch (Exception e) {
        e.printStackTrace();
      }
    }
  }

  public static class Peek {
    public final static String USAGE = "Usage: -m model_file";

    public static void main(String[] args) {

    // init configuration
    Configuration.init(args);

    try {

      // read in an ObjGMMAdapter object
      String model_file = ReadArgu.getOption(
          'm', args,
          "model0.ser"
          );

      GaussianMixtureModel gmm = (GaussianMixtureModel)GaussianMixtureModel.createFromFile(
          model_file);
      System.out.println(gmm);

    }
    catch (Exception e) {
      e.printStackTrace();
    }
  }

  }

  public static void main(String[] args) {

    String str_sun = "Sun";
    String str_tree = "Tree";
    String str_water = "Water";

    String[] objects = new String[]{str_sun, str_tree, str_water};

    int m = 4; // number of kernels
    int n = 2; // dimension of data vector
    int k = 200; // number of data

    // means for each object
    double[][] means = {
        {50, 50},
        {100, 50},
        {40, 120},
        {90, 70}
    };

    // variances
    double[][][] variances = {
        {{20,  0},
         { 0, 30}},
        {{50,  0},
         { 0, 30}},
         {{50,  0},
          { 0, 30}},
         {{25,  0},
          { 0, 30}},
    };

    // probability for each object
    double[] obj_probs = {0.2, 0.2, 0.3, 0.3};


//    int m = 3; // number of kernels
//    int n = 3; // dimension of data vector
//    int k = 300; // number of data
//
//    // means for each object
//    double[][] means = {
//        {180, 50, 50},  // sun
//        {20, 150, 50},  // tree
//        {20, 20, 200}   // water
//    };
//
//    // variances
//    double[][][] variances = {
//        {{20,  0,  0},   // sune
//         { 0, 30,  0},
//         { 0,  0, 30}},
//        {{15,  0,  0},   // tree
//         { 0, 50,  0},
//         { 0,  0, 30}},
//        {{15,  0,  0},   // water
//         { 0, 15,  0},
//         { 0,  0, 50}},
//    };
//
//    // probability for each object
//    double[] obj_probs = {0.2, 0.5, 0.3};


//    int m = 3; // number of kernels
//    int n = 1; // dimension of data vector
//    int k = 300; // number of data
//
//    // means for each object
//    double[][] means = {
//        {50},
//        {80},
//        {100}
//    };
//
//    // variances
//    double[][][] variances = {
//        {{40}},
//        {{20}},
//        {{50}}
//    };
//
//    // probability for each object
//    double[] obj_probs = {
//        0.4, 0.3, 0.3};


//    int m = 1; // number of kernels
//    int n = 1; // dimension of data vector
//    int k = 300; // number of data
//
//    // means for each object
//    double[][] means = {
//        {50}
//    };
//
//    // variances
//    double[][][] variances = {
//        {{40}}
//    };
//
//    // probability for each object
//    double[] obj_probs = {
//        1};


    double[] cum = new double[m];
    cum[0] = obj_probs[0];
    for (int i=1; i<m; i++) {
      cum[i] = cum[i-1] + obj_probs[i];
    }


    double data[][] = new double[k][n];
    Random wrandom = new Random();
    Random grandom = new Random();
    for (int i=0; i<k; i++) {
      double r = wrandom.nextDouble();
      int gaussian_index = 0;
      for (gaussian_index=0; gaussian_index<m; gaussian_index++) {
        if (r<=cum[gaussian_index]) break;
      }

      for (int j=0; j<n; j++) {
        data[i][j] = grandom.nextGaussian() * Math.sqrt(variances[gaussian_index][j][j]) + means[gaussian_index][j];
      }
    }

//    GaussianMixtureModel gmm = new AlignedGaussianMixtureModel(m, data, maximums);
    GaussianMixtureModel gmm = new GaussianMixtureModel(m, data, (double[])null);

    gmm.learnTheModel();
    System.out.println(gmm.toString());

    // get the distribution of current model
    int feature_len = n;
    int DISTRIBUTION_DRAW_SLICE_NUM = 100;
    int min =0, max=256;
    LinkedList[] ffs = new LinkedList[feature_len];
    try {
      gmm.calculateDistributions(gmm, ffs, min, max, DISTRIBUTION_DRAW_SLICE_NUM);
    } catch (Exception e) {
      e.printStackTrace();
    }

    // draw the distributions fro each dimension of the feature
    for (int i=0; i<ffs.length; i++) {
      Plot2DPanel pp = gmm.showDataModelOnPanel(ffs[i], i, min, max, DISTRIBUTION_DRAW_SLICE_NUM);
      new FrameView("ROC", pp);
    }
  }
}
