diff maths/pca/pca.c @ 483:fdaa63607c15

Untabify, indent, tidy

--- a/maths/pca/pca.c	Fri May 31 11:02:28 2019 +0100
+++ b/maths/pca/pca.c	Fri May 31 11:54:32 2019 +0100
@@ -32,58 +32,58 @@
 /* Create m * m covariance matrix from given n * m data matrix. */
 void covcol(double** data, int n, int m, double** symmat)
 {
-double *mean;
-int i, j, j1, j2;
+    double *mean;
+    int i, j, j1, j2;
 
 /* Allocate storage for mean vector */
 
-mean = (double*) malloc(m*sizeof(double));
+    mean = (double*) malloc(m*sizeof(double));
 
 /* Determine mean of column vectors of input data matrix */
 
-for (j = 0; j < m; j++)
+    for (j = 0; j < m; j++)
     {
-    mean[j] = 0.0;
-    for (i = 0; i < n; i++)
+        mean[j] = 0.0;
+        for (i = 0; i < n; i++)
         {
-        mean[j] += data[i][j];
+            mean[j] += data[i][j];
         }
-    mean[j] /= (double)n;
+        mean[j] /= (double)n;
     }
 
 /*
-printf("\nMeans of column vectors:\n");
-for (j = 0; j < m; j++)  {
-    printf("%12.1f",mean[j]);  }   printf("\n");
- */
+  printf("\nMeans of column vectors:\n");
+  for (j = 0; j < m; j++)  {
+  printf("%12.1f",mean[j]);  }   printf("\n");
+*/
 
 /* Center the column vectors. */
 
-for (i = 0; i < n; i++)
+    for (i = 0; i < n; i++)
     {
-    for (j = 0; j < m; j++)
+        for (j = 0; j < m; j++)
         {
-        data[i][j] -= mean[j];
+            data[i][j] -= mean[j];
         }
     }
 
 /* Calculate the m * m covariance matrix. */
-for (j1 = 0; j1 < m; j1++)
+    for (j1 = 0; j1 < m; j1++)
     {
-    for (j2 = j1; j2 < m; j2++)
+        for (j2 = j1; j2 < m; j2++)
         {
-        symmat[j1][j2] = 0.0;
-        for (i = 0; i < n; i++)
+            symmat[j1][j2] = 0.0;
+            for (i = 0; i < n; i++)
             {
-            symmat[j1][j2] += data[i][j1] * data[i][j2];
+                symmat[j1][j2] += data[i][j1] * data[i][j2];
             }
-        symmat[j2][j1] = symmat[j1][j2];
+            symmat[j2][j1] = symmat[j1][j2];
         }
     }
 
-free(mean);
+    free(mean);
 
-return;
+    return;
 
 }
 
