Graphviz  2.41.20170921.2350
pca.c
Go to the documentation of this file.
1 /* $Id$ $Revision$ */
2 /* vim:set shiftwidth=4 ts=8: */
3 
4 /*************************************************************************
5  * Copyright (c) 2011 AT&T Intellectual Property
6  * All rights reserved. This program and the accompanying materials
7  * are made available under the terms of the Eclipse Public License v1.0
8  * which accompanies this distribution, and is available at
9  * http://www.eclipse.org/legal/epl-v10.html
10  *
11  * Contributors: See CVS logs. Details at http://www.graphviz.org/
12  *************************************************************************/
13 
14 
15 #include "matrix_ops.h"
16 #include "pca.h"
17 #include "closest.h"
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21 
22 static int num_pairs = 4;
23 
24 void
25 PCA_alloc(DistType ** coords, int dim, int n, double **new_coords,
26  int new_dim)
27 {
28  double **DD = NULL; /* dim*dim matrix: coords*coords^T */
29  double sum;
30  int i, j, k;
31  double **eigs = NULL;
32  double *evals = NULL;
33  double *storage_ptr;
34 
35  eigs = N_GNEW(new_dim, double *);
36  for (i = 0; i < new_dim; i++)
37  eigs[i] = N_GNEW(dim, double);
38  evals = N_GNEW(new_dim, double);
39 
40  DD = N_GNEW(dim, double *);
41  storage_ptr = N_GNEW(dim * dim, double);
42  for (i = 0; i < dim; i++) {
43  DD[i] = storage_ptr;
44  storage_ptr += dim;
45  }
46 
47  for (i = 0; i < dim; i++) {
48  for (j = 0; j <= i; j++) {
49  /* compute coords[i]*coords[j] */
50  sum = 0;
51  for (k = 0; k < n; k++) {
52  sum += coords[i][k] * coords[j][k];
53  }
54  DD[i][j] = DD[j][i] = sum;
55  }
56  }
57 
58  power_iteration(DD, dim, new_dim, eigs, evals, TRUE);
59 
60  for (j = 0; j < new_dim; j++) {
61  for (i = 0; i < n; i++) {
62  sum = 0;
63  for (k = 0; k < dim; k++) {
64  sum += coords[k][i] * eigs[j][k];
65  }
66  new_coords[j][i] = sum;
67  }
68  }
69 
70  for (i = 0; i < new_dim; i++)
71  free(eigs[i]);
72  free(eigs);
73  free(evals);
74  free(DD[0]);
75  free(DD);
76 }
77 
78 boolean
79 iterativePCA_1D(double **coords, int dim, int n, double *new_direction)
80 {
81  vtx_data *laplacian;
82  float **mat1 = NULL;
83  double **mat = NULL;
84  double eval;
85 
86  /* Given that first projection of 'coords' is 'coords[0]'
87  compute another projection direction 'new_direction'
88  that scatters points that are close in 'coords[0]'
89  */
90 
91  /* find the nodes that were close in 'coords[0]' */
92  /* and construct appropriate Laplacian */
93  closest_pairs2graph(coords[0], n, num_pairs * n, &laplacian);
94 
95  /* Compute coords*Lap*coords^T */
96  mult_sparse_dense_mat_transpose(laplacian, coords, n, dim, &mat1);
97  mult_dense_mat_d(coords, mat1, dim, n, dim, &mat);
98  free(mat1[0]);
99  free(mat1);
100 
101  /* Compute direction */
102  return power_iteration(mat, dim, 1, &new_direction, &eval, TRUE);
103 /* ?? When is mat freed? */
104 }
105 
106 #ifdef UNUSED
107 
108 double dist(double **coords, int dim, int p1, int p2)
109 {
110  int i;
111  double sum = 0;
112 
113  for (i = 0; i < dim; i++) {
114  sum +=
115  (coords[i][p1] - coords[i][p2]) * (coords[i][p1] -
116  coords[i][p2]);
117  }
118  return sqrt(sum);
119 }
120 
121 
122 void weight_laplacian(double **X, int n, int dim, vtx_data * laplacian)
123 {
124  int i, j, neighbor;
125 
126  int *edges;
127  float *ewgts;
128  for (i = 0; i < n; i++) {
129  edges = laplacian[i].edges;
130  ewgts = laplacian[i].ewgts;
131  *ewgts = 0;
132  for (j = 1; j < laplacian[i].nedges; j++) {
133  neighbor = edges[j];
134  *ewgts -= ewgts[j] =
135  float (-1.0 / (dist(X, dim, i, neighbor) + 1e-10));
136  }
137  }
138 }
139 
140 #endif
boolean iterativePCA_1D(double **coords, int dim, int n, double *new_direction)
Definition: pca.c:79
int power_iteration(double **square_mat, int n, int neigs, double **eigs, double *evals, int initialize)
Definition: matrix_ops.c:24
void closest_pairs2graph(double *place, int n, int num_pairs, vtx_data **graph)
Definition: closest.c:354
int nedges
Definition: sparsegraph.h:80
void mult_sparse_dense_mat_transpose(vtx_data *A, double **B, int dim1, int dim2, float ***CC)
Definition: matrix_ops.c:214
int * edges
Definition: sparsegraph.h:81
int DistType
Definition: sparsegraph.h:92
void mult_dense_mat_d(double **A, float **B, int dim1, int dim2, int dim3, double ***CC)
Definition: matrix_ops.c:178
#define NULL
Definition: logic.h:39
float * ewgts
Definition: sparsegraph.h:82
void PCA_alloc(DistType **coords, int dim, int n, double **new_coords, int new_dim)
Definition: pca.c:25
double dist(Site *s, Site *t)
Definition: site.c:41
#define N_GNEW(n, t)
Definition: agxbuf.c:20
#define TRUE
Definition: cgraph.h:38