#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          nclust = 0;
    long          ndiff = 0;
    int           ctype = 0;
    long          k = 3;
    double       *dendist = 0;
    double       *mergedist = 0;
    long         *counts = 0;
    long         *ic = 0;
    long         *nic = 0;
    long         *indc = 0;
    long         *lmerge = 0;
    long         *umerge = 0;
    long         *denord = 0;
    long          bcat = 0;
    int           info = 0;
    FILE         *fp = 0;
    long          i, j;

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

    fclose(fp);

    /*
      Summary information.
    */
    printf(" Number of Observations = %ld \n",nrec);
    printf(" Number of Variables    = %ld \n",nvar);
    if (ctype == 0)
        printf(" Goup Average Method\n\n");
    else
        printf(" Minimum Variance Method\n\n");

    /*
      Compute the hierarchical clustering.
    */
    nagdmc_hclust(rec1,nvar,nrec,dblk,data,nxvar,xvar,ctype,lmerge,umerge,
                  mergedist,denord,dendist,&info);
    if (handle_info("nagdmc_hclust",info)) {
        free(data);
        free(dendist);
        free(mergedist);
        free(denord);
        free(lmerge);
        free(umerge);
        free(ic);
        free(indc);        
        free(counts); 
        return 2;
    }


    /*
      Print results of the hierarchical clustering.
    */
    printf(" No. of Clusters\t\tDistance\n\n");
    for (i=nrec-10; i<nrec-1; ++i) 
        printf("\t%li\t\t%f \n",nrec-i-1,mergedist[i]);


    /*
      Assign data records to a given number of groups.
    */
    nagdmc_cind(nrec,mergedist,denord,dendist,k,ic,nic,&nclust,&info);
    if (handle_info("nagdmc_cind",info)) {
        free(data);
        free(dendist);
        free(mergedist);
        free(denord);
        free(lmerge);
        free(umerge);
        free(ic);
        free(indc);        
        free(counts);
        free(nic);
        return 2;
    }

    
    /*
      Cross-tabulate the groupings and the classes.
    */
    if (nclust == k) {
        nagdmc_tab2(nrec,k,bcat,ic,indc,counts,&ndiff,&info);

        if (handle_info("nagdmc_tab2",info)) {
            free(data);
            free(dendist);
            free(mergedist);
            free(denord);
            free(lmerge);
            free(umerge);
            free(ic);
            free(indc);        
            free(counts); 
            free(nic);
            return 2;
        }
        
        printf("\nTable 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");
        }
    }
    else
        printf(" Different number of clusters found (%ld).\n\n",nclust);


    (void)printf("\n  %6s\t%6s\t%10s\t%9s\t%8s\n",
                 "LMerge","UMerge","Merge Dist","Den Order","Den dist");
    for (i=0; i<nrec-1; i++) {
        (void)printf("  %-6li\t%-6li\t%-10.4f\t%-9li\t%-8.4f\n",
                     lmerge[i],umerge[i],mergedist[i],denord[i],dendist[i]);
    }
    (void)printf("  %-6s\t%-6s\t%-10s\t%-9li\t%-8.4f\n\n",
                 "N/A","N/A","N/A",denord[i],dendist[i]);

    /*
      Return allocated memory to the operating system.
    */
    free(dendist);
    free(mergedist);
    free(data);
    free(denord);
    free(lmerge);
    free(umerge);
    free(ic);
    free(indc);
    free(counts);
    free(nic);
    
    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;
}