@@ -101,84 +101,84 @@
 /**  Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */
 
 /* Householder reduction of matrix a to tridiagonal form.
-Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
-Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
-Springer-Verlag, 1976, pp. 489-494.
-W H Press et al., Numerical Recipes in C, Cambridge U P,
-1988, pp. 373-374.  */
+   Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
+   Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
+   Springer-Verlag, 1976, pp. 489-494.
+   W H Press et al., Numerical Recipes in C, Cambridge U P,
+   1988, pp. 373-374.  */
 void tred2(double** a, int n, double* d, double* e)
 {
-	int l, k, j, i;
-	double scale, hh, h, g, f;
-	
-	for (i = n-1; i >= 1; i--)
+    int l, k, j, i;
+    double scale, hh, h, g, f;
+        
+    for (i = n-1; i >= 1; i--)
     {
-		l = i - 1;
-		h = scale = 0.0;
-		if (l > 0)
-		{
-			for (k = 0; k <= l; k++)
-				scale += fabs(a[i][k]);
-			if (scale == 0.0)
-				e[i] = a[i][l];
-			else
-			{
-				for (k = 0; k <= l; k++)
-				{
-					a[i][k] /= scale;
-					h += a[i][k] * a[i][k];
-				}
-				f = a[i][l];
-				g = f>0 ? -sqrt(h) : sqrt(h);
-				e[i] = scale * g;
-				h -= f * g;
-				a[i][l] = f - g;
-				f = 0.0;
-				for (j = 0; j <= l; j++)
-				{
-					a[j][i] = a[i][j]/h;
-					g = 0.0;
-					for (k = 0; k <= j; k++)
-						g += a[j][k] * a[i][k];
-					for (k = j+1; k <= l; k++)
-						g += a[k][j] * a[i][k];
-					e[j] = g / h;
-					f += e[j] * a[i][j];
-				}
-				hh = f / (h + h);
-				for (j = 0; j <= l; j++)
-				{
-					f = a[i][j];
-					e[j] = g = e[j] - hh * f;
-					for (k = 0; k <= j; k++)
-						a[j][k] -= (f * e[k] + g * a[i][k]);
-				}
-			}
-		}
-		else
-			e[i] = a[i][l];
-		d[i] = h;
+        l = i - 1;
+        h = scale = 0.0;
+        if (l > 0)
+        {
+            for (k = 0; k <= l; k++)
+                scale += fabs(a[i][k]);
+            if (scale == 0.0)
+                e[i] = a[i][l];
+            else
+            {
+                for (k = 0; k <= l; k++)
+                {
+                    a[i][k] /= scale;
+                    h += a[i][k] * a[i][k];
+                }
+                f = a[i][l];
+                g = f>0 ? -sqrt(h) : sqrt(h);
+                e[i] = scale * g;
+                h -= f * g;
+                a[i][l] = f - g;
+                f = 0.0;
+                for (j = 0; j <= l; j++)
+                {
+                    a[j][i] = a[i][j]/h;
+                    g = 0.0;
+                    for (k = 0; k <= j; k++)
+                        g += a[j][k] * a[i][k];
+                    for (k = j+1; k <= l; k++)
+                        g += a[k][j] * a[i][k];
+                    e[j] = g / h;
+                    f += e[j] * a[i][j];
+                }
+                hh = f / (h + h);
+                for (j = 0; j <= l; j++)
+                {
+                    f = a[i][j];
+                    e[j] = g = e[j] - hh * f;
+                    for (k = 0; k <= j; k++)
+                        a[j][k] -= (f * e[k] + g * a[i][k]);
+                }
+            }
+        }
+        else
+            e[i] = a[i][l];
+        d[i] = h;
     }
-	d[0] = 0.0;
-	e[0] = 0.0;
-	for (i = 0; i < n; i++)
+    d[0] = 0.0;
+    e[0] = 0.0;
+    for (i = 0; i < n; i++)
     {
-		l = i - 1;
-		if (d[i])
-		{
-			for (j = 0; j <= l; j++)
-			{
-				g = 0.0;
-				for (k = 0; k <= l; k++)
-					g += a[i][k] * a[k][j];
-				for (k = 0; k <= l; k++)
-					a[k][j] -= g * a[k][i];
-			}
-		}
-		d[i] = a[i][i];
-		a[i][i] = 1.0;
-		for (j = 0; j <= l; j++)
-			a[j][i] = a[i][j] = 0.0;
+        l = i - 1;
+        if (d[i])
+        {
+            for (j = 0; j <= l; j++)
+            {
+                g = 0.0;
+                for (k = 0; k <= l; k++)
+                    g += a[i][k] * a[k][j];
+                for (k = 0; k <= l; k++)
+                    a[k][j] -= g * a[k][i];
+            }
+        }
+        d[i] = a[i][i];
+        a[i][i] = 1.0;
+        for (j = 0; j <= l; j++)
+            a[j][i] = a[i][j] = 0.0;
     }
 }
 
