package yi.stat.em;

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

import yi.ds.*;
import yi.iu.processor.color.*;
import yi.util.io.*;
import yi.util.io.output.*;
import java.text.NumberFormat;
import yi.stat.*;

/**
 * Mixture model supported by EM algorithm
 */

public abstract class MixtureModel implements Serializable{

  /**
   * Store the init components for the objects.  The key is the name of the
   * object and the value is the corresponding component (DoubleArraycomponent)
   */
  protected ComponentHashMap m_initComponents = new ComponentHashMap();

  /**
   * Store the final components for the objects after iteration step.  The key
   * is the name of the object and the value is the corresponding component
   * (DoubleArraycomponent)
   */
  protected ComponentHashMap m_rstComponents = null;

  /**
   * The weight of each component
   */
  protected HashMap m_weights = new HashMap();

//  public void setEqualWeights() {
//  HashMap w = m_weights;
//  for (Iterator i=w.keySet().iterator(); i.hasNext(); ) {
//    w.put(i.next(), new Double(1.0/w.size()));
//  }
//}
//

  /**
   * The number of component used
   */
  protected int m_numComponents;


  /**
   * The observed data.  The first dimension is the number of data; the second
   * is the dimension of each data
   */
  transient protected double[][] m_data = null;

  /**
   * The break threshold for E-M method.  If
   * (current_likelihood - last_likelihood) / last_likelihood < BREAK_THRESHOLD,
   * the EM iterations will stop.
   */
  public double BREAK_THRESHOLD = 1e-7;

  public int MAX_ITERATION_STEP = 500;

  public boolean SPLITTING_AND_VANISHING = false;
  public double MCVR = 3.0;
  public double MCSR = 3.0;
  public int m_splitCount = 0;

  public boolean SHOW_PROGRESS = false;
  public boolean DRAW_COMPONET = false;

  static final long serialVersionUID =  -6426713474372283031L;

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

  /**
   * Create mixture model component.  This method can be overrided.  For
   * example, if return "AlignedGaussianComponent"s, aligned Gaussians will be
   * used.
   * @return
   */
  abstract MixtureComponential createEMComponent();

  /**
   * Return the number of components.  After learning, some components may be
   * removed.  This returned number may be the initial number of kernel, not the
   * number after training/learning/modeling.
   * @return
   */
  public int getNumComponents() {
    return m_numComponents;
  }

  /**
   * Return of the number of resultant components (after training)
   * @return
   */
  public int getRstNumComponents() {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    return m_rstComponents.size();
  }

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

