#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          chunksize = 0;
    long          iwts = -1;
    long          nscores = 2;
    long          nzw;
    int           pcatype = 2;
    long          std = 0;
    double       *loadings = 0;
    double       *results = 0;
    double       *s = 0;
    double       *scores = 0;
    double       *xbar = 0;
    double        wsum = 0.0;
    int           info = 0;
    FILE         *fp = 0;
    long          i, j;

    /*
      Allocate memory.
    */
    if (!(data = (double *)malloc(dblk*nvar * sizeof(double))) ||
        !(loadings = (double *)malloc(nvar*nvar * sizeof(double))) ||
        !(results = (double *)malloc(6*nvar * sizeof(double))) ||
        !(scores = (double *)malloc(nvar * sizeof(double))) ||
        !(xbar = (double *)malloc(nvar * sizeof(double)))) {
        free(data);
        free(loadings);
        free(results);
        free(scores);
        printf (" Memory allocation failure.\n\n");
        return 2;
    }
    
    if (pcatype == 2) {
        if (!(s = (double *)malloc(nvar * sizeof(double)))) {
            free(data);
            free(loadings);
            free(results);
            free(scores);
            free(xbar);
            printf (" Memory allocation failure.\n\n");
            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(loadings);
        free(results);
        free(scores);
        free(xbar);
        if (s) free(s);
        return 2;
    }
    
    for (i=0; i<dblk; ++i) {
        for (j=0; j<nvar; ++j) 
            fscanf(fp,"%lf ",&data[i*nvar+j]);
        fscanf(fp,"%*d");
    }

    fclose(fp);

    /*
      Print task summary information.
    */
    printf("\n Number of Observations = %ld.\n",nrec);
    printf(" Number of Variables    = %ld.\n",nvar);
    if (pcatype < 1)
        printf(" Using Sum of Squares Matrix.\n");
    else if (pcatype == 1)
         printf(" Using Covariance Matrix.\n");
    else if (pcatype == 2)
        printf(" Using Correlation Matrix.\n");
    else
        printf(" Using User Standardisation.\n");
    if (iwts >= 0)
        printf(" Using Weights.\n\n");
    

    /*
      Compute means and standard deviations of data values.
    */
    nagdmc_dsu(rec1,nvar,nrec,dblk,data,0,0,chunksize,iwts,xbar,s,&wsum,&nzw,
               &info);
    if (handle_info("nagdmc_dsu",info)) {
        free(data);
        free(loadings);
        free(results);
        free(scores);
        free(xbar);
        if (s) free(s);
        return 2;
    }

    /*
      Compute the PCA.
    */
    nagdmc_pca(rec1,nvar,nrec,dblk,data,0,0,chunksize,iwts,pcatype,xbar,s,
               loadings,results,&info);
    if (handle_info("nagdmc_pca",info)) {
        free(data);
        free(loadings);
        free(results);
        free(scores);
        free(xbar);
        if (s) free(s);
        return 2;
    }

    /*
      Print the loadings and summary results from the PCA.
    */
    printf(" Loadings \n\n");
    for (i=0; i<nvar; ++i) 
    {
        for (j=0; j<nvar; ++j) 
            printf("%10f ",loadings[i*nvar+j]);
        printf(" \n");
    }
    
    printf(" Results \n\n");
    for (i=0; i<nvar; ++i)
    {
        for (j=0; j<6; ++j)
            printf("%10f ",results[6*i+j]);
        printf("\n");
    }

    /*
      Compute and print PCA scores using "nscores" components.
    */
    printf(" PCA Scores for %li Component Directions\n\n",nscores);
    for (i=0; i<10; ++i) {
        nagdmc_pca_score(nvar,nscores,std,&data[i*nvar],xbar,s,loadings,
                         results,scores,&info);

        if (handle_info("nagdmc_pca_score",info)) {
            free(data);
            free(loadings);
            free(results);
            free(scores);
            free(xbar);
            if (s) free(s);
            return 2;
        }

        for (j=0; j<nscores; ++j)
            printf(" %8.4f ",scores[j]);
        printf("\n");
    }

    /*
      Return allocated memory to operating system.
    */
    free(data);
    free(loadings);
    free(results);
    free(scores);
    free(xbar);
    if (s) free(s);
        
    return 0;
}

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

    return 0;
}