@@ -186,171 +186,168 @@
 
 void tqli(double* d, double* e, int n, double** z)
 {
-	int m, l, iter, i, k;
-	double s, r, p, g, f, dd, c, b;
-	
-	for (i = 1; i < n; i++)
-		e[i-1] = e[i];
-	e[n-1] = 0.0;
-	for (l = 0; l < n; l++)
+    int m, l, iter, i, k;
+    double s, r, p, g, f, dd, c, b;
+        
+    for (i = 1; i < n; i++)
+        e[i-1] = e[i];
+    e[n-1] = 0.0;
+    for (l = 0; l < n; l++)
     {
-		iter = 0;
-		do
-		{
-			for (m = l; m < n-1; m++)
-			{
-				dd = fabs(d[m]) + fabs(d[m+1]);
-				if (fabs(e[m]) + dd == dd) break;
-			}
-			if (m != l)
-			{
-				if (iter++ == 30) erhand("No convergence in TLQI.");
-				g = (d[l+1] - d[l]) / (2.0 * e[l]);
-				r = sqrt((g * g) + 1.0);
-				g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
-				s = c = 1.0;
-				p = 0.0;
-				for (i = m-1; i >= l; i--)
-				{
-					f = s * e[i];
-					b = c * e[i];
-					if (fabs(f) >= fabs(g))
+        iter = 0;
+        do
+        {
+            for (m = l; m < n-1; m++)
+            {
+                dd = fabs(d[m]) + fabs(d[m+1]);
+                if (fabs(e[m]) + dd == dd) break;
+            }
+            if (m != l)
+            {
+                if (iter++ == 30) erhand("No convergence in TLQI.");
+                g = (d[l+1] - d[l]) / (2.0 * e[l]);
+                r = sqrt((g * g) + 1.0);
+                g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
+                s = c = 1.0;
+                p = 0.0;
+                for (i = m-1; i >= l; i--)
+                {
+                    f = s * e[i];
+                    b = c * e[i];
+                    if (fabs(f) >= fabs(g))
                     {
-						c = g / f;
-						r = sqrt((c * c) + 1.0);
-						e[i+1] = f * r;
-						c *= (s = 1.0/r);
+                        c = g / f;
+                        r = sqrt((c * c) + 1.0);
+                        e[i+1] = f * r;
+                        c *= (s = 1.0/r);
                     }
-					else
+                    else
                     {
-						s = f / g;
-						r = sqrt((s * s) + 1.0);
-						e[i+1] = g * r;
-						s *= (c = 1.0/r);
+                        s = f / g;
+                        r = sqrt((s * s) + 1.0);
+                        e[i+1] = g * r;
+                        s *= (c = 1.0/r);
                     }
-					g = d[i+1] - p;
-					r = (d[i] - g) * s + 2.0 * c * b;
-					p = s * r;
-					d[i+1] = g + p;
-					g = c * r - b;
-					for (k = 0; k < n; k++)
-					{
-						f = z[k][i+1];
-						z[k][i+1] = s * z[k][i] + c * f;
-						z[k][i] = c * z[k][i] - s * f;
-					}
-				}
-				d[l] = d[l] - p;
-				e[l] = g;
-				e[m] = 0.0;
-			}
-		}  while (m != l);
-	}
+                    g = d[i+1] - p;
+                    r = (d[i] - g) * s + 2.0 * c * b;
+                    p = s * r;
+                    d[i+1] = g + p;
+                    g = c * r - b;
+                    for (k = 0; k < n; k++)
+                    {
+                        f = z[k][i+1];
+                        z[k][i+1] = s * z[k][i] + c * f;
+                        z[k][i] = c * z[k][i] - s * f;
+                    }
+                }
+                d[l] = d[l] - p;
+                e[l] = g;
+                e[m] = 0.0;
+            }
+        }  while (m != l);
+    }
 }
 
 /* In place projection onto basis vectors */
 void pca_project(double** data, int n, int m, int ncomponents)
 {
-	int  i, j, k, k2;
-	double  **symmat, /* **symmat2, */ *evals, *interm;
-	
-	//TODO: assert ncomponents < m
-	
-	symmat = (double**) malloc(m*sizeof(double*));
-	for (i = 0; i < m; i++)
-		symmat[i] = (double*) malloc(m*sizeof(double));
-		
-	covcol(data, n, m, symmat);
-	
-	/*********************************************************************
-		Eigen-reduction
-		**********************************************************************/
-	
+    int  i, j, k, k2;
+    double  **symmat, /* **symmat2, */ *evals, *interm;
+        
+    //TODO: assert ncomponents < m
+        
+    symmat = (double**) malloc(m*sizeof(double*));
+    for (i = 0; i < m; i++)
+        symmat[i] = (double*) malloc(m*sizeof(double));
+                
+    covcol(data, n, m, symmat);
+        
+    /*********************************************************************
+                Eigen-reduction
+    **********************************************************************/
+        
     /* Allocate storage for dummy and new vectors. */
     evals = (double*) malloc(m*sizeof(double));     /* Storage alloc. for vector of eigenvalues */
     interm = (double*) malloc(m*sizeof(double));    /* Storage alloc. for 'intermediate' vector */
     //MALLOC_ARRAY(symmat2,m,m,double);    
-	//for (i = 0; i < m; i++) {
-	//	for (j = 0; j < m; j++) {
-	//		symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */
-	//	}
-	//}
+    //for (i = 0; i < m; i++) {
+    //      for (j = 0; j < m; j++) {
+    //              symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */
+    //      }
+    //}
     tred2(symmat, m, evals, interm);  /* Triangular decomposition */
-tqli(evals, interm, m, symmat);   /* Reduction of sym. trid. matrix */
+    tqli(evals, interm, m, symmat);   /* Reduction of sym. trid. matrix */
 /* evals now contains the eigenvalues,
-columns of symmat now contain the associated eigenvectors. */	
+   columns of symmat now contain the associated eigenvectors. */   
 
 /*
-	printf("\nEigenvalues:\n");
-	for (j = m-1; j >= 0; j--) {
-		printf("%18.5f\n", evals[j]); }
-	printf("\n(Eigenvalues should be strictly positive; limited\n");
-	printf("precision machine arithmetic may affect this.\n");
-	printf("Eigenvalues are often expressed as cumulative\n");
-	printf("percentages, representing the 'percentage variance\n");
-	printf("explained' by the associated axis or principal component.)\n");
-	
-	printf("\nEigenvectors:\n");
-	printf("(First three; their definition in terms of original vbes.)\n");
-	for (j = 0; j < m; j++) {
-		for (i = 1; i <= 3; i++)  {
-			printf("%12.4f", symmat[j][m-i]);  }
-		printf("\n");  }
- */
+  printf("\nEigenvalues:\n");
+  for (j = m-1; j >= 0; j--) {
+  printf("%18.5f\n", evals[j]); }
+  printf("\n(Eigenvalues should be strictly positive; limited\n");
+  printf("precision machine arithmetic may affect this.\n");
+  printf("Eigenvalues are often expressed as cumulative\n");
+  printf("percentages, representing the 'percentage variance\n");
+  printf("explained' by the associated axis or principal component.)\n");
+        
+  printf("\nEigenvectors:\n");
+  printf("(First three; their definition in terms of original vbes.)\n");
+  for (j = 0; j < m; j++) {
+  for (i = 1; i <= 3; i++)  {
+  printf("%12.4f", symmat[j][m-i]);  }
+  printf("\n");  }
+*/
 
 /* Form projections of row-points on prin. components. */
 /* Store in 'data', overwriting original data. */
-for (i = 0; i < n; i++) {
-	for (j = 0; j < m; j++) {
-		interm[j] = data[i][j]; }   /* data[i][j] will be overwritten */
+    for (i = 0; i < n; i++) {
+        for (j = 0; j < m; j++) {
+            interm[j] = data[i][j]; }   /* data[i][j] will be overwritten */
         for (k = 0; k < ncomponents; k++) {
-			data[i][k] = 0.0;
-			for (k2 = 0; k2 < m; k2++) {
-				data[i][k] += interm[k2] * symmat[k2][m-k-1]; }
+            data[i][k] = 0.0;
+            for (k2 = 0; k2 < m; k2++) {
+                data[i][k] += interm[k2] * symmat[k2][m-k-1]; }
         }
-}
+    }
 
-/*	
-printf("\nProjections of row-points on first 3 prin. comps.:\n");
- for (i = 0; i < n; i++) {
-	 for (j = 0; j < 3; j++)  {
-		 printf("%12.4f", data[i][j]);  }
-	 printf("\n");  }
- */
+/*      
+        printf("\nProjections of row-points on first 3 prin. comps.:\n");
+        for (i = 0; i < n; i++) {
+        for (j = 0; j < 3; j++)  {
+        printf("%12.4f", data[i][j]);  }
+        printf("\n");  }
+*/
 
 /* Form projections of col.-points on first three prin. components. */
 /* Store in 'symmat2', overwriting what was stored in this. */
 //for (j = 0; j < m; j++) {
-//	 for (k = 0; k < m; k++) {
-//		 interm[k] = symmat2[j][k]; }  /*symmat2[j][k] will be overwritten*/
+//       for (k = 0; k < m; k++) {
+//               interm[k] = symmat2[j][k]; }  /*symmat2[j][k] will be overwritten*/
 //  for (i = 0; i < 3; i++) {
-//	symmat2[j][i] = 0.0;
-//		for (k2 = 0; k2 < m; k2++) {
-//			symmat2[j][i] += interm[k2] * symmat[k2][m-i-1]; }
-//		if (evals[m-i-1] > 0.0005)   /* Guard against zero eigenvalue */
-//			symmat2[j][i] /= sqrt(evals[m-i-1]);   /* Rescale */
-//		else
-//			symmat2[j][i] = 0.0;    /* Standard kludge */
+//      symmat2[j][i] = 0.0;
+//              for (k2 = 0; k2 < m; k2++) {
+//                      symmat2[j][i] += interm[k2] * symmat[k2][m-i-1]; }
+//              if (evals[m-i-1] > 0.0005)   /* Guard against zero eigenvalue */
+//                      symmat2[j][i] /= sqrt(evals[m-i-1]);   /* Rescale */
+//              else
+//                      symmat2[j][i] = 0.0;    /* Standard kludge */
 //    }
 // }
 
 /*
- printf("\nProjections of column-points on first 3 prin. comps.:\n");
- for (j = 0; j < m; j++) {
-	 for (k = 0; k < 3; k++)  {
-		 printf("%12.4f", symmat2[j][k]);  }
-	 printf("\n");  }
-	*/
+  printf("\nProjections of column-points on first 3 prin. comps.:\n");
+  for (j = 0; j < m; j++) {
+  for (k = 0; k < 3; k++)  {
+  printf("%12.4f", symmat2[j][k]);  }
+  printf("\n");  }
+*/
 
 
-for (i = 0; i < m; i++)
-	free(symmat[i]);
-free(symmat);
+    for (i = 0; i < m; i++)
+        free(symmat[i]);
+    free(symmat);
 //FREE_ARRAY(symmat2,m);
-free(evals);
-free(interm);
+    free(evals);
+    free(interm);
 
 }
-
-
-