/* nag_mv_kmeans_cluster_analysis (g03efc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 * Mark 7, revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <nagx04.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagg03.h>

#define CMEANS(I, J) cmeans[(I) *tdcmeans + J]
#define X(I, J)      x[(I) *tdx + J]
int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
  Integer  exit_status = 0, i, *inc = 0, *isx = 0, j, k, m, maxit, n, *nic = 0,
           nvar;
  Integer  tdcmeans, tdx;
  NagError fail;
  char     weight[2];
  double   *cmeans = 0, *css = 0, *csw = 0, *wt = 0, *wtptr, *x = 0;

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  fprintf(fpout,
          "nag_mv_kmeans_cluster_analysis (g03efc) Example Program Results"
          "\n\n");

  /* Skip heading in the data file */
  fscanf(fpin, "%*[^\n]");

  fscanf(fpin, "%s", weight);
  fscanf(fpin, "%ld", &n);
  fscanf(fpin, "%ld", &m);
  fscanf(fpin, "%ld", &nvar);
  fscanf(fpin, "%ld", &k);
  fscanf(fpin, "%ld", &maxit);

  if (n >= 2 && nvar >= 1 && m >= nvar && k >= 2)
    {
      if (!(cmeans = NAG_ALLOC((k)*(nvar), double)) ||
          !(css = NAG_ALLOC(k, double)) ||
          !(csw = NAG_ALLOC(k, double)) ||
          !(wt = NAG_ALLOC(n, double)) ||
          !(x = NAG_ALLOC((n)*(m), double)) ||
          !(inc = NAG_ALLOC(n, Integer)) ||
          !(isx = NAG_ALLOC(m, Integer)) ||
          !(nic = NAG_ALLOC(k, Integer)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tdx = m;
      tdcmeans = nvar;
    }
  else
    {
      fprintf(fpout, "Invalid n or nvar or m or k.\n");
      exit_status = 1;
      return exit_status;
    }
  if (*weight == 'W')
    {
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            fscanf(fpin, "%lf", &X(i, j));
          fscanf(fpin, "%lf", &wt[i]);
        }
      wtptr = wt;
    }
  else
    {
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            fscanf(fpin, "%lf", &X(i, j));
        }
      wtptr = 0;
    }
  for (i = 0; i < k; ++i)
    {
      for (j = 0; j < nvar; ++j)
        fscanf(fpin, "%lf", &CMEANS(i, j));
    }
  for (j = 0; j < m; ++j)
    fscanf(fpin, "%ld", &isx[j]);

  /* nag_mv_kmeans_cluster_analysis (g03efc).
   * K-means
   */
  nag_mv_kmeans_cluster_analysis(n, m, x, tdx, isx, nvar, k, cmeans,
                                 tdcmeans, wtptr, inc, nic, css, csw, maxit,
                                 &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout,
              "Error from nag_mv_kmeans_cluster_analysis (g03efc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  fprintf(fpout, "\nThe cluster each point belongs to\n");
  for (i = 0; i < n; ++i)
    fprintf(fpout, " %6ld%s", inc[i], (i+1)%10?"":"\n");

  fprintf(fpout, "\n\nThe number of points in each cluster\n");
  for (i = 0; i < k; ++i)
    fprintf(fpout, " %6ld", nic[i]);

  fprintf(fpout, "\n\nThe within-cluster sum of weights of each cluster\n");
  for (i = 0; i < k; ++i)
    fprintf(fpout, " %9.2f", csw[i]);

  fprintf(fpout, "\n\nThe within-cluster sum of squares of each cluster\n\n");
  for (i = 0; i < k; ++i)
    fprintf(fpout, " %13.4f", css[i]);

  fprintf(fpout, "\n\nThe final cluster centres\n");
  fprintf(fpout, "                  1       2       3       4       5\n");
  for (i = 0; i < k; ++i)
    {
      fprintf(fpout, " %5ld     ", i+1);
      for (j = 0; j < nvar; ++j)
        fprintf(fpout, "%8.4f", CMEANS(i, j));
      fprintf(fpout, "\n");
    }
 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (cmeans) NAG_FREE(cmeans);
  if (css) NAG_FREE(css);
  if (csw) NAG_FREE(csw);
  if (wt) NAG_FREE(wt);
  if (x) NAG_FREE(x);
  if (inc) NAG_FREE(inc);
  if (isx) NAG_FREE(isx);
  if (nic) NAG_FREE(nic);
  return exit_status;
}