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 |
141 | |