#include <nagdmc.h>

/*
  handle_info() prints to screen information based on the value of the info
  parameter.
*/
int
handle_info(const char func[], int info);

int
main(void) {
    const char    file[] = {"iris.dat"};
    long          rec1 = 0;
    long          nvar = 4;
    long          nrec = 150;
    long          dblk = 150;
    double       *data = 0;
    long          nxvar = 0;
    long         *xvar = 0;
    long          iwts = -1;
    long          k = 3;
    long          maxit = 1000;
    int           info = 0;
    long          chunksize = 30;
    long          ndiff = 0;
    double       *cmeans = 0;
    double       *css = 0;
    double        tss = 0.0;
    long         *ic = 0;
    long         *nic = 0;
    long         *indc = 0;
    long         *counts = 0;
    long          seed[5];
    long          bcat = 0;
    FILE         *fp = 0;
    long          i, j, l;

    /*
      Allocate memory.
    */
    if (!(data = (double *)malloc(dblk*nvar * sizeof(double))) ||
        !(cmeans = (double *)malloc(k*nvar * sizeof(double))) ||
        !(counts = (long *)malloc(k*k * sizeof(long))) ||
        !(css = (double *)malloc(sizeof(double)*(k))) ||
        !(ic = (long *)malloc(nrec * sizeof(long))) ||
        !(indc = (long *)malloc(dblk * sizeof(long))) ||
        !(nic = (long *)malloc(k * sizeof(long)))) {
        printf(" Memory allocation failure.\n\n");
        free(data);
        free(cmeans);
        free(counts);
        free(css);
        free(ic);
        free(indc);
        return 2;
    }
    
    /*
      Read data values.
    */
    if ((fp = fopen(file,"r")) == 0) {
        printf(" Data file named %s was not found.\n\n",file);
        free(data);
        free(cmeans);
        free(counts);
        free(css);
        free(ic);
        free(indc);
        free(nic);
        return 2;
    }
    
    for (i=0; i<dblk; ++i) {
        for (j=0; j<nvar; ++j) 
            fscanf(fp,"%lf ",&data[i*nvar+j]);
        fscanf(fp,"%li\n",&indc[i]);
    }

    fclose(fp);

    /*
      Print summary information.
    */
    printf("\n Number of Observations = %ld \n", nrec);
    printf(" Number of Variables    = %ld \n", nvar);
    printf(" Number of Clusters     = %ld \n", k);

    /*
      Set initial clusters by using a repeatable random seed.
    */
    nagdmc_srs(125793,seed);

    nagdmc_rints(0,k,rec1,rec1+nrec,ic,seed,&info);

    if (handle_info("nagdmc_rints",info)) {
        free(data);
        free(cmeans);
        free(css);
        free(nic);
        free(ic);
        free(indc);
        free(counts);
        return 2;
    }
    
    for (j=0; j<k; ++j) {
        l = ic[j];
        
        for (i=0; i<nvar; ++i) 
            cmeans[j*nvar+i] = data[l*nvar+i];
    }
    
    nagdmc_nrgp(rec1,nvar,nrec,dblk,data,0,0,chunksize,nxvar,xvar,k,cmeans,
                ic,&info);

    if (handle_info("nagdmc_nrgp",info)) {
        free(data);
        free(cmeans);
        free(css);
        free(nic);
        free(ic);
        free(indc);
        free(counts);
        return 2;
    }

    nagdmc_kmeans(rec1,nvar,nrec,dblk,data,nxvar,xvar,iwts,k,ic,cmeans,maxit,
                  &info);

    if (handle_info("nagdmc_kmeans",info)) {
        free(data);
        free(cmeans);
        free(css);
        free(nic);
        free(ic);
        free(indc);
        free(counts);
        return 2;
    }

    nagdmc_wcss(rec1,nvar,nrec,dblk,data,0,0,chunksize,nxvar,xvar,iwts,k,cmeans,
                ic,css,nic,&tss,&info);

    if (handle_info("nagdmc_wcss",info)) {
        free(data);
        free(cmeans);
        free(css);
        free(nic);
        free(ic);
        free(indc);
        free(counts);
        return 2;
    }

    /*
      Print results of clustering.
    */
    printf(" Within Cluster Sum of Squares = %f.\n",tss);
    printf(" The Number of Data Records in Clusters:\n\n");
    for (i=0; i<k; ++i) 
        printf(" %li ",nic[i]);
    printf("\n\n");

    /*
      Cross-tabulate groupings with the observed values.
    */
    nagdmc_tab2(nrec,k,bcat,ic,indc,counts,&ndiff,&info);

    if (handle_info("nagdmc_tab2",info)) {
        free(data);
        free(cmeans);
        free(css);
        free(nic);
        free(ic);
        free(indc);
        free(counts);
        return 2;
    }

    /*
      Print tabulation.
    */
    printf(" Table of Counts\n\n");
    for (i=0; i<k; ++i) {
        for (j=0; j<k; j++)
            printf("\t%ld",counts[i*k+j]);
        printf("\n");
    }
    
    /*
      Return allocated memory to the operating system.
    */
    free(data);
    free(cmeans);
    free(css);
    free(nic);
    free(ic);
    free(indc);
    free(counts);

    return 0;
}

int
handle_info(const char func[], int info) {
    if (info == -999)
    {
        printf(" Invalid licence, please contact NAG.\n\n");
        return 2;
    }
    else if (info > 0)
    {
        printf(" Error code %i from %s.\n\n",info,func);
        return 1;
    }
    else if (info < 0)
        printf (" Information code %i from %s.\n\n",info,func);

    return 0;
}