最大信息系数（maximal information coefficient，MIC）核心程序之OptimizeXAxis

/*
 * Returns the normalized MI scores.
 *
 * Parameters
 *   dx (IN) : x-data sorted in increasing order by dx-values
 *   dy (IN) : y-data sorted in increasing order by dx-values
 *   n (IN) : number of elements in dx and dy
 *   Q_map (IN) : the map Q computed by EquipartitionYAxis() sorted in
 *                increasing order by dx-values
 *   q (IN) : number of partitions in Q_map
 *   P_map (IN) : the map P computed by GetSuperclumpsPartition() sorted
 *                in increasing order by Dx-values
 *   p (IN) : number of partitions in P_map
 *   x (IN) : maximum grid size on dx-values
 *   score (OUT) : mutual information scores. score must be a
 *                 preallocated array of dimension x-1
 * Returns
 *   0 on success, 1 if an error occurs
 */
int OptimizeXAxis(double *dx, double *dy, int n, int *Q_map, int q,
      int *P_map, int p, int x, double *score)
{
  int i, s, t, l;
  int *c;
  int **cumhist;
  double **I, **HP2Q;
  double F, F_max, HQ, ct, cs;
  double **cumhist_log, *c_log;

  /* return score=0 if p=1 */
  if (p == 1)
    {
      for (i=0; i<x-1; i++)
        score[i] = 0.0;
      return 0;
    }

  /* compute c */
  c = compute_c(P_map, p, n);
  if (c == NULL)
    goto error_c;

  /* precompute log(c) (log(c)=0 when c=0) */
  c_log = compute_c_log(c, p);
  if (c_log == NULL)
    goto error_c_log;

  /* compute the cumulative histogram matrix along P_map */
  cumhist = compute_cumhist(Q_map, q, P_map, p, n);
  if (cumhist == NULL)
    goto error_cumhist;

  /* precompute log(cumhist) (log(cumhist)=0 when cumhist=0) */
  cumhist_log = compute_cumhist_log(cumhist, q, p);
  if (cumhist == NULL)
    goto error_cumhist_log;

  /* I matrix initialization */
  I = init_I(p, x);
  if (I == NULL)
    goto error_I;

  /* Precomputes the HP2Q matrix */
  HP2Q = compute_HP2Q(cumhist, c, q, p);
  if (HP2Q == NULL)
    goto error_HP2Q;

  /* compute H(Q) */
  HQ = hq(cumhist, cumhist_log, q, p, n);

  /* Find the optimal partitions of size 2, Algorithm 2 in SOM, lines 3-8 */
  for (t=2; t<=p; t++)
    {
      F_max = -DBL_MAX;
      for (s=1; s<=t; s++)
        {
          F = hp3(c, c_log, s, t) - hp3q(cumhist, cumhist_log, c, q, p, s, t);
          if (F > F_max)
            {
              I[t][2] = HQ + F;
              F_max = F;
            }
        }
    }

  /*
   * Inductively build the rest of the table of optimal partitions,
   * Algorithm 2 in SOM, lines 10-17
   */
  for (l=3; l<=x; l++)
    {
      for (t=l; t<=p; t++)
        {
          ct = (double) c[t-1];
          F_max = -DBL_MAX;
          for (s=l-1; s<=t; s++)
            {
              cs = (double) c[s-1];
              F = ((cs/ct) * (I[s][l-1]-HQ)) - (((ct-cs)/ct) * HP2Q[s][t]);

              if (F > F_max)
                {
                  I[t][l] = HQ + F;
                  F_max = F;
                }
            }
        }
    }

  /* Algorithm 2 in SOM, line 19 */
  for (i=p+1; i<=x; i++)
    I[p][i] = I[p][p];

  /* score */
  for (i=2; i<=x; i++)
    score[i-2] = I[p][i] / MIN(log(i), log(q));

  /* start frees */
  for (i=0; i<=p; i++)
    free(HP2Q[i]);
  free(HP2Q);

  for (i=0; i<=p; i++)
    free(I[i]);
  free(I);

  for (i=0; i<q; i++)
    free(cumhist_log[i]);
  free(cumhist_log);

  for (i=0; i<q; i++)
    free(cumhist[i]);
  free(cumhist);

  free(c_log);
  free (c);
  /* end frees */

  return 0;

  /* gotos */
  error_HP2Q:
    for (i=0; i<=p; i++)
      free(I[i]);
    free(I);
  error_I:
    for (i=0; i<q; i++)
      free(cumhist_log[i]);
    free(cumhist_log);
  error_cumhist_log:
    for (i=0; i<q; i++)
      free(cumhist[i]);
    free(cumhist);
  error_cumhist:
    free(c_log);
  error_c_log:
    free(c);
  error_c:
    return 1;
}