1 /* Last changed Time-stamp: <95/07/12 19:49:53 ivo>*/
3 Split Decomposition of Distance Matrices
4 as described by H.J.Bandelt and A.W.M.Dress
5 Adv Math, 92 (1992) p47
6 c Peter Stadler and Ivo L. Hofacker
14 #define PRIVATE static
16 #define SHORT_OUTPUT 1
25 double isolation_index; } Split;
27 PUBLIC Split *split_decomposition(float **dist);
28 PUBLIC void free_Split(Split *x);
29 PUBLIC void print_Split(Split *x);
30 PUBLIC void sort_Split(Split *x);
32 PUBLIC Split *split_decomposition(float **dist)
37 int i,j,sp,new_sp,spp;
40 double alpha,beta,tmp;
45 number_of_points = (int) dist[0][0];
50 maxlen = (number_of_points+1)*(number_of_points+1);
51 full_slots = (short *) space(2*(maxlen+1)*sizeof(short));
54 SD = space((maxlen+1)*sizeof(Split));
56 SD[0].splitlist[0] = NULL;
57 SD[0].splitlist[1] = NULL;
59 SD[0].isolation_index = 0.0;
61 SD[1].splitlist[0] = (short *) space((number_of_points+1)*sizeof(short));
62 SD[1].splitlist[0][0] = number_of_points;
63 SD[1].splitlist[0][1] = 1;
64 SD[1].splitlist[1] = (short *) space((number_of_points+1)*sizeof(short));
65 SD[1].splitlist[1][1] = 2;
67 SD[1].isolation_index = dist[1][2];
72 for( elm=3; elm <= number_of_points; elm++){
74 for (sp=1, spp=n_of_splits;sp<=n_of_splits; sp++){
76 SD[sp].splitlist[0][0] = elm;
78 /* 1. Split = {old_A | old_B+{elm}} */
80 alpha = 2.0*SD[sp].isolation_index;
81 for(i2=1; i2<=(elm - SD[sp].splitsize); i2++){
82 if(i2==(elm-SD[sp].splitsize)) {
86 x= SD[sp].splitlist[1][i2];
88 for(j1=1; j1 <= SD[sp].splitsize; j1++){
89 y= SD[sp].splitlist[0][j1];
90 for(j2=1; j2 <= (SD[sp].splitsize); j2++){
91 z= SD[sp].splitlist[0][j2];
93 /* calculate the value beta = beta(elm,x; y,z) */
95 beta = dist[elm][y] + dist[x][z];
96 tmp = dist[elm][z] + dist[x][y];
97 if(tmp>beta) beta=tmp;
98 tmp = dist[elm][x] + dist[y][z];
99 if(tmp>beta) beta=tmp;
100 beta -= ( dist[elm][x] + dist[y][z] );
102 if(beta<alpha) alpha=beta;
109 /* add {old_A | old_B + {elm}} to the split-list */
112 SD[spp].splitsize = SD[sp].splitsize;
113 SD[spp].isolation_index = alpha;
114 SD[spp].splitlist[0] = (short *) space((number_of_points+1)*sizeof(short));
115 SD[spp].splitlist[0][0] = elm;
116 SD[spp].splitlist[1] = (short *) space((number_of_points+1)*sizeof(short));
117 for(i=1; i<=SD[spp].splitsize; i++)
118 SD[spp].splitlist[0][i] = SD[sp].splitlist[0][i];
119 for(i=1; i<=(elm-1-SD[spp].splitsize); i++)
120 SD[spp].splitlist[1][i] = SD[sp].splitlist[1][i];
121 SD[spp].splitlist[1][elm-SD[spp].splitsize] = elm;
124 /* 2. Split = {old_A+{elm} | old_B} */
126 alpha = 2.0*SD[sp].isolation_index;
127 for(i2=1; i2<= (SD[sp].splitsize+1); i2++){
128 if(i2==(SD[sp].splitsize+1)){
133 x= SD[sp].splitlist[0][i2];
135 for(j1=1; j1 <= (elm-1-SD[sp].splitsize); j1++){
136 y= SD[sp].splitlist[1][j1];
137 for(j2=1; j2 <= (elm-1-SD[sp].splitsize); j2++){
138 z= SD[sp].splitlist[1][j2];
140 /* calculate the value beta = beta(elm,x; y,z) */
142 beta = dist[elm][y] + dist[x][z];
143 tmp = dist[elm][z] + dist[x][y];
144 if(tmp>beta) beta=tmp;
145 tmp = dist[elm][x] + dist[y][z];
146 if(tmp>beta) beta=tmp;
147 beta -= ( dist[elm][x] + dist[y][z] );
149 if(beta<alpha) alpha=beta;
156 /* replace {old_A | old_B} by {old_A+{elm} | old_B} in the splitlist */
158 SD[sp].splitlist[0][SD[sp].splitsize] = elm;
159 SD[sp].isolation_index = alpha;
162 /* remove {old_A | old_B} from the split list */
164 free(SD[sp].splitlist[0]);
165 SD[sp].splitlist[0] = NULL;
166 free(SD[sp].splitlist[1]);
167 SD[sp].splitlist[1] = NULL;
168 SD[sp].splitsize = 0;
169 SD[sp].isolation_index = 0.0;
173 /* 3. Split = { {elm} | {1,...,elm-1} } */
176 for(i=1;i<=elm-1;i++){
177 for(j=1; j<=elm-1;j++){
178 tmp = dist[elm][i]+dist[elm][j] - dist[i][j];
179 if( tmp < alpha) alpha = tmp;
186 SD[spp].splitsize = 1;
187 SD[spp].isolation_index = alpha;
188 SD[spp].splitlist[0] = (short *) space((number_of_points+1)*sizeof(short));
189 SD[spp].splitlist[0][0] = elm;
190 SD[spp].splitlist[1] = (short *) space((number_of_points+1)*sizeof(short));
191 SD[spp].splitlist[0][1] = elm;
192 for (i=1; i<= elm-1; i++)
193 SD[spp].splitlist[1][i] = i;
195 /* note that spp now points to the last entry in the splitlist !!! */
197 /* garbage collection .... unfortunately a bit boring */
199 for ( sp=1, new_sp=0; sp<= spp; sp++ ){
203 /* copy all the junk from sp to new_sp */
204 full_slots[new_sp]=1; full_slots[sp]=0;
205 SD[new_sp].splitsize = SD[sp].splitsize;
206 SD[new_sp].isolation_index = SD[sp].isolation_index;
207 if(SD[new_sp].splitlist[0]==NULL)
208 SD[new_sp].splitlist[0] = (short *) space((number_of_points+1)*sizeof(short));
209 if(SD[new_sp].splitlist[1]==NULL)
210 SD[new_sp].splitlist[1] = (short *) space((number_of_points+1)*sizeof(short));
212 SD[new_sp].splitlist[0][i] = SD[sp].splitlist[0][i];
213 SD[new_sp].splitlist[1][i] = SD[sp].splitlist[1][i];
215 free(SD[sp].splitlist[0]);
216 SD[sp].splitlist[0] = NULL;
217 free(SD[sp].splitlist[1]);
218 SD[sp].splitlist[1] = NULL;
219 SD[sp].splitsize = 0;
220 SD[sp].isolation_index = 0.0;
224 n_of_splits = new_sp;
225 SD[0].splitsize = n_of_splits;
227 for(test1=0, i=2; i<=elm; i++) for( j=1; j<i; j++) test1+=dist[i][j];
228 for(test2=0, i=1; i<= n_of_splits; i++)
229 test2 += (elm-SD[i].splitsize)*SD[i].splitsize*SD[i].isolation_index;
230 SD[0].isolation_index = (test1 - test2)/test1;
233 } /* End of iteration */
236 /* Calculate fraction of split-prime part for the full matrix
237 if not already done */
239 for(test1=0.0, i=2; i<=number_of_points; i++)
240 for( j=1; j<i; j++) test1+=dist[i][j];
241 for(test2=0.0, i=1; i<= n_of_splits; i++)
242 test2 += (number_of_points-SD[i].splitsize)*
243 SD[i].splitsize*SD[i].isolation_index;
244 SD[0].isolation_index = (test1 - test2)/test1;
249 /* free superfluous space */
251 S = space((n_of_splits+1)*sizeof(Split));
252 for (i=0; i<= n_of_splits; i++) S[i] = SD[i];
261 /* -------------------------------------------------------------------------- */
263 PUBLIC void free_Split(Split *x)
266 for (i=0; i<=x[0].splitsize; i++) {
267 if( x[i].splitlist[0] != NULL ) free( x[i].splitlist[0] );
268 if( x[i].splitlist[1] != NULL ) free( x[i].splitlist[1] );
273 /* -------------------------------------------------------------------------- */
275 PUBLIC void print_Split(Split *x)
278 printf("> %d Split Decomposition",x[0].splitsize);
279 for (k=1; k<=x[0].splitsize; k++) {
280 if( (x[k].splitlist[0])&&(x[k].splitlist[1]) ) {
281 printf("\n%3d %8.4f : {",k,x[k].isolation_index);
282 for(j=1;j<=x[k].splitsize;j++)
283 printf(" %3d",x[k].splitlist[0][j]);
288 for(j=1;j<=(x[k].splitlist[0][0]-x[k].splitsize);j++)
289 printf(" %3d",x[k].splitlist[1][j]);
294 printf("\n %8.4f : { [Split prime fraction] }\n", x[0].isolation_index);
298 /* -------------------------------------------------------------------------- */
300 int CompareSplit(const void *v1 , const void *v2)
304 S1 = (Split *) v1; S2 = (Split *) v2;
305 if(S1->isolation_index < S2->isolation_index) return +1;
306 if(S2->isolation_index < S1->isolation_index) return -1;
310 /* -------------------------------------------------------------------------- */
312 PUBLIC void sort_Split(Split *x)
317 qsort((void *) &(x[1]), x[0].splitsize ,sizeof(Split),CompareSplit );
319 for(i=1; i<=x[0].splitsize; i++) {
320 if(x[i].splitsize > x[i].splitlist[0][0]-x[i].splitsize) {
321 for(j=1; j<=x[i].splitsize; j++) {
322 t1 = x[i].splitlist[0][j];
323 x[i].splitlist[0][j] = x[i].splitlist[1][j];
324 x[i].splitlist[1][j] = t1;
326 x[i].splitsize = x[i].splitlist[0][0]-x[i].splitsize;
332 /* -------------------------------------------------------------------------- */