1 /**********************************************************************
2 ga_similarity.c
3 **********************************************************************
4
5 ga_similarity - Genetic algorithm genome/chromosome comparison routines.
6 Copyright ©2001-2004, Stewart Adcock <stewart@linux-domain.com>
7 All rights reserved.
8
9 The latest version of this program should be available at:
10 http://gaul.sourceforge.net/
11
12 This program is free software; you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation; either version 2 of the License, or
15 (at your option) any later version. Alternatively, if your project
16 is incompatible with the GPL, I will probably agree to requests
17 for permission to use the terms of any other license.
18
19 This program is distributed in the hope that it will be useful, but
20 WITHOUT ANY WARRANTY WHATSOEVER.
21
22 A full copy of the GNU General Public License should be in the file
23 "COPYING" provided with this distribution; if not, see:
24 http://www.gnu.org/
25
26 **********************************************************************
27
28 Synopsis: Routines for comparing genomes/chromosomes.
29
30 Definitions:
31
32 I define a pair of alleles to be matching if they are
33 identical and not matching if the are different - even
34 if the differing values are equivalent (e.g. if the
35 modulus is significant rather than the actual value.)
36
37 It is probably best not to use the following functions,
38 i.e. they should be deprecated and are likely to be
39 deleted soon:
40 boolean ga_compare_genome()
41 int ga_count_match_alleles()
42 double ga_genome_hamming_similarity()
43 double ga_genome_euclidian_similarity()
44
45 References: A general reference that I've found useful in the past
46 but maybe is tricky to obtain, but is:
47 Bradshaw J., "Introduction to Tversky similarity
48 measure", MUG '97 - 11th annual Daylight user group
49 meeting.
50
51 Interested parties are additionally directed toward
52 the documentation:
53 http://gaul.sourceforge.net/similarity.pdf
54
55 To do: Equivalent functions for the other chromosome types.
56
57 **********************************************************************/
58
59 #include "gaul/ga_core.h"
60
61 /**********************************************************************
62 ga_compare_genome()
63 synopsis: Compares two genotypes. Simple memcmp() comparison of
64 the chromosomes.
65 parameters: population *pop Population of entities (you may use
66 differing populations if they are "compatible")
67 entity *alpha Test entity.
68 entity *beta Test entity.
69 return: Returns TRUE if all chromosomes are identical.
70 last updated: 19/12/00
71 **********************************************************************/
72
73 #ifdef COMPILE_DEPRECATED_FUNCTIONS
ga_compare_genome(population * pop,entity * alpha,entity * beta)74 boolean ga_compare_genome(population *pop, entity *alpha, entity *beta)
75 {
76 int i; /* Loop variable over all chromosomes. */
77
78 /* Checks */
79 if (!alpha || !beta) die("Null pointer to entity structure passed");
80
81 for (i=0; i<pop->num_chromosomes; i++)
82 {
83 /*
84 * FIXME: Here is a kludge.
85 * This version is more efficient, but invalid when chromosomes have
86 * differing lengths. When differing length chromosomes are properly handled
87 * this code should be changed.
88 */
89 if ( memcmp( alpha->chromosome[i], beta->chromosome[i], pop->len_chromosomes*sizeof(int) ) !=0 )
90 return FALSE;
91 }
92
93 return TRUE;
94 }
95 #endif
96
97
98 /**********************************************************************
99 ga_count_match_alleles()
100 synopsis: Compares two chromosomes and counts matching alleles.
101 parameters: const int length Chromosome length.
102 const int *alpha Alpha chromosome.
103 const int *beta Beta chromosome.
104 return: Returns number of matching alleles.
105 last updated: 13/02/01
106 **********************************************************************/
107
108 #ifdef COMPILE_DEPRECATED_FUNCTIONS
109
ga_count_match_alleles(const int length,const int * alpha,const int * beta)110 int ga_count_match_alleles(const int length, const int *alpha, const int *beta)
111 {
112 int i; /* Loop variable over all alleles. */
113 int count=0; /* Number of matching alleles. */
114
115 /* Checks. */
116 if (!alpha || !beta) die("Null pointer to chromosome passed");
117
118 for (i=0; i<length; i++)
119 if (alpha[i] == beta[i]) count++;
120
121 return count;
122 }
123 #endif
124
125
126 /**********************************************************************
127 ga_genome_hamming_similarity()
128 synopsis: Measures the Hamming coefficient of the genomes of
129 a pair of entities.
130 parameters: population *pop Population of entities (you may use
131 differing populations if they are "compatible")
132 entity *alpha Test entity.
133 entity *beta Test entity.
134 return: double similarity
135 last updated: 13/02/01
136 **********************************************************************/
137
138 #ifdef COMPILE_DEPRECATED_FUNCTIONS
139
ga_genome_hamming_similarity(population * pop,entity * alpha,entity * beta)140 double ga_genome_hamming_similarity(population *pop, entity *alpha, entity *beta)
141 {
142 int i; /* Loop over chromosomes. */
143 int match=0; /* Count of matching alleles. */
144 int length, total_length=0; /* Number of alleles. */
145
146 /* Checks */
147 if (!alpha || !beta) die("Null pointer to entity structure passed");
148
149 for (i=0; i<pop->num_chromosomes; i++)
150 {
151 /*
152 * FIXME: Here is a kludge.
153 * This version is more efficient, but invalid when chromosomes have
154 * differing lengths. When differing length chromosomes are properly handled
155 * this code should be changed.
156 */
157 length = pop->len_chromosomes;
158
159 total_length += length;
160 match += ga_count_match_alleles(length, alpha->chromosome[i], beta->chromosome[i]);
161 }
162
163 return (double) match/total_length;
164 }
165 #endif
166
167
168 /**********************************************************************
169 ga_genome_euclidian_similarity()
170 synopsis: Measures the Euclidean coefficient of the genomes of
171 a pair of entities.
172 parameters: population *pop Population of entities (you may use
173 differing populations if they are "compatible")
174 entity *alpha Test entity.
175 entity *beta Test entity.
176 return: double similarity
177 last updated: 13/02/01
178 **********************************************************************/
179
180 #ifdef COMPILE_DEPRECATED_FUNCTIONS
181
ga_genome_euclidian_similarity(population * pop,entity * alpha,entity * beta)182 double ga_genome_euclidian_similarity(population *pop, entity *alpha, entity *beta)
183 {
184 return sqrt(ga_genome_hamming_similarity(pop, alpha, beta));
185 }
186 #endif
187
188
189 /**********************************************************************
190 ga_similarity_bitstring_count_1_alleles()
191 synopsis: Counts "1"s in a bitstring chromosome.
192 parameters: const population *pop Population.
193 const entity *alpha entity containing alpha chromosome.
194 const int chromosomeid Index of chromosome to consider.
195 return: Returns number of alleles with value "1".
196 last updated: 16 May 2002
197 **********************************************************************/
198
ga_similarity_bitstring_count_1_alleles(const population * pop,const entity * alpha,const int chromosomeid)199 int ga_similarity_bitstring_count_1_alleles( const population *pop,
200 const entity *alpha, const int chromosomeid )
201 {
202 int i; /* Loop variable over all alleles. */
203 int count=0; /* Number of matching alleles. */
204 byte *a; /* Comparison bitstring. */
205
206 /* Checks. */
207 if (!pop) die("Null pointer to population structure passed");
208 if (!alpha) die("Null pointer to entity structure passed");
209 if (chromosomeid<0 || chromosomeid>=pop->num_chromosomes)
210 die("Invalid chromosome index passed");
211
212 a = (byte*)(alpha->chromosome[chromosomeid]);
213
214 for ( i=0; i<pop->len_chromosomes; i++ )
215 {
216 if (ga_bit_get( a, i ) == 1) count++;
217 }
218
219 return count;
220 }
221
222
223 /**********************************************************************
224 ga_similarity_bitstring_count_match_alleles()
225 synopsis: Compares two bitstring chromosomes and counts matching alleles.
226 parameters: const population *pop Population.
227 const entity *alpha entity containing alpha chromosome.
228 const entity *beta entity containing beta chromosome.
229 const int chromosomeid Index of chromosome to consider.
230 return: Returns number of matching alleles.
231 last updated: 14 May 2002
232 **********************************************************************/
233
ga_similarity_bitstring_count_match_alleles(const population * pop,const entity * alpha,const entity * beta,const int chromosomeid)234 int ga_similarity_bitstring_count_match_alleles( const population *pop,
235 const entity *alpha, const entity *beta,
236 const int chromosomeid )
237 {
238 int i; /* Loop variable over all alleles. */
239 int count=0; /* Number of matching alleles. */
240 byte *a, *b; /* Comparison bitstrings. */
241
242 /* Checks. */
243 if (!pop) die("Null pointer to population structure passed");
244 if (!alpha || !beta) die("Null pointer to entity structure passed");
245 if (chromosomeid<0 || chromosomeid>=pop->num_chromosomes) die("Invalid chromosome index passed");
246
247 a = (byte*)(alpha->chromosome[chromosomeid]);
248 b = (byte*)(beta->chromosome[chromosomeid]);
249
250 for ( i=0; i<pop->len_chromosomes; i++ )
251 {
252 if (ga_bit_get( a, i ) == ga_bit_get( b, i )) count++;
253 }
254
255 return count;
256 }
257
258
259 /**********************************************************************
260 ga_similarity_bitstring_count_and_alleles()
261 synopsis: Compares two bitstring chromosomes and counts alleles
262 which both have value of "1". i.e. Count bits set
263 after AND operation.
264 parameters: const population *pop Population.
265 const entity *alpha entity containing alpha chromosome.
266 const entity *beta entity containing beta chromosome.
267 const int chromosomeid Index of chromosome to consider.
268 return: Returns number of alleles set in both bitstrings.
269 last updated: 14 May 2002
270 **********************************************************************/
271
ga_similarity_bitstring_count_and_alleles(const population * pop,const entity * alpha,const entity * beta,const int chromosomeid)272 int ga_similarity_bitstring_count_and_alleles( const population *pop,
273 const entity *alpha, const entity *beta,
274 const int chromosomeid )
275 {
276 int i; /* Loop variable over all alleles. */
277 int count=0; /* Number of matching alleles. */
278 byte *a, *b; /* Comparison bitstrings. */
279
280 /* Checks. */
281 if (!pop) die("Null pointer to population structure passed");
282 if (!alpha || !beta) die("Null pointer to entity structure passed");
283 if (chromosomeid<0 || chromosomeid>=pop->num_chromosomes) die("Invalid chromosome index passed");
284
285 a = (byte*)(alpha->chromosome[chromosomeid]);
286 b = (byte*)(beta->chromosome[chromosomeid]);
287
288 for ( i=0; i<pop->len_chromosomes; i++ )
289 {
290 if (ga_bit_get( a, i ) && ga_bit_get( b, i )) count++;
291 }
292
293 return count;
294 }
295
296
297 /**********************************************************************
298 ga_similarity_integer_count_match_alleles()
299 synopsis: Compares two chromosomes and counts matching alleles.
300 parameters: const population *pop Population.
301 const entity *alpha entity containing alpha chromosome.
302 const entity *beta entity containing beta chromosome.
303 const int chromosomeid Index of chromosome to consider.
304 return: Returns number of matching alleles.
305 last updated: 14 May 2002
306 **********************************************************************/
307
ga_similarity_integer_count_match_alleles(const population * pop,const entity * alpha,const entity * beta,const int chromosomeid)308 int ga_similarity_integer_count_match_alleles( const population *pop,
309 const entity *alpha, const entity *beta,
310 const int chromosomeid )
311 {
312 int i; /* Loop variable over all alleles. */
313 int count=0; /* Number of matching alleles. */
314 int *a, *b; /* Comparison integer arrays. */
315
316 /* Checks. */
317 if (!pop) die("Null pointer to population structure passed");
318 if (!alpha || !beta) die("Null pointer to entity structure passed");
319 if (chromosomeid<0 || chromosomeid>=pop->num_chromosomes) die("Invalid chromosome index passed");
320
321 a = (int*)(alpha->chromosome[chromosomeid]);
322 b = (int*)(beta->chromosome[chromosomeid]);
323
324 for (i=0; i<pop->len_chromosomes; i++)
325 if (a[i] == b[i]) count++;
326
327 return count;
328 }
329
330
331 /**********************************************************************
332 ga_similarity_bitstring_tanimoto()
333 synopsis: Compares the chromosomes of two entities.
334 parameters: const population *pop Population.
335 const entity *alpha entity containing alpha chromosome.
336 const entity *beta entity containing beta chromosome.
337 return: Returns Tanimoto similarity coefficient.
338 last updated: 16 May 2002
339 **********************************************************************/
340
ga_similarity_bitstring_tanimoto(const population * pop,const entity * alpha,const entity * beta)341 double ga_similarity_bitstring_tanimoto(const population *pop,
342 const entity *alpha, const entity *beta)
343 {
344 int i; /* Loop variable over all chromosomes. */
345 int a=0, b=0; /* Number of ones in the individual entities' chromosomes. */
346 int n=0; /* Number of ones both entities' chromosomes. */
347
348 /* Checks. */
349 if (!pop) die("Null pointer to population structure passed");
350 if (!alpha || !beta) die("Null pointer to entity structure passed");
351
352 for (i=0; i<pop->num_chromosomes; i++)
353 {
354 n += ga_similarity_bitstring_count_and_alleles( pop, alpha, beta, i );
355 a += ga_similarity_bitstring_count_1_alleles( pop, alpha, i );
356 b += ga_similarity_bitstring_count_1_alleles( pop, beta, i );
357 }
358
359 return (double) n/(a+b-n);
360 }
361
362
363 /**********************************************************************
364 ga_similarity_bitstring_dice()
365 synopsis: Compares the chromosomes of two entities.
366 parameters: const population *pop Population.
367 const entity *alpha entity containing alpha chromosome.
368 const entity *beta entity containing beta chromosome.
369 return: Returns Dice similarity coefficient.
370 last updated: 16 May 2002
371 **********************************************************************/
372
ga_similarity_bitstring_dice(const population * pop,const entity * alpha,const entity * beta)373 double ga_similarity_bitstring_dice(const population *pop,
374 const entity *alpha, const entity *beta)
375 {
376 int i; /* Loop variable over all chromosomes. */
377 int a=0, b=0; /* Number of ones in the individual entities' chromosomes. */
378 int n=0; /* Number of ones both entities' chromosomes. */
379
380 /* Checks. */
381 if (!pop) die("Null pointer to population structure passed");
382 if (!alpha || !beta) die("Null pointer to entity structure passed");
383
384 for (i=0; i<pop->num_chromosomes; i++)
385 {
386 n += ga_similarity_bitstring_count_and_alleles( pop, alpha, beta, i );
387 a += ga_similarity_bitstring_count_1_alleles( pop, alpha, i );
388 b += ga_similarity_bitstring_count_1_alleles( pop, beta, i );
389 }
390
391 return (double) 2*n/(a+b);
392 }
393
394
395 /**********************************************************************
396 ga_similarity_bitstring_euclidean()
397 synopsis: Compares the chromosomes of two entities.
398 parameters: const population *pop Population.
399 const entity *alpha entity containing alpha chromosome.
400 const entity *beta entity containing beta chromosome.
401 return: Returns Euclidean similarity coefficient.
402 last updated: 16 May 2002
403 **********************************************************************/
404
ga_similarity_bitstring_euclidean(const population * pop,const entity * alpha,const entity * beta)405 double ga_similarity_bitstring_euclidean(const population *pop,
406 const entity *alpha, const entity *beta)
407 {
408 int i; /* Loop variable over all chromosomes. */
409 int a=0, b=0; /* Number of ones in the individual entities' chromosomes. */
410 int n=0; /* Number of ones both entities' chromosomes. */
411
412 /* Checks. */
413 if (!pop) die("Null pointer to population structure passed");
414 if (!alpha || !beta) die("Null pointer to entity structure passed");
415
416 for (i=0; i<pop->num_chromosomes; i++)
417 {
418 n += ga_similarity_bitstring_count_and_alleles( pop, alpha, beta, i );
419 a += ga_similarity_bitstring_count_1_alleles( pop, alpha, i );
420 b += ga_similarity_bitstring_count_1_alleles( pop, beta, i );
421 }
422
423 return 1.0 - sqrt((double)(a+b-2*n))/pop->len_chromosomes;
424 }
425
426
427 /**********************************************************************
428 ga_similarity_bitstring_hamming()
429 synopsis: Compares the chromosomes of two entities.
430 parameters: const population *pop Population.
431 const entity *alpha entity containing alpha chromosome.
432 const entity *beta entity containing beta chromosome.
433 return: Returns Hamming similarity coefficient.
434 last updated: 16 May 2002
435 **********************************************************************/
436
ga_similarity_bitstring_hamming(const population * pop,const entity * alpha,const entity * beta)437 double ga_similarity_bitstring_hamming(const population *pop,
438 const entity *alpha, const entity *beta)
439 {
440 int i; /* Loop variable over all chromosomes. */
441 int a=0, b=0; /* Number of ones in the individual entities' chromosomes. */
442 int n=0; /* Number of ones both entities' chromosomes. */
443
444 /* Checks. */
445 if (!pop) die("Null pointer to population structure passed");
446 if (!alpha || !beta) die("Null pointer to entity structure passed");
447
448 for (i=0; i<pop->num_chromosomes; i++)
449 {
450 n += ga_similarity_bitstring_count_and_alleles( pop, alpha, beta, i );
451 a += ga_similarity_bitstring_count_1_alleles( pop, alpha, i );
452 b += ga_similarity_bitstring_count_1_alleles( pop, beta, i );
453 }
454
455 return 1.0 - (a+b-2*n)/pop->len_chromosomes;
456 }
457
458
459 /**********************************************************************
460 ga_similarity_bitstring_cosine()
461 synopsis: Compares the chromosomes of two entities.
462 parameters: const population *pop Population.
463 const entity *alpha entity containing alpha chromosome.
464 const entity *beta entity containing beta chromosome.
465 return: Returns Cosine similarity coefficient.
466 last updated: 16 May 2002
467 **********************************************************************/
468
ga_similarity_bitstring_cosine(const population * pop,const entity * alpha,const entity * beta)469 double ga_similarity_bitstring_cosine(const population *pop,
470 const entity *alpha, const entity *beta)
471 {
472 int i; /* Loop variable over all chromosomes. */
473 int a=0, b=0; /* Number of ones in the individual entities' chromosomes. */
474 int n=0; /* Number of ones both entities' chromosomes. */
475
476 /* Checks. */
477 if (!pop) die("Null pointer to population structure passed");
478 if (!alpha || !beta) die("Null pointer to entity structure passed");
479
480 for (i=0; i<pop->num_chromosomes; i++)
481 {
482 n += ga_similarity_bitstring_count_and_alleles( pop, alpha, beta, i );
483 a += ga_similarity_bitstring_count_1_alleles( pop, alpha, i );
484 b += ga_similarity_bitstring_count_1_alleles( pop, beta, i );
485 }
486
487 if (a==0 || b==0) return 0;
488
489 return n/sqrt((double)(a*b));
490 }
491
492
493 /**********************************************************************
494 ga_similarity_double_count_match_alleles()
495 synopsis: Compares two "double" chromosomes and counts matching
496 alleles. A match is defined to be when
497 x+GA_TINY_DOUBLE>y>x-GA_TINY_DOUBLE
498 parameters: const population *pop Population.
499 const entity *alpha entity containing alpha chromosome.
500 const entity *beta entity containing beta chromosome.
501 const int chromosomeid Index of chromosome to consider.
502 return: Returns number of matching alleles.
503 last updated: 23 May 2002
504 **********************************************************************/
505
ga_similarity_double_count_match_alleles(const population * pop,const entity * alpha,const entity * beta,const int chromosomeid)506 int ga_similarity_double_count_match_alleles( const population *pop,
507 const entity *alpha, const entity *beta,
508 const int chromosomeid )
509 {
510 int i; /* Loop variable over all alleles. */
511 int count=0; /* Number of matching alleles. */
512 double *a, *b; /* Comparison chromosomes. */
513
514 /* Checks. */
515 if (!pop) die("Null pointer to population structure passed");
516 if (!alpha || !beta) die("Null pointer to entity structure passed");
517 if (chromosomeid<0 || chromosomeid>=pop->num_chromosomes) die("Invalid chromosome index passed");
518
519 a = (double*)(alpha->chromosome[chromosomeid]);
520 b = (double*)(beta->chromosome[chromosomeid]);
521
522 for ( i=0; i<pop->len_chromosomes; i++ )
523 {
524 if (a[i]+GA_TINY_DOUBLE>b[i] && b[i]>a[i]-GA_TINY_DOUBLE) count++;
525 }
526
527 return count;
528 }
529
530
531 /**********************************************************************
532 ga_similarity_double_tanimoto()
533 synopsis: Compares the chromosomes of two entities.
534 parameters: const population *pop Population.
535 const entity *alpha entity containing alpha chromosome.
536 const entity *beta entity containing beta chromosome.
537 return: Returns Tanimoto similarity coefficient.
538 Range: -1/3 to +1
539 last updated: 01 Apr 2004
540 **********************************************************************/
541
ga_similarity_double_tanimoto(const population * pop,const entity * alpha,const entity * beta)542 double ga_similarity_double_tanimoto(const population *pop,
543 const entity *alpha, const entity *beta)
544 {
545 int i, j; /* Loop variable over all chromosomes, alleles. */
546 double ab=0.0, aa=0.0, bb=0.0; /* Components of the similarity equation. */
547 double *a, *b; /* Comparison chromosomes. */
548
549 /* Checks. */
550 if (!pop) die("Null pointer to population structure passed");
551 if (!alpha || !beta) die("Null pointer to entity structure passed");
552
553 for (i=0; i<pop->num_chromosomes; i++)
554 {
555 a = (double*)(alpha->chromosome[i]);
556 b = (double*)(beta->chromosome[i]);
557
558 for (j=0; j<pop->len_chromosomes; j++)
559 {
560 aa += a[j]*a[j];
561 ab += a[j]*b[j];
562 bb += b[j]*b[j];
563 }
564 }
565
566 return ab/(aa+bb-ab);
567 }
568
569
570 /**********************************************************************
571 ga_similarity_double_dice()
572 synopsis: Compares the chromosomes of two entities.
573 parameters: const population *pop Population.
574 const entity *alpha entity containing alpha chromosome.
575 const entity *beta entity containing beta chromosome.
576 return: Returns Dice similarity coefficient.
577 Range: -1 to +1
578 last updated: 01 Apr 2004
579 **********************************************************************/
580
ga_similarity_double_dice(const population * pop,const entity * alpha,const entity * beta)581 double ga_similarity_double_dice(const population *pop,
582 const entity *alpha, const entity *beta)
583 {
584 int i, j; /* Loop variable over all chromosomes, alleles. */
585 double ab=0.0, aa=0.0, bb=0.0; /* Components of the similarity equation. */
586 double *a, *b; /* Comparison chromosomes. */
587
588 /* Checks. */
589 if (!pop) die("Null pointer to population structure passed");
590 if (!alpha || !beta) die("Null pointer to entity structure passed");
591
592 for (i=0; i<pop->num_chromosomes; i++)
593 {
594 a = (double*)(alpha->chromosome[i]);
595 b = (double*)(beta->chromosome[i]);
596
597 for (j=0; j<pop->len_chromosomes; j++)
598 {
599 aa += a[j]*a[j];
600 ab += a[j]*b[j];
601 bb += b[j]*b[j];
602 }
603 }
604
605 return (2*ab)/(aa+bb);
606 }
607
608
609 /**********************************************************************
610 ga_similarity_double_cosine()
611 synopsis: Compares the chromosomes of two entities.
612 parameters: const population *pop Population.
613 const entity *alpha entity containing alpha chromosome.
614 const entity *beta entity containing beta chromosome.
615 return: Returns Cosine similarity coefficient.
616 Range: -1 to +1
617 last updated: 01 Apr 2004
618 **********************************************************************/
619
ga_similarity_double_cosine(const population * pop,const entity * alpha,const entity * beta)620 double ga_similarity_double_cosine(const population *pop,
621 const entity *alpha, const entity *beta)
622 {
623 int i, j; /* Loop variable over all chromosomes, alleles. */
624 double ab=0.0, aa=0.0, bb=0.0; /* Components of the similarity equation. */
625 double *a, *b; /* Comparison chromosomes. */
626
627 /* Checks. */
628 if (!pop) die("Null pointer to population structure passed");
629 if (!alpha || !beta) die("Null pointer to entity structure passed");
630
631 for (i=0; i<pop->num_chromosomes; i++)
632 {
633 a = (double*)(alpha->chromosome[i]);
634 b = (double*)(beta->chromosome[i]);
635
636 for (j=0; j<pop->len_chromosomes; j++)
637 {
638 aa += a[j]*a[j];
639 ab += a[j]*b[j];
640 bb += b[j]*b[j];
641 }
642 }
643
644 return ab/sqrt(aa+bb);
645 }
646
647
648