ARB
ali_profile.cxx
Go to the documentation of this file.
1 // =============================================================== //
2 // //
3 // File : ali_profile.cxx //
4 // Purpose : //
5 // //
6 // Institute of Microbiology (Technical University Munich) //
7 // http://www.arb-home.de/ //
8 // //
9 // =============================================================== //
10 
11 #include "ali_profile.hxx"
12 
13 #include <BI_helix.hxx>
14 #include <cctype>
15 
16 
17 inline int ALI_PROFILE::is_binding_marker(char c) {
18  return c == '~' || c == 'x';
19 }
20 
21 
22 ALI_TLIST<ali_family_member *> *ALI_PROFILE::find_family(ALI_SEQUENCE *Sequence, ALI_PROFILE_CONTEXT *context) {
23  // find the family in the pt server
24  char message_buffer[100];
25  ALI_PT &pt = (ALI_PT &) *(context->pt);
26  ALI_ARBDB &arbdb = (ALI_ARBDB &) *(context->arbdb);
28  ali_family_member *family_member;
29  ALI_TLIST<ali_family_member *> *family_list;
30  ALI_TLIST<ali_pt_member *> *pt_fam_list;
31  ALI_TLIST<ali_pt_member *> *pt_ext_list;
32  ali_pt_member *pt_member;
33  float weight, d;
34  unsigned long number;
35 
36  // Initialization
37  family_list = new ALI_TLIST<ali_family_member *>;
38 
39  ali_message("Searching for the family");
40  pt.find_family(Sequence, context->find_family_mode);
41  ali_message("Family found");
42 
43  pt_fam_list = pt.get_family_list();
44  pt_ext_list = pt.get_extension_list();
45 
46  ali_message("Reading family:");
47 
48  d = (context->ext_max_weight - 1.0) / (float) pt_fam_list->cardinality();
49 
50  arbdb.begin_transaction();
51 
52  // calculate the real family members
53  number = 0;
54  while (!pt_fam_list->is_empty()) {
55  pt_member = pt_fam_list->first();
56  seq = arbdb.get_sequence(pt_member->name, context->mark_family_flag);
57  if (seq) {
58  weight = 1 + d * number;
59  sprintf(message_buffer, "%s (weight = %f, matches = %d)",
60  pt_member->name, weight, pt_member->matches);
61  ali_message(message_buffer);
62  family_member = new ali_family_member(seq,
63  (float) pt_member->matches,
64  weight);
65  family_list->append_end(family_member);
66  number++;
67  }
68  else {
69  ali_warning("Sequence not found in Database (Sequence ignored)");
70  }
71  pt_fam_list->delete_element();
72  }
73  delete pt_fam_list;
74 
75  ali_message("Reading family extension:");
76 
77  d = -1.0 * context->ext_max_weight / (float) pt_ext_list->cardinality();
78 
79  // calculate the extension of the family
80  number = 0;
81  while (!pt_ext_list->is_empty()) {
82  pt_member = pt_ext_list->first();
83  seq = arbdb.get_sequence(pt_member->name,
85  if (seq) {
86  weight = context->ext_max_weight + d * number;
87  sprintf(message_buffer, "%s (weight = %f, matches = %d)",
88  pt_member->name, weight, pt_member->matches);
89  ali_message(message_buffer);
90  family_member = new ali_family_member(seq,
91  (float) pt_member->matches,
92  weight);
93  family_list->append_end(family_member);
94  number++;
95  }
96  else {
97  ali_warning("Sequence not found in Database (Sequence ignored)");
98  }
99  pt_ext_list->delete_element();
100  }
101  delete pt_ext_list;
102 
103  arbdb.commit_transaction();
104 
105  return family_list;
106 }
107 
108 void ALI_PROFILE::calculate_costs(ALI_TLIST<ali_family_member *> *family_list, ALI_PROFILE_CONTEXT *context) {
109  // calculate the costs for aligning against a family
110  ali_family_member *family_member;
111  float a[7], w[7], w_sum, sm[7][7];
112  float base_gap_weights[5], w_bg_sum;
113  long members;
114  size_t p;
115  int i;
116  size_t j;
117  unsigned long *l;
118  float *g;
119  unsigned char **seq;
120  long *seq_len;
121  float (*w_Del)[], (*percent)[];
122 
123  // allocate temporary memory
124  members = family_list->cardinality();
125  l = (unsigned long *) CALLOC((unsigned int) members, sizeof(long));
126  g = (float *) CALLOC((unsigned int) members, sizeof(float));
127  seq = (unsigned char **) CALLOC((unsigned int) members, sizeof(char *));
128  seq_len = (long *) CALLOC((unsigned int) members, sizeof(long));
129  ali_out_of_memory_if(!l || !g || !seq || !seq_len);
130 
131  // initialize arrays
132  family_member = family_list->first();
133  prof_len = family_member->sequence->length();
134  seq[0] = family_member->sequence->sequence();
135  seq_len[0] = family_member->sequence->length();
136  g[0] = family_member->weight;
137  i = 1;
138  sub_costs_maximum = 0.0;
139 
140  while (family_list->has_next()) {
141  family_member = family_list->next();
142  seq[i] = family_member->sequence->sequence();
143  seq_len[i] = family_member->sequence->length();
144  g[i] = family_member->weight;
145  i++;
146  if (prof_len < family_member->sequence->length()) {
147  ali_warning("Family members have different length");
148  prof_len = family_member->sequence->length();
149  }
150  }
151 
152  // Calculate the substitution cost matrix
153  for (i = 0; i < 5; i++)
154  for (j = 0; j < 5; j++)
155  sm[i][j] = context->substitute_matrix[i][j];
156 
157  // Initialize l-array (position of last base)
158  for (i = 0; i < members; i++)
159  l[i] = prof_len + 1;
160 
161  // allocate memory for costs
162 
163  base_weights = (float (**) [4]) CALLOC((unsigned int) prof_len, sizeof(float [4]));
164  sub_costs = (float (**) [6]) CALLOC((unsigned int) prof_len, sizeof(float [6]));
165  binding_costs = (float (*) [5][5]) CALLOC((unsigned int) 5, sizeof(float [5]));
166  lmin = (unsigned long *) CALLOC((unsigned int) prof_len, sizeof(unsigned long));
167  lmax = (unsigned long *) CALLOC((unsigned int) prof_len, sizeof(unsigned long));
168  gap_costs = (float ***) CALLOC((unsigned int) prof_len, sizeof(float *));
169  gap_percents = (float***) CALLOC((unsigned int) prof_len, sizeof(float *));
170 
171  ali_out_of_memory_if(!binding_costs || !sub_costs || !lmin || !lmax || !gap_costs || !gap_percents || !base_weights);
172 
173  // Copy the binding costs matrix
174  w_bind_maximum = context->binding_matrix[0][0];
175  for (i = 0; i < 5; i++)
176  for (j = 0; j < 5; j++) {
177  (*binding_costs)[i][j] = context->binding_matrix[i][j];
178  if ((*binding_costs)[i][j] > w_bind_maximum)
179  w_bind_maximum = (*binding_costs)[i][j];
180  }
181 
182  // calculate the costs for EVERY position
183  ali_message("Calculating costs for substitution");
184  for (p = 0; p < prof_len; p++) {
185  // Initialization
186  for (i = 0; i < 7; i++)
187  a[i] = w[i] = sm[5][i] = sm[i][5] = sm[6][i] = sm[i][6] = 0.0;
188  for (i = 0; i < 6; i++)
189  (*sub_costs)[p][i] = 0.0;
190  w_sum = 0.0;
191  w_bg_sum = 0.0;
192 
193  // Statistic consensus
194  for (i = 0; i < members; i++) {
195  if (p < size_t(seq_len[i])) {
196  a[seq[i][p]]++;
197  w[seq[i][p]] += g[i];
198  if (ali_is_real_base(seq[i][p]))
199  w_sum += g[i];
200  if (ali_is_real_base(seq[i][p]) || ali_is_gap(seq[i][p]))
201  w_bg_sum += g[i];
202  }
203  else {
204  a[ALI_DOT_CODE]++;
205  w[ALI_DOT_CODE] += g[i];
206  }
207  }
208 
209  // Relative weight of bases
210  if (w_sum != 0)
211  for (i = 0; i < 4; i++)
212  (*base_weights)[p][i] = w[i] / w_sum;
213  else
214  for (i = 0; i < 4; i++)
215  (*base_weights)[p][i] = 0.25;
216 
217  // Relative weight of bases and gaps
218  if (w_bg_sum != 0)
219  for (i = 0; i < 5; i++)
220  base_gap_weights[i] = w[i] / w_bg_sum;
221  else
222  for (i = 0; i < 5; i++)
223  base_gap_weights[i] = 0.2;
224 
225  // Expandation of substitute matrix (for 'n')
226  for (j = 0; j < 5; j++) {
227  for (i = 0; i < 4; i++) {
228  sm[5][j] += (*base_weights)[p][i] * sm[i][j];
229  sm[j][5] += (*base_weights)[p][i] * sm[j][i];
230  }
231  }
232  for (i = 0; i < 4; i++)
233  sm[5][5] += (*base_weights)[p][i] * sm[i][i];
234 
235  // Expandation of substitute matrix (for '.')
236  for (j = 0; j < 6; j++)
237  for (i = 0; i < 5; i++) {
238  sm[6][j] += base_gap_weights[i] * sm[i][j];
239  sm[j][6] += base_gap_weights[i] * sm[j][i];
240  }
241  for (i = 0; i < 5; i++)
242  sm[6][6] += base_gap_weights[i] * sm[i][i];
243 
244  // Substitution costs
245  for (i = 0; i < members; i++) {
246  if (p < size_t(seq_len[i])) {
247  for (j = 0; j < 6; j++) {
248  (*sub_costs)[p][j] += g[i] * sm[seq[i][p]][j];
249  }
250  }
251  else {
252  for (j = 0; j < 6; j++) {
253  (*sub_costs)[p][j] += g[i] * sm[ALI_DOT_CODE][j];
254  }
255  }
256  }
257  for (j = 0; j < 6; j++) {
258  (*sub_costs)[p][j] /= members;
259  if ((*sub_costs)[p][j] > sub_costs_maximum)
260  sub_costs_maximum = (*sub_costs)[p][j];
261  }
262 
263  // Calculate dynamic deletion costs and percents of real gaps
264  lmax[p] = 0;
265  lmin[p] = p;
266  for (i = 0; i < members; i++)
267  if (l[i] < p) {
268  if (lmin[p] > l[i])
269  lmin[p] = l[i];
270  if (lmax[p] < l[i])
271  lmax[p] = l[i];
272  }
273  if (lmin[p] == p && lmax[p] == 0) {
274  lmin[p] = lmax[p] = p;
275  }
276 
277  w_Del = (float (*) []) CALLOC((unsigned int) (lmax[p]-lmin[p]+2), sizeof(float));
278  percent = (float (*) []) CALLOC((unsigned int) (lmax[p]-lmin[p]+2), sizeof(float));
279  ali_out_of_memory_if(!w_Del || !percent);
280  (*gap_costs)[p] = (float *) w_Del;
281  (*gap_percents)[p] = (float *) percent;
282 
283  // Calculate dynamic deletion costs
284  for (j = 0; j <= lmax[p] - lmin[p] + 1; j++) {
285  (*w_Del)[j] = 0;
286  for (i = 0; i < members; i++) {
287  // Normal case
288  if (p < size_t(seq_len[i])) {
289  if (l[i] == prof_len + 1 || l[i] >= j + lmin[p]) {
290  (*w_Del)[j] += g[i] * sm[seq[i][p]][4] * context->multi_gap_factor;
291  }
292  else {
293  (*w_Del)[j] += g[i] * sm[seq[i][p]][4];
294  }
295  }
296  // expand sequence with dots
297  else {
298  if (l[i] >= j + lmin[p] && l[i] < prof_len+1) {
299  (*w_Del)[j] += g[i] * sm[ALI_DOT_CODE][4] * context->multi_gap_factor;
300  }
301  else {
302  (*w_Del)[j] += g[i] * sm[ALI_DOT_CODE][4];
303  }
304  }
305  }
306  (*w_Del)[j] /= members;
307  }
308 
309  // Update the l-array
310  for (i = 0; i < members; i++) {
311  if (!ali_is_gap(seq[i][p]))
312  l[i] = p;
313  }
314 
315  // Calculate percents of real gaps
316  for (j = 0; j <= lmax[p] - lmin[p] + 1; j++) {
317  (*percent)[j] = 0;
318  for (i = 0; i < members; i++) {
319  if (l[i] >= j + lmin[p] && l[i] < prof_len+1) {
320  (*percent)[j] += 1.0;
321  }
322  }
323  (*percent)[j] /= members;
324  }
325  }
326 
327  ali_message("Calculation finished");
328 
329  free(l);
330  free(g);
331  free(seq);
332  free(seq_len);
333 }
334 
335 int ALI_PROFILE::find_next_helix(char h[], unsigned long h_len,
336  unsigned long pos,
337  unsigned long *helix_nr,
338  unsigned long *start, unsigned long *end)
339 {
340  // find the next helix
341  unsigned long i;
342 
343  for (i = pos; i < h_len && !isdigit(h[i]); i++) ;
344  if (i >= h_len) return -1;
345 
346  *start = i;
347  sscanf(&h[i], "%ld", helix_nr);
348  for (; i < h_len && isdigit(h[i]); i++) ;
349  for (; i < h_len && !isdigit(h[i]); i++) ;
350  *end = i - 1;
351 
352  return 0;
353 }
354 
355 int ALI_PROFILE::find_comp_helix(char h[], unsigned long h_len,
356  unsigned long pos, unsigned long helix_nr,
357  unsigned long *start, unsigned long *end)
358 {
359  // find the complementary part of a helix
360  unsigned long nr, i;
361 
362  i = pos;
363  do {
364  for (; i < h_len && !isdigit(h[i]); i++) ;
365  if (i >= h_len) return -1;
366  *start = i;
367  sscanf(&h[i], "%ld", &nr);
368  for (; i < h_len && isdigit(h[i]); i++) ;
369  } while (helix_nr != nr);
370 
371  for (; i < h_len && !isdigit(h[i]); i++) ;
372  *end = i - 1;
373 
374  return 0;
375 }
376 
377 void ALI_PROFILE::delete_comp_helix(char h1[], char h2[], unsigned long h_len,
378  unsigned long start, unsigned long end)
379 {
380  unsigned long i;
381 
382  for (i = start; i < h_len && i <= end; i++) {
383  h1[i] = '.';
384  h2[i] = '.';
385  }
386 }
387 
388 
389 void ALI_PROFILE::initialize_helix(ALI_PROFILE_CONTEXT *context) {
390  // initialize the array, representing the helix
391  const char *error_string;
392  BI_helix bi_helix;
393 
394  unsigned long i;
395 
396  // read helix
397  if ((error_string = bi_helix.init(context->arbdb->gb_main)))
398  ali_warning(error_string);
399 
400  helix_len = bi_helix.size();
401  helix = (long **) CALLOC((unsigned int) helix_len, sizeof(long));
402  helix_borders = (char **) CALLOC((unsigned int) helix_len, sizeof(long));
403  ali_out_of_memory_if(!helix || !helix_borders);
404 
405  // convert helix for internal use
406  for (i = 0; i < helix_len; i++)
407  if (bi_helix.pairtype(i) == HELIX_PAIR)
408  (*helix)[i] = bi_helix.opposite_position(i);
409  else
410  (*helix)[i] = -1;
411 }
412 
413 
415  char message_buffer[100];
416  ali_family_member *family_member;
417  ALI_TLIST<ali_family_member *> *family_list;
418 
419  norm_sequence = new ALI_NORM_SEQUENCE(seq);
420 
421  multi_gap_factor = context->multi_gap_factor;
422 
423  initialize_helix(context);
424 
425  family_list = find_family(seq, context);
426  if (family_list->is_empty()) {
427  ali_error("Family not found (maybe incompatible PT and DB Servers)");
428  }
429 
430  calculate_costs(family_list, context);
431 
432  insert_cost = sub_costs_maximum * context->insert_factor;
433  multi_insert_cost = insert_cost * context->multi_insert_factor;
434 
435  sprintf(message_buffer, "Multi gap factor = %f", multi_gap_factor);
436  ali_message(message_buffer);
437  sprintf(message_buffer, "Maximal substitution costs = %f", sub_costs_maximum);
438  ali_message(message_buffer);
439  sprintf(message_buffer, "Normal insert costs = %f", insert_cost);
440  ali_message(message_buffer);
441  sprintf(message_buffer, "Multiple insert costs = %f", multi_insert_cost);
442  ali_message(message_buffer);
443 
444  // Delete the family list
445  family_member = family_list->first();
446  delete family_member->sequence;
447  delete family_member;
448  while (family_list->has_next()) {
449  family_member = family_list->next();
450  delete family_member->sequence;
451  delete family_member;
452  }
453  delete family_list;
454 }
455 
457  size_t i;
458 
459  free(helix);
460  free(helix_borders);
461  free(binding_costs);
462  free(sub_costs);
463  if (gap_costs) {
464  for (i = 0; i < prof_len; i++)
465  if ((*gap_costs)[i])
466  free((*gap_costs)[i]);
467  free(gap_costs);
468  }
469  if (gap_percents) {
470  for (i = 0; i < prof_len; i++)
471  if ((*gap_percents)[i])
472  free((*gap_percents)[i]);
473  free(gap_percents);
474  }
475  free(lmin);
476  free(lmax);
477  delete norm_sequence;
478 }
479 
480 
481 int ALI_PROFILE::is_in_helix(unsigned long pos, unsigned long *first, unsigned long *last) {
482  // test whether a position is inside a helix
483  long i;
484 
485  if (pos > helix_len)
486  return 0;
487 
488  switch ((*helix_borders)[pos]) {
490  *first = pos;
491  for (i = (long) pos + 1; i < (long) prof_len; i++)
492  if ((*helix_borders)[i] == ALI_PROFILE_BORDER_RIGHT) {
493  *last = (unsigned long) i;
494  return 1;
495  }
496  ali_warning("Helix borders inconsistent (1)");
497  return 0;
499  *last = pos;
500  for (i = (long) pos - 1; i >= 0; i--)
501  if ((*helix_borders)[i] == ALI_PROFILE_BORDER_LEFT) {
502  *first = (unsigned long) i;
503  return 1;
504  }
505  ali_warning("Helix borders inconsistent (2)");
506  return 0;
507  default:
508  for (i = (long) pos - 1; i >= 0; i--)
509  switch ((*helix_borders)[i]) {
511  return 0;
513  *first = (unsigned long) i;
514  for (i = (long) pos + 1; i < (long) prof_len; i++)
515  switch ((*helix_borders)[i]) {
517  ali_warning("Helix borders inconsistent (3)");
518  printf("pos = %ld\n", pos);
519  return 0;
521  *last = (unsigned long) i;
522  return 1;
523  }
524  }
525  }
526  return 0;
527 }
528 
529 int ALI_PROFILE::is_outside_helix(unsigned long pos, unsigned long *first, unsigned long *last) {
530  // test, whether a position is outside a helix
531  long i;
532 
533  switch ((*helix_borders)[pos]) {
535  return 0;
537  return 0;
538  default:
539  for (i = (long) pos - 1; i >= 0; i--)
540  switch ((*helix_borders)[i]) {
542  return 0;
544  *first = (unsigned long) i + 1;
545  for (i = (long) pos + 1; i < (long) prof_len; i++)
546  switch ((*helix_borders)[i]) {
548  *last = (unsigned long) i - 1;
549  return 1;
551  ali_warning("Helix borders inconsistent [1]");
552  return 0;
553  }
554  *last = prof_len - 1;
555  return 1;
556  }
557  *first = 0;
558  for (i = (long) pos + 1; i < (long) prof_len; i++)
559  switch ((*helix_borders)[i]) {
561  *last = (unsigned long) i - 1;
562  return 1;
564  ali_warning("Helix borders inconsistent [2]");
565  return 0;
566  }
567  *last = prof_len - 1;
568  return 1;
569  }
570 }
571 
572 
574  // generate a 'consensus sequence'
575 
576  char *seq;
577  size_t p;
578  int i, min_i;
579  float min;
580 
581 
582  seq = (char *) CALLOC((unsigned int) prof_len + 1, sizeof(char));
583  ali_out_of_memory_if(!seq);
584 
585  for (p = 0; p < prof_len; p++) {
586  min = (*sub_costs)[p][0];
587  min_i = 0;
588  for (i = 1; i < 5; i++) {
589  if (min > (*sub_costs)[p][i]) {
590  min = (*sub_costs)[p][i];
591  min_i = i;
592  }
593  else {
594  if (min == (*sub_costs)[p][i])
595  min_i = -1;
596  }
597  }
598  if (min_i >= 0)
599  seq[p] = ali_number_to_base(min_i);
600  else {
601  if (min > 0)
602  seq[p] = '*';
603  else
604  seq[p] = '.';
605  }
606  }
607  seq[prof_len] = '\0';
608 
609  return seq;
610 }
611 
612 float ALI_PROFILE::w_binding(unsigned long first_seq_pos, ALI_SEQUENCE *seq) {
613  // calculate the costs of a binding
614  unsigned long pos_1_seq, pos_2_seq, last_seq_pos;
615  long pos_compl;
616  float costs = 0;
617 
618  last_seq_pos = first_seq_pos + seq->length() - 1;
619  for (pos_1_seq = first_seq_pos; pos_1_seq <= last_seq_pos; pos_1_seq++) {
620  pos_compl = (*helix)[pos_1_seq];
621  if (pos_compl >= 0) {
622  pos_2_seq = (unsigned long) pos_compl;
623  if (pos_2_seq > pos_1_seq && pos_2_seq <= last_seq_pos)
624  costs += w_bind(pos_1_seq, seq->base(pos_1_seq),
625  pos_2_seq, seq->base(pos_2_seq));
626  else
627  if (pos_2_seq < first_seq_pos || pos_2_seq > last_seq_pos)
628  costs += w_bind_maximum;
629  }
630  }
631 
632  return costs;
633 }
634 
635 
float w_bind(unsigned long pos_1, unsigned char base_1, unsigned long pos_2, unsigned char base_2)
int ali_is_real_base(char c)
Definition: ali_misc.hxx:76
void ali_out_of_memory_if(bool cond)
Definition: ali_misc.hxx:52
unsigned char base(unsigned long position)
ALI_PROFILE(ALI_SEQUENCE *sequence, ALI_PROFILE_CONTEXT *context)
int has_next()
Definition: ali_tlist.hxx:168
long
Definition: AW_awar.cxx:154
char ali_number_to_base(unsigned char n)
Definition: ali_misc.hxx:134
int find_family(ALI_SEQUENCE *sequence, int find_type=1)
Definition: ali_pt.cxx:177
float w_binding(unsigned long first_position_of_sequence, ALI_SEQUENCE *sequence)
int is_in_helix(unsigned long pos, unsigned long *first, unsigned long *last)
size_t opposite_position(size_t pos) const
Definition: BI_helix.hxx:96
int cardinality()
Definition: ali_tlist.hxx:162
void delete_element()
Definition: ali_tlist.hxx:355
unsigned char * sequence()
#define ALI_PROFILE_BORDER_RIGHT
Definition: ali_profile.hxx:22
FILE * seq
Definition: rns.c:46
static HelixNrInfo * start
char * cheapest_sequence()
#define ALI_DOT_CODE
Definition: ali_misc.hxx:40
int ali_is_gap(char c)
Definition: ali_misc.hxx:112
void * CALLOC(long i, long j)
Definition: ali_misc.hxx:56
static int weight[maxsites+1]
ALI_NORM_SEQUENCE * sequence()
GBDATA * gb_main
Definition: ali_arbdb.hxx:24
ALI_SEQUENCE * sequence
Definition: ali_profile.hxx:55
char * name
Definition: ali_pt.hxx:43
ALI_SEQUENCE * get_sequence(char *name, int and_mark=0)
Definition: ali_arbdb.cxx:69
void ali_error(const char *message, const char *func)
Definition: ali_main.cxx:90
ALI_TLIST< ali_pt_member * > * get_family_list()
Definition: ali_pt.cxx:302
int is_empty()
Definition: ali_tlist.hxx:165
void begin_transaction()
Definition: ali_arbdb.hxx:35
size_t size() const
Definition: BI_helix.hxx:88
unsigned long length()
ALI_TLIST< ali_pt_member * > * get_extension_list()
Definition: ali_pt.cxx:312
BI_PAIR_TYPE pairtype(size_t pos) const
Definition: BI_helix.hxx:101
void ali_warning(const char *message)
Definition: ali_misc.hxx:47
void ali_message(const char *message)
Definition: ali_misc.hxx:46
void commit_transaction()
Definition: ali_arbdb.hxx:38
int is_outside_helix(unsigned long pos, unsigned long *first, unsigned long *last)
void append_end(T &a)
Definition: ali_tlist.hxx:198
#define min(a, b)
Definition: f2c.h:153
double binding_matrix[5][5]
Definition: ali_profile.hxx:49
unsigned long length()
double substitute_matrix[5][5]
Definition: ali_profile.hxx:48
#define ALI_PROFILE_BORDER_LEFT
Definition: ali_profile.hxx:21