    return m_rstComponents;
  }

  /**
   * Return the list of result components
   * @return
   */
  public Collection getRstComponentKeys() {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    return m_rstComponents.keySet();
  }


  /**
   * Represent the component means in color.  Those deleted component will
   * be return as white (0xffffffff)
   * @param color_space can be "lab", "hsi"
   * @return
   */
  public int[] getRstMeansInColor(String color_space) {
    HashMap means = getComponentMeans();
    ArrayList sorted_key = new ArrayList(m_initComponents.keySet());
    Collections.sort(sorted_key);
    int[] colors = new int[sorted_key.size()];

    boolean use_int = (sorted_key.get(0) instanceof Integer);
    if (use_int) {
      int[] keys = yi.util.data.Conversion.toIntArray(sorted_key);
      int max_key = StatUtil.maximum(keys);
      colors = new int[max_key + 1];
      Arrays.fill(colors, 0xffffffff);
      for (int i = 0; i<keys.length; i++) {
        Object key = new Integer(keys[i]);
        double[] data = (double[]) means.get(key);
        if (data == null) {
          colors[keys[i]] = 0xffffffff;
          continue;
        }
        if (color_space.compareToIgnoreCase("lab") == 0) {
          double[] cie = new double[] {
              data[0], data[1], data[2]
          };
          colors[keys[i]] = yi.iu.processor.color.CIELab.lab2rgb(cie);
        }
        else if (color_space.compareToIgnoreCase("hsi") == 0) {
          colors[keys[i]] = java.awt.Color.HSBtoRGB( (float) data[0],
                                                    (float) data[1],
                                                    (float) data[2]);
        } else {
          return null;
        }
      }
      return colors;
    }

    int i_index = 0;
    for (Iterator i = sorted_key.iterator(); i.hasNext(); i_index++) {
      Object key = i.next();
      double[] data = (double[]) means.get(key);
      if (data == null) {
        colors[i_index] = 0xffffffff;
        continue;
      }
      if (color_space.compareToIgnoreCase("lab") == 0) {
        double[] cie = new double[] {
            data[0], data[1], data[2]
        };
        colors[i_index] = yi.iu.processor.color.CIELab.lab2rgb(cie);
      }
      else if (color_space.compareToIgnoreCase("hsi") == 0) {
        colors[i_index] = java.awt.Color.HSBtoRGB( (float) data[0],
                                                  (float) data[1],
                                                  (float) data[2]);
      } else {
        return null;
      }
    }
    return colors;
  }

  /**
   * Return component weights (including those removed)
   * @return
   */
  public double[] getRstComponentWeightsInSortedArray() {
    int n = getNumComponents();
    ArrayList sorted_key = new ArrayList(m_weights.keySet());
    Collections.sort(sorted_key);
    double[] prob_array = new double[sorted_key.size()];

    HashMap probs = getWeights();
    int i_index = 0;
    for (Iterator i = sorted_key.iterator(); i.hasNext(); i_index++) {
      Integer key = (Integer) i.next();
      Double data = (Double) probs.get(key);
      if (data == null) {
        prob_array[i_index] = 0;
        continue;
      }
      else  {
        prob_array[i_index] = data.doubleValue();
      }
    }
    return prob_array;
  }

  /**
   * This method triggers the iteration steps, which really do E-M.
   */
  public void learnTheModel() {
    System.out.println("Gaussian Mixture Model Learning ... " +
                       this.getClass().getName() +
                       ", M = " + this.m_numComponents +
                       ", break threshold = " + this.BREAK_THRESHOLD +
                       ", aligned = " + (this.createEMComponent() instanceof AlignedGaussianComponent) +
                       ", split&vanish = " + SPLITTING_AND_VANISHING);
   if (SPLITTING_AND_VANISHING) {
     System.out.println("vanish ratio = " + MCVR + ", split ratio = " + this.MCSR);
   }
   // go iteration steps
    iterationSteps();
  }

  /**
   * After learning, some original Components may be removed.  So the index
   * in the resultant model is not continous anymore.  This method will make
   * the index again continous.
   */
  public void removeEmptyClusterIndex() {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    System.out.println("removeEmptyClusterIndex");

    ComponentHashMap new_components = new ComponentHashMap();
    HashMap new_weight = new HashMap();

    ArrayList keys = new ArrayList(m_rstComponents.keySet());
    java.util.Collections.sort(keys);
    int j_index = 0;
    for (Iterator j = keys.iterator(); j.hasNext(); j_index++) {
      Object key = j.next();
      Object value = m_rstComponents.get(key);
      new_components.put(new Integer(j_index), value);
      value = m_weights.get(key);
      new_weight.put(new Integer(j_index), value);
      System.out.println(((Integer)key).toString() + " -> " + j_index);
    }
    m_rstComponents.clear();
    m_rstComponents = new_components;
    m_weights.clear();
    m_weights = new_weight;
    System.out.println("number component = " + m_rstComponents.size());
  }

  public void removeComponent(int index) {
    m_rstComponents.remove(new Integer(index));
    m_weights.remove(new Integer(index));
  }

  /**
   * Return the value on that data from density function, p(data|gmm).
   * @param data
   * @return
   */
  public double getDensity(double[] data)  throws ModelSingularException{

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

    // use the calculated model to get the density
    return getDensity(m_rstComponents, m_weights, data);
  }

  /**
   * Return the density of the given data under the specified component
   * @param component
   * @param data
   * @return
   * @throws EmModelSingularException
   */
  public double getDensity(Object obj, double[] data) throws ModelSingularException{
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }

    MixtureComponential component = (MixtureComponential) m_rstComponents.get(obj);
    return component.getDensity(data);
  }

  /**
   * Get the weight of the specified component.  The returned weight is not
   * necessarily the weight after learning.  So it can be the initial weight
   * @param obj
   * @return
   */
  public double getWeight(Object component) {
    return ((Double)m_weights.get(component)).doubleValue();
  }


  /**
   * Return the learned weights of the components
   * @return
   */
  public HashMap getWeights() {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }
    return m_weights;
  }

  /**
   * Return p(data|model)
   * @param data
   * @return
   */
  public HashMap getDensitiesUnderComponents(double[] data) {
    return getDensitiesUnderComponents(data, true, true);
  }

  public HashMap getDensitiesUnderComponents(double[] data, boolean component_weighted, boolean normalize) {
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }
    return getDensitiesUnderComponents(m_rstComponents, component_weighted? m_weights : null, data, normalize);
  }

  /**
   * Return the mean of each component in the given HashMap.  It is not
   * necessary that the mean are those after learning.
   * @param clusters
   * @return
   */
  protected HashMap getComponentMeans(ComponentHashMap clusters) {
    HashMap means = new HashMap();
    Iterator iter = clusters.keySet().iterator();
    while (iter.hasNext()) {
      Object key = iter.next();
      MixtureComponential component =  (MixtureComponential)clusters.get(key);
      try {
        double[] mean = component.getMeans();
        means.put(key, mean);
      } catch (ModelSingularException e) {
        e.printStackTrace();
      }
    }

    return means;
  }


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

    return getComponentMeans(m_rstComponents);
  }

  public MixtureComponential getComponent(Object index) {

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

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


  public double getParameterLogLikelihood(double[][] data) throws ModelSingularException{
    // calculate the model if it is not calculated yet
    if (m_rstComponents == null) {
      learnTheModel();
    }
    return getParameterLogLikelihood(m_rstComponents, m_weights, data);
  }


  /**
   * Save the mean of each cluster in each iteration
   */
  LinkedListHashMap mean_hist = new LinkedListHashMap();
  protected void addToMeanHist(LinkedListHashMap mean_hist, HashMap means) {
    for (Iterator i=means.keySet().iterator(); i.hasNext(); ) {
      Object key = i.next();
      mean_hist.addToList(key, means.get(key));
    }
  }
  public boolean OUTPUT_MEANS = false;
  private void showMeanHist(LinkedListHashMap mean_hist, String img_name) {
    // cluster name (sorted)
    ArrayList cluster_names = new ArrayList(mean_hist.keySet());
    Collections.sort(cluster_names);

    // for each cluster
    int[][] demo_colors = new int[cluster_names.size()][];
    int i_index = 0;
    int max_hist_len = 0;
    for (Iterator i = cluster_names.iterator(); i.hasNext(); i_index++) {
      Object cluster_name = i.next();
      if (OUTPUT_MEANS) System.out.println(cluster_name);

      LinkedList hist = mean_hist.getList(cluster_name);
      demo_colors[i_index] = new int[hist.size()]; // use the first 3 for color demo
      if (hist.size() > max_hist_len) max_hist_len = hist.size();
      int j_index = 0;
      for (Iterator j= hist.iterator(); j.hasNext(); j_index++) {
        double[] mean = (double[]) j.next();
        if (OUTPUT_MEANS) OutputResult.print(System.out, mean);
        if (mean.length >= 3) {
          demo_colors[i_index][j_index] = java.awt.Color.HSBtoRGB( (float) mean[0],
              (float) mean[1], (float) mean[2]);
        }
      }
    }
    if (img_name != null) {
      int[] demo_img = ColorUtil.createSpectrum(demo_colors, max_hist_len, 10, 20, 0xff000000);
      OutputResult.outputColorImage(img_name, demo_img, 10 * max_hist_len, 20 * cluster_names.size());
    }

  }

  /**
   * maximum interation step for mixture of gmm
   */
  int MAXIMUM_GMMMIXTURE_ITER = 200;
  /**
   * remember the gaussian component means (in gmm mixture model) for iterations
   * the key is the object, the value is double[][], first dimension is iteration
   * second dimension is the mean colors for each Gaussian component
   */
  transient HashMap m_gmmmm_gc_color_means = new HashMap();
  /**
   * remember the gaussian component weights (in gmm mixture model) for iterations
   * the key is the object, the value is double[][], first dimension is iteration
   * second dimension is the weights for each Gaussian component
   */
  transient HashMap m_gmmmm_gc_weights = new HashMap();

  /**
   * The iteration step is made of E-step and M-step
   */
  protected void iterationSteps() {
    try{

      // iteration statistical info
      LinkedList likelihoods = new LinkedList();

      learnComponentModels(m_initComponents);
      ComponentHashMap last_components = (ComponentHashMap) m_initComponents.
          clone();
      if (this instanceof GMMMixtureModel) {
       int a=0;
      }
      double last_likelihood = getParameterLogLikelihood(last_components,
          m_weights, m_data);

      double best_likelihood = last_likelihood;
      ComponentHashMap best_components = last_components;
      if (this instanceof GMMMixtureModel) {
        System.out.println("iter 0, likelihood = " + last_likelihood);
        for (Iterator o = m_initComponents.keySet().iterator(); o.hasNext(); ) {
          Object obj = o.next();
          MixtureModel mm = (MixtureModel)last_components.get(obj);
          int[][] gmmmm_gc_color_means = new int[MAXIMUM_GMMMIXTURE_ITER][];
          gmmmm_gc_color_means[0] = mm.getRstMeansInColor("hsi");
          m_gmmmm_gc_color_means.put(obj, gmmmm_gc_color_means);
          double[][] gmmmm_gc_weights = new double[MAXIMUM_GMMMIXTURE_ITER][];
          gmmmm_gc_weights[0] = mm.getRstComponentWeightsInSortedArray();
          m_gmmmm_gc_weights.put(obj, gmmmm_gc_weights);
        }
      }
      else {
        likelihoods.addLast(new Double(last_likelihood));
      }
      double like = 0;

      // remember mean history
      addToMeanHist(mean_hist, getComponentMeans(m_initComponents));
      // see the histogram of each component (Gaussian)
      if (DRAW_COMPONET) {
//      drawComponentHistogram(m_initComponents, 0);
      }

      /** @todo timing */
      // timing
      long time_e = 0, time_m = 0, time_i = 0, time_l = 0;

      // init a empty component set
      ComponentHashMap new_components = createNewComponents(last_components);
      int iteration_step = 1;
      do {

        long time1 = System.currentTimeMillis();

        // E-step
        for (int i = 0; i < m_data.length; i++) {
          double[] xk = m_data[i];
          // p(wi | xk)
          HashMap probs = getDensitiesUnderComponents(last_components,
              m_weights, xk);

          Iterator prob_iter = probs.keySet().iterator();
          while (prob_iter.hasNext()) {
            Object key = prob_iter.next();
            double prob = ( (Double) probs.get(key)).doubleValue(); // p(wi|xk)
            if (Double.isNaN(prob)) {
              System.err.println("data xk got NaN prob for component \"" + key +
                                 "\"");
              getDensitiesUnderComponents(last_components,
                                          m_weights, xk);
            }
            new_components.add(key, xk, prob);
          }
        }

        long time2 = System.currentTimeMillis();

        // M-step
        // recaculate the obj_probs (weights for each component)
        recalculateWeights(new_components);

        long time3 = System.currentTimeMillis();

        // check the integrity of the new models
        learnComponentModels(new_components);

        // remember the mean
        addToMeanHist(mean_hist, getComponentMeans(new_components));

        long time4 = System.currentTimeMillis();

        // check the likelihood of the models
        // break if it it the time
        like = getParameterLogLikelihood(new_components, m_weights, m_data);

        System.err.println("iter = " + iteration_step + ", likelihood = " + like);
//        for (Iterator j = new_components.keySet().iterator(); j.hasNext();) {
//          Object obj = j.next();
//          System.out.println("component = " + obj);
//          System.out.println((MixtureComponential)new_components.get(obj));
//        }

        if (like < -1e10 || Double.isNaN( (like))) {
          System.err.println("likelihood is undesirable: " + like);
//          {
//            double like2 = 0;
//            for (int i = 0; i < m_data.length; i++) {
//              double[] xk = m_data[i];
//              double d = getDensity(new_components, m_weights, xk);
//              OutputResult.print(System.err, xk);
//              System.err.println("density = " + d);
//              like2 += Math.log(d);
//              System.err.println("like = " + like2);
//              if (like2 < -1e6 || Double.isNaN( (like2))) {
//                getDensity(new_components, m_weights, xk);
//              }
//            }
//          }
//          break;
        }

        if (this instanceof GMMMixtureModel) {

          if (like > best_likelihood) {
            best_likelihood = like;
            best_components = new_components;
          }

          getParameterLogLikelihood(new_components, m_weights, m_data);
          System.out.println("outer iter = " + iteration_step + ", likelihood = " +
                             like);
          for (Iterator o = m_initComponents.keySet().iterator(); o.hasNext(); ) {
            Object obj = o.next();
            MixtureModel mm = (MixtureModel)new_components.get(obj);
            int[][] gmmmm_gc_color_means = (int[][]) m_gmmmm_gc_color_means.get(
                obj);
            if (iteration_step >= gmmmm_gc_color_means.length) break;
            gmmmm_gc_color_means[iteration_step] = mm.getRstMeansInColor("hsi");
            double[][] gmmmm_gc_weights = (double[][]) m_gmmmm_gc_weights.get(
                obj);
            gmmmm_gc_weights[iteration_step] =
                mm.getRstComponentWeightsInSortedArray();
          }
          if (iteration_step >= MAXIMUM_GMMMIXTURE_ITER && like >= best_likelihood ||
              iteration_step >= MAXIMUM_GMMMIXTURE_ITER * 2) {
//            // USE THE ONE WITH BEST LIKELIHOOD
//            new_components = best_components;
//            System.out.println("use the one with best likelihood of " + best_likelihood);

            for (Iterator o = m_initComponents.keySet().iterator(); o.hasNext(); ) {
              Object obj = o.next();
              int[][] gmmmm_gc_color_means = (int[][]) m_gmmmm_gc_color_means.
                  get(obj);
              double[][] gmmmm_gc_weights = (double[][]) m_gmmmm_gc_weights.get(
                  obj);
//              OutputResult.print(obj.toString() + ".gc.color.means",
//                                 gmmmm_gc_color_means);
              int[] color_means_demo = ColorUtil.createSpectrum(
                  gmmmm_gc_color_means,
                  gmmmm_gc_color_means[0].length, 50, 50, 0xffffffff);
              OutputResult.outputGIFImage(obj.toString() + ".gc.color.means.gif",
                                          color_means_demo,
                                          50 * gmmmm_gc_color_means[0].length,
                                          50 * gmmmm_gc_color_means.length);
//              OutputResult.print(obj.toString() + ".gc.weights",
//                                 gmmmm_gc_weights);
            }
            System.out.println("GMM mixture iteration maximum set to " +
                               MAXIMUM_GMMMIXTURE_ITER + ", so jump out");
            break;
          }
        }
        else {
          likelihoods.addLast(new Double(like));
        }

        if (Math.abs(like - last_likelihood) <
            BREAK_THRESHOLD * Math.abs(last_likelihood)) {

          if (SPLITTING_AND_VANISHING) {
            Random r = new Random();
            boolean changed = false;
            // use a splitting and vanishing algorithm during training to learn
            // the number of mixture components.  Two thresholds are defined,
            // a mixture-coefficient vanishing ratio (MCVR), and a
            // mixture-coefficient splitting ration (MCSR).  Under a M-component
            // mixture, with component probabilities p_m, if p_m < 1/(M*MCVR),
            // then the m-th component will vanish.   If p_m>MCSR/M, that
            // component will split.

            int current_comp_num = new_components.size();
            // first, vanishing
            // second, splitting
            for (Iterator c=new_components.keySet().iterator(); c.hasNext();) {
              Object comp_key = c.next();
              // get the weight of the component
              double weight = ((Double)m_weights.get(comp_key)).doubleValue();
              if (weight < 1 / (current_comp_num * MCVR)) {
                System.out.println("vanish component " + comp_key + " , weight = " + weight);
                c.remove();
                m_weights.remove(comp_key);
                changed = true;
                continue;
              }
              if (weight > MCSR / current_comp_num) {
                int new_comp_index = getNumComponents() + (m_splitCount++);
                MixtureComponential comp = (MixtureComponential)new_components.get(comp_key);
                MixtureComponential new_comp = createEMComponent();
                if (new_comp instanceof GaussianComponent) {
                  GaussianComponent gc = (GaussianComponent)new_comp;
                  double[] new_means = (double[])comp.getMeans();
                  double[][] cov = (double[][])((GaussianComponent)comp).getCovariance().clone();
                  for (int d=0; d<new_means.length; d++) {
                    new_means[d] += r.nextGaussian() * Math.sqrt(cov[d][d]);
                  }
                  gc.setMeans(new_means);
                  gc.setCovariance(cov);
                  // add the new component
                  System.out.println("split component " + comp_key + " to " + new_comp_index + ", weight = " + weight/2);
                  new_components.put(new Integer(new_comp_index), new_comp);
                  m_weights.put(new Integer(new_comp_index), new Double(weight/2));
                  m_weights.put(comp_key, new Double(weight/2));
                  break;
                } else {
                  System.err.println("splitting and vanishing algorithm only works for Gaussian components");
                  break;
                }
              }
            }
            if (!changed) break;
          } else {
              break;
          }
        }
        last_likelihood = like;

        long time5 = System.currentTimeMillis();

        iteration_step++;
        if (iteration_step > MAX_ITERATION_STEP) {
          System.err.println("iteration steps greater than MAX_ITERATION_STEP " +
                             MAX_ITERATION_STEP);
          break;
        }
        last_components = new_components;
        new_components = new ComponentHashMap();

        if (SHOW_PROGRESS) {
          System.out.println(like);
        }

        time_e += time2 - time1;
        time_m += time3 - time2;
        time_i += time4 - time3;
        time_l += time5 - time4;

      }
      while (true);
//      System.out.println("likelihood for the parameters = " + like);
//
//      System.out.println("timing e=:" + (time_e / (double) iteration_step) +
//                         ", m=" + (time_m / (double) iteration_step) +
//                         ", checkModelsIntegrity = " +
//                         (time_i / (double) iteration_step) +
//                         ", getParameterLikelihood = " +
//                         (time_l / (double) iteration_step));
//
//      // show the mean hist
//      showMeanHist(mean_hist, "mean_hist.ppm");

      // keep the stable means
      m_rstComponents = new_components;

//      if (SHOW_PROGRESS)
//        System.out.println();
      double[] likelihood_iters = yi.util.data.Conversion.toDoubleArray(
          likelihoods);
      OutputResult.print(System.out, likelihood_iters);
    }
    catch (ModelSingularException e) {
      e.printStackTrace();
    }
  }

  /**
   * When training by iterations, a new set of components are created in each
   * iteration.  This is for create it
   * @param last_components
   * @return
   */
  protected ComponentHashMap createNewComponents(ComponentHashMap last_components) {
    return new ComponentHashMap();
  }

  /**
   * Recalculate the weights of each component and update m_weights
   * @param new_components
   */
  protected void recalculateWeights(ComponentHashMap components) {
    Iterator iter = components.keySet().iterator();
    while (iter.hasNext()) {
      Object key = iter.next();
      MixtureComponential component = (MixtureComponential) components.get(key);
      double prob = component.getProbabilitySum() / m_data.length;
      m_weights.put(key, new Double(prob));
    }

  }

  /**
   * Calculate the density given the components and their weights
   * @param components
   * @param weights
   * @param data
   * @return
   */
  protected static double getDensity(ComponentHashMap components, HashMap weights, double[] data) throws ModelSingularException{

    double prob_sum = 0;
    for (Iterator iter = components.keySet().iterator(); iter.hasNext();) {
      Object key = iter.next();
      MixtureComponential component = (MixtureComponential) components.get(key);
      double prob = component.getDensity(data);
      double model_prob = ( (Double) weights.get(key)).doubleValue();
      prob_sum += prob * model_prob;
    }

    return prob_sum;
  }


  /**
   * Return the probabilities of the data under the different models represented
   * by the GMM.  This is the E-Step in EM algorithm
   * @param components
   * @param model_probs
   * @param data
   * @return
   */
  protected HashMap getDensitiesUnderComponents(HashMap components,
                                         HashMap model_probs, double[] data) {
    return getDensitiesUnderComponents(components, model_probs, data, true);
  }

  /**
   * Return the probabilities of the data under the different models represented
   * by the GMM.  This is the E-Step in EM algorithm
   * @param components
   * @param model_probs
   * @param data
   * @param normalize whether normalize the probabilities
   * @return
   */
  protected HashMap getDensitiesUnderComponents(HashMap components,
                                         HashMap model_probs, double[] data, boolean normalize) {
    Iterator iter = components.keySet().iterator();
    HashMap probs = new HashMap();
    double prob_sum = 0;
    while (iter.hasNext()) {
      Object key = iter.next();
      MixtureComponential component = (MixtureComponential) components.get(key);
      try {
        double prob = component.getDensity(data);
        double model_prob = (model_probs == null) ? 1 : ( (Double) model_probs.get(key)).doubleValue();
        prob *= model_prob;
        if (prob > 0) {
          probs.put(key, new Double(prob));
          prob_sum += prob;
        } else {
          probs.put(key, new Double(0));
        }
      }
      catch (ModelSingularException e) {
        e.printStackTrace();
      }
    }

    // it is possible that a data point gets probabilities of zeros for all the
    // components.  In that case, assgin it to the components with equal prob.
    if (prob_sum == 0 && normalize) {
      System.err.println("data got 0 probs from all the components, assume equal probs");
      iter = probs.keySet().iterator();
      while (iter.hasNext()) {
        Object key = iter.next();
        probs.put(key, new Double(1.0 / probs.size()));
      }
      return probs;
    }

    // normalize
    if (normalize) {
      iter = probs.keySet().iterator();
      while (iter.hasNext()) {
        Object key = iter.next();
        double prob = ( (Double) probs.get(key)).doubleValue();
        probs.put(key, new Double(prob / prob_sum));
      }
    }
    return probs;
  }

  /**
   * Calculate the likelihood function of the model parameters
   * @param components
   * @param weights
   * @param data
   * @return
   */
  protected double getParameterLogLikelihood(ComponentHashMap components, HashMap weights, double[][] data) throws ModelSingularException{
    double like = 0;
    for (int i = 0; i < data.length; i++) {
      double[] xk = data[i];
      double d = getDensity(components, weights, xk);
      // if density = 0, log will can -infinity.  Since likelihood is related to
      // stopping the iteration, a NaN liklihood is anyway not good.  A
      // temporary here is ignore the likelihood of that data.  NOTE: ignoring
      // it doesn't give the right answer.  Ignoring it is equivallent to
      // log(p) = 0 or p = 1.  p=1 is a very high density here vs. p=0 is the
      // lowest
      if (d <1e-300) d=1e-300;
      like += Math.log(d);
    }
    return like;
  }

  /**
   * Check the integrity of each components.  Fix if possible
   * @param new_components
   * @param old_components
   */
  protected void learnComponentModels(HashMap components) {
    int i_index = 0;
    int index_to_monitor = -1;
    for (Iterator i = components.keySet().iterator(); i.hasNext(); i_index++) {
      Object key = i.next();
      if (i_index == index_to_monitor) {
        i_index = index_to_monitor;
      }

      MixtureComponential component = (MixtureComponential)components.get(key);
      if (component instanceof ComponentialGMM) {
        System.out.println("learn Component model on " + key);
      }
      try {
        int pre_num_comp = 0;
        if (component instanceof ComponentialGMM) {
          pre_num_comp = ((ComponentialGMM)component).getRstNumComponents();
        }
        component.learnTheModel();
//        System.out.println(component);
        if (component instanceof ComponentialGMM) {
          int post_num_comp = ((ComponentialGMM)component).getRstNumComponents();
          if (pre_num_comp != post_num_comp) {
            System.out.println("component " + key + ", index = " + i_index + " has components removed");
          }
        }
      } catch (ModelSingularException e) {
        // e.printStackTrace();
        double weight = ((Double)m_weights.get(key)).doubleValue();
        double prob_sum = component.getProbabilitySum();
        System.err.print("\"" + key + "\"" + " singular ");
        final NumberFormat nf = NumberFormat.getInstance();
        nf.setMaximumFractionDigits(2);
        System.err.println("(weight = " + nf.format(weight * 100) + "%, prob_sum = " + nf.format(prob_sum) + ")" + ", removed");
        i.remove();
        continue;
      }
    }
  }

  public String toString() {
    StringBuffer str = new StringBuffer();
    // only generate the string decription when the model is calculated
    if (m_rstComponents != null) {
      // sort the object list to print the object in alphabetic order
      ArrayList keys = new ArrayList(m_rstComponents.keySet());
      Collections.sort(keys);
      for (Iterator i = keys.iterator(); i.hasNext(); ) {
        Object key = i.next();
        MixtureComponential component = (MixtureComponential) m_rstComponents.get(key);
        try {
          str.append("<<" + key + ">> " + (((Double)m_weights.get(key)).doubleValue() * 100) + "%\n");
        }
        catch (Exception e) {
          e.printStackTrace();
        }
        str.append(component.toString());
      }
    }
    return str.toString();
  }


  /**
   * ComponentHashMap contains a list of mixture components and its keys are
   * the name of those components.
   */
  public class ComponentHashMap extends HashMap {
    public void add(Object key, double[] value, double prob) {
      MixtureComponential component = (MixtureComponential)get(key);
      if (component == null) {  // if there is no corresponding component yet
        component = createEMComponent();
        put(key, component);
      }
      component.add(value, prob);
    }
  }


//  /** @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];
//
//      /** @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;
//
//  }

  /**
   * Classify the given data to the Gaussian model with highest p(data|g).
   * Consider the component weight by default.
   * @param data
   * @return
   */
  public Object classify(double[] data) {
    return classify(data, true);
  }

  /**
   * Classify the given data to the Gaussian model with highest p(data|g)
   * @param data to classified
   * @param comp_weighted take the component weight in account
   * @return
   */
  public Object classify(double[] data, boolean comp_weighted) {

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

    HashMap probs = getDensitiesUnderComponents(data, comp_weighted, false);

    Iterator prob_iter = probs.keySet().iterator();
    double max_prob = 0;
    Object max_prob_obj = null;
    while (prob_iter.hasNext()) {
      Object key = prob_iter.next();
      double prob = ( (Double) probs.get(key)).doubleValue(); // p(wi|xk)
      if (prob > max_prob) {
        max_prob = prob;
        max_prob_obj = key;
      }
    }

    return max_prob_obj;
  }

  /**
   * Convert the component (given by "comp") to the index in the resultant
   * components (beginning from 0)
   * @param comp
   * @return
   */
  public int objectToRstComponentIndex(Object comp) {
    ArrayList sorted_key = new ArrayList(m_rstComponents.keySet());
    Collections.sort(sorted_key);
    return sorted_key.indexOf(comp);
  }

  /**
   * Write the mixture model as a java object to file
   * @param file
   */
  public void writeToFile(String file) throws java.io.IOException {
    DOU_Ser.output(file, this);
  }

  /**
   * create a mixture model from a file
   *
   * @param file
   * @return ImgDataCenter
   */
  public static MixtureModel createFromFile(String file) throws Exception {
    return (MixtureModel) yi.util.io.input.DIU_Ser.input(file);
  }


}
