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
22static int num_pairs = 4;
23
24void
25PCA_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
78boolean
79iterativePCA_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
108double 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
122void 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