source: trunk/greenstone3-extensions/vishnu/src/cluster/mysammon.c@ 8189

Last change on this file since 8189 was 8189, checked in by kjdon, 20 years ago

first version of Imperial College's Visualiser code
  • Property svn:keywords set to Author Date Id Revision
File size: 6.0 KB
Line 
1/* - cannot have duplicated input values
2 * - check the init pt for Ys (cannot have duplicated pts)
3 * - err checking for opening files and mem alloc
4 * - check the input arg is less than MAX_N and MAX_L
5 */
6
7
8
9#include <stdio.h>
10#include <stdlib.h>
11#include <math.h>
12#include <assert.h>
13
14
15#define MF 0.3 /* Magic Factor */
16#define _MAX_Err 0.00005 /* Max Mapping Error */
17#define _MAX_M 1000 /* Max Iteration Number */
18#define _MAX_N 1000 /* Max Number of Entries */
19#define _MAX_L 128 /* Max Dimension for X */
20
21
22
23
24/* SS298: Added function for scaling values between 0 and 1
25 * n : number of entries
26 * mc99 fixed the errors that occured when maxX or MaxY ended up 0;
27 * fp libs and ops don't always know that 0/0 is 0
28 * More strictly this scales values to between .1 and 0.9
29 */
30void scale_values(int n, double **Ys)
31{
32
33 // get the maximum and the minimum of Ys array
34
35 // get the maximum and the minimum of Ys array
36
37 double maxX = Ys[0][0];
38 double minX = Ys[0][0];
39 double maxY = Ys[0][1];
40 double minY = Ys[0][1];
41
42 int i;
43 for (i = 1 ; i < n ; i++)
44 {
45
46 if (maxX < Ys[i][0]) maxX = Ys[i][0];
47 if (minX > Ys[i][0]) minX = Ys[i][0];
48 if (maxY < Ys[i][1]) maxY = Ys[i][1];
49 if (minY > Ys[i][1]) minY = Ys[i][1];
50 }
51 maxY -= minY;
52 maxX -= minX;
53 for (i = 0 ; i < n ; i++)
54 {
55 if (maxX > 0)
56 {
57 Ys[i][0] -= minX;
58 Ys[i][0] /= maxX;
59 Ys[i][0] *=0.8;
60 Ys[i][0] +=0.1;
61 }
62 else Ys[i][0]=0.5;
63 if (maxY > 0)
64 {
65 Ys[i][1] -= minY;
66 Ys[i][1] /= maxY;
67 Ys[i][1] *= 0.8;
68 Ys[i][1] += 0.1;
69 }
70 else Ys[i][1]=0.5;
71 }
72
73}
74
75
76/* return a Euclidean distance between two 'd'-dimensional vectors 'xs' and 'ys'
77 */
78double eucl_dist(int d, double *xs, double *ys)
79{
80 double ans=0.0;
81 double tmp;
82 int i;
83 for(i=0;i<d;i++)
84 {
85 tmp = xs[i]-ys[i];
86 ans += tmp * tmp;
87 }
88 return sqrt(ans);
89}
90
91
92/* n - number of entries
93 * d - dimension
94 * post:
95 * return an array of eucidean distances for X
96 * and assign its total distance (will be used in the mapping error) to 'total'
97 */
98double *compute_Xs_eucl_dists(double *total,int n, int d, double **Xs)
99{
100 double *tmp;
101 int i,j,k;
102
103 tmp = (double *) calloc(n*(n-1)/2, sizeof(double));
104
105
106 *total=0.0;
107 k=0;
108 for(i=1; i<n; i++)
109 for(j=0;j<i;j++)
110 {
111 tmp[k] = eucl_dist(d, Xs[i], Xs[j]);
112 *total+=tmp[k];
113 k++;
114 }
115 return tmp;
116}
117
118
119
120
121
122
123/*----------------------------------------------------------------*/
124void mysammon(double **Xs, int noe, int dim, char *buffer, int *cursize, int *curptr)
125{
126
127 int i, j, m, p, n;
128 char bbuf[20];
129 double *Xs_dists, total_Xs_dists;
130
131 double dist_star,dist, alpha, beta, gamma, zeta, err;
132 double uppDer0, uppDer1, lowDer0, lowDer1;
133 double **Ys;
134 double maxX=0.0;
135 double minX=0.0;
136 double maxY=0.0;
137 double minY=0.0;
138
139 Ys = calloc(noe,sizeof(double *));
140 assert(Ys);
141 for (i = 0;i<noe;i++)
142 {
143 Ys[i]=calloc(2,sizeof(double));
144 assert(Ys[i]);
145 }
146
147 Xs_dists = compute_Xs_eucl_dists(&total_Xs_dists, noe, dim, Xs);
148
149 p=0;
150 for (i=0;i<noe;i++)
151 {
152 Ys[i][0] = Xs[i][0];
153 Ys[i][1] = Xs[i][0]+Xs_dists[(i==(noe-1))?i:p];
154 p+=(noe-(i+1));;
155 }
156
157 m=0; /* iteration counter */
158 do{
159
160 for(p=0;p<noe;p++)
161 {
162 uppDer0=uppDer1=lowDer0=lowDer1=0.0;
163 /* calcuate partial derivatives */
164 for(j=0;j<noe;j++)
165 {
166 if (p==j) continue;
167 dist_star = Xs_dists[(j>p) ? (j*(j-1)/2+p) : (p*(p-1)/2+j)];
168
169 dist = eucl_dist(2,Ys[p],Ys[j]);
170 if (dist == 0) dist = 0.00001;
171 alpha = dist_star-dist;
172 beta = dist_star*dist;
173 if (beta == 0) beta = 0.00000000000001;
174 gamma = alpha/beta;
175
176 /* 1st dimension */
177 zeta = Ys[p][0]-Ys[j][0];
178 uppDer0 += gamma * zeta;
179 lowDer0 += (alpha-(zeta*zeta/dist) * (1+alpha/dist))/beta;
180
181 /* 2nd dimension */
182 zeta = Ys[p][1]-Ys[j][1];
183 uppDer1 += gamma * zeta;
184 lowDer1 += (alpha-(zeta*zeta/dist) * (1+alpha/dist))/beta;
185 } /* for */
186
187 Ys[p][0] = Ys[p][0]+MF*uppDer0/fabs(lowDer0);
188 Ys[p][1] = Ys[p][1]+MF*uppDer1/fabs(lowDer1);
189 if(p==0)
190 {
191 maxX = Ys[0][0];
192 minX = Ys[0][0];
193 maxY = Ys[0][1];
194 minY = Ys[0][1];
195 }
196 else
197 {
198 if (maxX < Ys[p][0]) maxX = Ys[p][0];
199 if (minX > Ys[p][0]) minX = Ys[p][0];
200 if (maxY < Ys[p][1]) maxY = Ys[p][1];
201 if (minY > Ys[p][1]) minY = Ys[p][1];
202 }
203
204 } /* for */
205
206 /* calcuate the mapping error */
207 err=0.0;
208 p=0;
209 for(j=1;j<noe;j++)
210 for(i=0;i<j;i++)
211 {
212 dist_star = Xs_dists[p++];
213 alpha = dist_star - eucl_dist(2,Ys[j],Ys[i]);
214 err += (alpha*alpha/dist_star);
215 }
216 err/=total_Xs_dists;
217
218 m++;
219 } while ((err>_MAX_Err) && (m!=_MAX_M));
220 // This is ugly desperate measures
221 srand(noe);
222 if(maxY==minY)
223 for(p=0;p<noe;p++)
224 Ys[p][1]=(double)rand()/(double)RAND_MAX;
225 if(maxX==minX)
226 for(p=0;p<noe;p++)
227 Ys[p][0]=(double)rand()/(double)RAND_MAX;;
228
229
230 scale_values(noe,Ys);
231
232 n = sprintf(bbuf,"%d\n", noe);
233 appendString(buffer, cursize, curptr, bbuf, strlen(bbuf));
234
235 for(i=0;i<noe;i++)
236 {
237 n = sprintf(bbuf,"%f,%f\n", Ys[i][0], Ys[i][1]);
238 appendString(buffer, cursize, curptr, bbuf, strlen(bbuf));
239 }
240
241 free(Xs_dists);
242 for (i = 0;i<noe;i++)
243 {
244 free(Ys[i]);
245 }
246 free(Ys);
247}
248
249
250
251
252
253
254
255
256
257
258
259
260
Note: See TracBrowser for help on using the repository browser.