sources/sov/sov.c

   1 /*-----------------------------------------------------------
   2 /
   3 /   Program:      sov.c
   4 /
   5 /   Secondary structure prediction accuracy evaluation
   6 /
   7 /   SOV (Segment OVerlap) measure
   8 /
   9 /   Copyright by Adam Zemla (11/16/1996)
  10 /   Email: adamz@llnl.gov
  11 /
  12 /------------------------------------------------------------
  13 /
  14 /   Compile:      cc sov.c -o sov -lm
  15 /
  16 /------------------------------------------------------------*/
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <string.h>
  20 #include <math.h>
  21
  22 #define MAXRES            2000
  23
  24 typedef struct
  25 {
  26   int input;
  27   int order;
  28   int q3_what;
  29   int sov_what;
  30   int sov_method;
  31   float sov_delta;
  32   float sov_delta_s;
  33   int sov_out;
  34   char fname[80];
  35 } parameters;
  36
  37 char *letter_AA = "ARNDCQEGHILKMFPSTWYV-?X";    /* 23 chars */
  38
  39 void default_parameters (parameters *);
  40 int read_aa_osec_psec (char[MAXRES], char[MAXRES], char[MAXRES], parameters *, char *);
  41 float sov (int, char[MAXRES], char[MAXRES], parameters *);
  42 float q3 (int, char[MAXRES], char[MAXRES], parameters *);
  43 int check_aa (char, char *, int);
  44
  45 int
  46 main (int argc, char *argv[])
  47 {
  48   int i, n_aa, sov_method;
  49   char c, aa[MAXRES], osec[MAXRES], psec[MAXRES];
  50   parameters pdata;
  51   float out0, out1, out2, out3;
  52
  53   if (argc < 2) {
  54     printf (" Usage: sov <input_data>\n");
  55     printf (" HELP:  sov -h\n");
  56     exit (0);
  57   }
  58   if (!strncmp (argv[1], "-h\0", 2) || !strncmp (argv[1], "help\0", 5) || !strncmp (argv[1], "-help\0", 6)) {
  59     system ("more ./README.sov");
  60     printf ("\n");
  61     exit (0);
  62   }
  63
  64   default_parameters (&pdata);
  65
  66   strcpy (pdata.fname, argv[1]);
  67
  68   n_aa = read_aa_osec_psec (aa, osec, psec, &pdata, letter_AA);
  69
  70   if (pdata.input == 1) {
  71     n_aa = read_aa_osec_psec (aa, osec, psec, &pdata, letter_AA);
  72   }
  73
  74   if (pdata.order == 1) {
  75     for (i = 0; i < n_aa; i++) {
  76       c = osec[i];
  77       osec[i] = psec[i];
  78       psec[i] = c;
  79     }
  80   }
  81
  82   if (n_aa <= 0) {
  83     printf ("\n ERROR! There is no 'AA OSEC PSEC' data in submited prediction.");
  84     printf ("\n        The submission should contain an observed and predicted");
  85     printf ("\n        secondary structure in COLUMN format.\n");
  86     exit (0);
  87   }
  88
  89   printf ("\n\n SECONDARY STRUCTURE PREDICTION");
  90   printf ("\n NUMBER OF RESIDUES PREDICTED: LENGTH = %d", n_aa);
  91   printf ("\n AA  OSEC  PSEC     NUM");
  92   for (i = 0; i < n_aa; i++) {
  93     printf ("\n  %1c   %1c     %1c      %4d", aa[i], osec[i], psec[i], i + 1);
  94   }
  95   printf ("\n -----------------------\n");
  96   printf ("\n SECONDARY STRUCTURE PREDICTION ACCURACY EVALUATION.  N_AA = %4d\n", n_aa);
  97   if (pdata.sov_out >= 1) {
  98     printf ("\n SOV parameters:   DELTA = %5.2f   DELTA-S = %5.2f\n", pdata.sov_delta, pdata.sov_delta_s);
  99   }
 100
 101   printf ("\n                                   ALL    HELIX   STRAND     COIL\n");
 102
 103   pdata.q3_what = 0;
 104   out0 = q3 (n_aa, osec, psec, &pdata);
 105   pdata.q3_what = 1;
 106   out1 = q3 (n_aa, osec, psec, &pdata);
 107   pdata.q3_what = 2;
 108   out2 = q3 (n_aa, osec, psec, &pdata);
 109   pdata.q3_what = 3;
 110   out3 = q3 (n_aa, osec, psec, &pdata);
 111   printf ("\n Q3                         :   %6.1f   %6.1f   %6.1f   %6.1f", out0 * 100.0, out1 * 100.0, out2 * 100.0, out3 * 100.0);
 112   printf ("\n");
 113
 114   sov_method = pdata.sov_method;
 115
 116   if (sov_method != 0)
 117     pdata.sov_method = 1;
 118
 119   if (pdata.sov_method == 1) {
 120     pdata.sov_what = 0;
 121     out0 = sov (n_aa, osec, psec, &pdata);
 122     pdata.sov_what = 1;
 123     out1 = sov (n_aa, osec, psec, &pdata);
 124     pdata.sov_what = 2;
 125     out2 = sov (n_aa, osec, psec, &pdata);
 126     pdata.sov_what = 3;
 127     out3 = sov (n_aa, osec, psec, &pdata);
 128     printf ("\n SOV                        :   %6.1f   %6.1f   %6.1f   %6.1f", out0 * 100.0, out1 * 100.0, out2 * 100.0, out3 * 100.0);
 129     printf ("\n");
 130   }
 131
 132   if (sov_method != 1)
 133     pdata.sov_method = 0;
 134
 135   if (pdata.sov_method == 0) {
 136     pdata.sov_delta = 1.0;
 137
 138     pdata.sov_what = 0;
 139     out0 = sov (n_aa, osec, psec, &pdata);
 140     pdata.sov_what = 1;
 141     out1 = sov (n_aa, osec, psec, &pdata);
 142     pdata.sov_what = 2;
 143     out2 = sov (n_aa, osec, psec, &pdata);
 144     pdata.sov_what = 3;
 145     out3 = sov (n_aa, osec, psec, &pdata);
 146     printf ("\n SOV (1994 JMB. [delta=50%%]):   %6.1f   %6.1f   %6.1f   %6.1f", out0 * 100.0, out1 * 100.0, out2 * 100.0, out3 * 100.0);
 147
 148     pdata.sov_delta = 0.0;
 149
 150     pdata.sov_what = 0;
 151     out0 = sov (n_aa, osec, psec, &pdata);
 152     pdata.sov_what = 1;
 153     out1 = sov (n_aa, osec, psec, &pdata);
 154     pdata.sov_what = 2;
 155     out2 = sov (n_aa, osec, psec, &pdata);
 156     pdata.sov_what = 3;
 157     out3 = sov (n_aa, osec, psec, &pdata);
 158     printf ("\n SOV (1994 JMB. [delta=0])  :   %6.1f   %6.1f   %6.1f   %6.1f", out0 * 100.0, out1 * 100.0, out2 * 100.0, out3 * 100.0);
 159
 160     printf ("\n");
 161   }
 162
 163   printf ("\n -----------------------\n");
 164
 165   exit (0);
 166 }
 167
 168 /*-----------------------------------------------------------
 169 /
 170 /   check_aa - checks an amino acid
 171 /
 172 /------------------------------------------------------------*/
 173 int
 174 check_aa (char token, char *letter, int n)
 175 {
 176   int i;
 177
 178   for (i = 0; i < n; i++) {
 179     if (letter[i] == token)
 180       return i;
 181   }
 182   return n;
 183 }
 184
 185 /*-----------------------------------------------------------
 186 /
 187 /   read_aa_osec_psec - read secondary structure segments file
 188 /
 189 /------------------------------------------------------------*/
 190 int
 191 read_aa_osec_psec (char aa[MAXRES], char sss1[MAXRES], char sss2[MAXRES], parameters * pdata, char *letter)
 192 {
 193   int i, j, coil, n_aa, n_aa_1, n_aa_2, n_aa_3, f_seq;
 194   float x;
 195   char line[1000], keyword[250], first[250], second[250], third[250], junk[250];
 196   FILE *fp;
 197
 198   if ((fp = fopen (pdata->fname, "r")) == NULL) {
 199     printf ("\n# error opening file %s for read\n\n", pdata->fname);
 200     exit (0);
 201   }
 202
 203   f_seq = 0;
 204   pdata->input = 0;
 205   n_aa = 0;
 206   n_aa_1 = 0;
 207   n_aa_2 = 0;
 208   n_aa_3 = 0;
 209   coil = 0;
 210
 211   while (fgets (line, 1000, fp) != NULL) {
 212     strcpy (keyword, "   ");
 213     strcpy (first, "   ");
 214     strcpy (second, "   ");
 215     strcpy (third, "   ");
 216     strcpy (junk, "   ");
 217     i = 0;
 218     j = 0;
 219     while (line[i] == ' ' && line[i] != '\n' && line[i] != '\0' && i < 250)
 220       i++;
 221     if (i < 250) {
 222       j = i;
 223       while (line[i] != ' ' && line[i] != '\n' && line[i] != '\0' && i < 250)
 224         i++;
 225     }
 226     j = i - j;
 227     if (j < 250 && j > 0) {
 228       sscanf (line, "%s", keyword);
 229     }
 230     if (!strncmp (keyword, "#", 1)) {
 231     } else if (!strncmp (keyword, "-----", 5)) {
 232     } else if (!strncmp (keyword, "NUMBER\0", 7)) {
 233     } else if (!strncmp (keyword, "SECONDARY\0", 10)) {
 234     } else if (!strncmp (keyword, "END\0", 4) && f_seq == 0) {
 235       fclose (fp);
 236       return n_aa;
 237     } else if (!strncmp (keyword, "AA-OSEC-PSEC\0", 13)) {
 238       printf ("%s", line);
 239       sscanf (line, "%s %s", keyword, first);
 240       strcpy (pdata->fname, first);
 241       pdata->input = 1;
 242     } else if (line[0] == '\n' || !strncmp (keyword, "   \0", 4)) {
 243     } else if (!strncmp (keyword, "AA\0", 3) && f_seq == 0) {
 244       sscanf (line, "%s %s %s", keyword, first, second);
 245       if (!strncmp (keyword, "AA\0", 3) && !strncmp (first, "PSEC\0", 5) && !strncmp (second, "OSEC\0", 5)) {
 246         pdata->order = 1;
 247       }
 248     } else if (!strncmp (keyword, "SOV-DELTA\0", 10)) {
 249       printf ("%s", line);
 250       sscanf (line, "%s %f", keyword, &x);
 251       pdata->sov_delta = x;
 252     } else if (!strncmp (keyword, "SOV-DELTA-S\0", 12)) {
 253       printf ("%s", line);
 254       sscanf (line, "%s %f", keyword, &x);
 255       pdata->sov_delta_s = x;
 256     } else if (!strncmp (keyword, "SOV-METHOD\0", 9)) {
 257       printf ("%s", line);
 258       sscanf (line, "%s %d", keyword, &i);
 259       pdata->sov_method = i;
 260     } else if (!strncmp (keyword, "SOV-OUTPUT\0", 9)) {
 261       printf ("%s", line);
 262       sscanf (line, "%s %d", keyword, &i);
 263       pdata->sov_out = i;
 264     } else if (line[0] == '>') {
 265       printf ("%s", line);
 266       if (f_seq < 2)
 267         n_aa = 0;
 268       f_seq++;
 269     } else if (f_seq == 0) {
 270       if (j > 1) {
 271         if (!strncmp (keyword, "SSP\0", 4)) {
 272           sscanf (line, "%s %s %s %s %s", keyword, junk, first, second, third);
 273         } else {
 274           printf ("\n ERROR! (line: %d) Check COLUMN format of your prediction!\n", n_aa + 1);
 275           fclose (fp);
 276           exit (0);
 277         }
 278       } else {
 279         sscanf (line, "%s %s %s", first, second, third);
 280       }
 281       aa[n_aa] = first[0];
 282       sss1[n_aa] = second[0];
 283       sss2[n_aa] = third[0];
 284       if (check_aa (aa[n_aa], letter, 23) == 23) {
 285         printf ("\n# ERROR!\n%s", line);
 286         printf ("\n# ERROR! (line: %d) Check amino acid code  %c\n", n_aa + 1, aa[n_aa]);
 287         fclose (fp);
 288         exit (0);
 289       }
 290       if (sss1[n_aa] == ' ' || sss2[n_aa] == ' ') {
 291         printf ("\n# ERROR!\n%s", line);
 292         printf ("\n# ERROR! (line: %d) Check secondary structure code\n", n_aa + 1);
 293         fclose (fp);
 294         exit (0);
 295       }
 296       if (sss1[n_aa] == 'L') {
 297         sss1[n_aa] = 'C';
 298         if (coil == 0) {
 299           printf ("\n# WARNING! (line: %d) The 'L' characters are interpreted as 'C' (coil)", n_aa + 1);
 300           coil = 1;
 301         }
 302       }
 303       if (sss2[n_aa] == 'L') {
 304         sss2[n_aa] = 'C';
 305         if (coil == 0) {
 306           printf ("\n# WARNING! (line: %d) The 'L' characters are interpreted as 'C' (coil)", n_aa + 1);
 307           coil = 1;
 308         }
 309       }
 310       if (sss1[n_aa] != 'C' && sss1[n_aa] != 'E' && sss1[n_aa] != 'H') {
 311         printf ("\n# ERROR!\n%s", line);
 312         printf ("\n# ERROR! (line: %d) Check secondary structure code  %c\n", n_aa + 1, sss1[n_aa]);
 313         fclose (fp);
 314         exit (0);
 315       }
 316       if (sss2[n_aa] != 'C' && sss2[n_aa] != 'E' && sss2[n_aa] != 'H') {
 317         printf ("\n# ERROR!\n%s", line);
 318         printf ("\n# ERROR! (line: %d) Check secondary structure code  %c\n", n_aa + 1, sss2[n_aa]);
 319         fclose (fp);
 320         exit (0);
 321       }
 322       n_aa++;
 323       if (n_aa >= MAXRES) {
 324         printf ("\n# ERROR! Check number of amino acid lines. (MAX = %d)\n\n", MAXRES);
 325         fclose (fp);
 326         exit (0);
 327       }
 328     } else if (f_seq == 1) {
 329       i = 0;
 330       while (line[i] != '\n') {
 331         if (line[i] != ' ' && line[i] != '\t' && line[i] != '\0' && line[i] != '\a' && line[i] != '\b' && line[i] != '\f' && line[i] != '\r' && line[i] != '\v' && i < 1000) {
 332           aa[n_aa] = 'X';
 333           sss1[n_aa] = line[i];
 334           if (sss1[n_aa] == 'L') {
 335             sss1[n_aa] = 'C';
 336             if (coil == 0) {
 337               printf ("\n# WARNING! The 'L' characters are interpreted as 'C' (coil)");
 338               coil = 1;
 339             }
 340           }
 341           if (sss1[n_aa] != 'C' && sss1[n_aa] != 'E' && sss1[n_aa] != 'H') {
 342             printf ("\n# ERROR!\n%s", line);
 343             printf ("\n# ERROR! Check secondary structure code: %c\n", sss1[n_aa]);
 344             fclose (fp);
 345             exit (0);
 346           }
 347           n_aa++;
 348           if (n_aa >= MAXRES) {
 349             printf ("\n# ERROR! Check number of residues. (MAX = %d)\n\n", MAXRES);
 350             fclose (fp);
 351             exit (0);
 352           }
 353         }
 354         i++;
 355       }
 356       n_aa_1 = n_aa;
 357     } else if (f_seq == 2) {
 358       i = 0;
 359       while (line[i] != '\n') {
 360         if (line[i] != ' ' && line[i] != '\t' && line[i] != '\0' && line[i] != '\a' && line[i] != '\b' && line[i] != '\f' && line[i] != '\r' && line[i] != '\v' && i < 1000) {
 361           aa[n_aa] = 'X';
 362           sss2[n_aa] = line[i];
 363           if (sss2[n_aa] == 'L') {
 364             sss2[n_aa] = 'C';
 365             if (coil == 0) {
 366               printf ("\n# WARNING! The 'L' characters are interpreted as 'C' (coil)");
 367               coil = 1;
 368             }
 369           }
 370           if (sss2[n_aa] != 'C' && sss2[n_aa] != 'E' && sss2[n_aa] != 'H') {
 371             printf ("\n# ERROR!\n%s", line);
 372             printf ("\n# ERROR! Check secondary structure code: %c\n", sss2[n_aa]);
 373             fclose (fp);
 374             exit (0);
 375           }
 376           n_aa++;
 377           if (n_aa >= MAXRES) {
 378             printf ("\n# ERROR! Check number of residues. (MAX = %d)\n\n", MAXRES);
 379             fclose (fp);
 380             exit (0);
 381           }
 382         }
 383         i++;
 384       }
 385       n_aa_2 = n_aa;
 386     } else if (f_seq == 3) {
 387       i = 0;
 388       while (line[i] != '\n') {
 389         if (line[i] != ' ' && line[i] != '\t' && line[i] != '\0' && line[i] != '\a' && line[i] != '\b' && line[i] != '\f' && line[i] != '\r' && line[i] != '\v' && i < 1000) {
 390           aa[n_aa_3] = line[i];
 391           if (check_aa (aa[n_aa_3], letter, 23) == 23) {
 392             printf ("\n# ERROR!\n%s", line);
 393             printf ("\n# ERROR! (N_res: %d) Check amino acid code  %c\n", n_aa_3 + 1, aa[n_aa_3]);
 394             fclose (fp);
 395             exit (0);
 396           }
 397           n_aa_3++;
 398           if (n_aa_3 >= MAXRES) {
 399             printf ("\n# ERROR! Check number of residues. (MAX = %d)\n\n", MAXRES);
 400             fclose (fp);
 401             exit (0);
 402           }
 403         }
 404         i++;
 405       }
 406     }
 407   }
 408   if (n_aa_1 != n_aa_2) {
 409     printf ("\n# ERROR! Check format of your submission.\n");
 410     fclose (fp);
 411     exit (0);
 412   }
 413   return n_aa;
 414 }
 415
 416 /*-------------------------------------------------------------
 417 /
 418 /  default_parameters - default parameters for SOV program
 419 /
 420 /------------------------------------------------------------*/
 421 void
 422 default_parameters (parameters * pdata)
 423 {
 424   pdata->input = 0;
 425   pdata->order = 0;
 426   pdata->sov_method = 1;
 427   pdata->sov_delta = 1.0;
 428   pdata->sov_delta_s = 0.5;
 429   pdata->sov_out = 0;
 430
 431   return;
 432 }
 433
 434 /*-----------------------------------------------------------
 435 /
 436 /    sov - evaluate SSp by the Segment OVerlap quantity (SOV)
 437 /                Input: secondary structure segments
 438 /
 439 /------------------------------------------------------------*/
 440 float
 441 sov (int n_aa, char sss1[MAXRES], char sss2[MAXRES], parameters * pdata)
 442 {
 443   int i, k, length1, length2, beg_s1, end_s1, beg_s2, end_s2;
 444   int j1, j2, k1, k2, minov, maxov, d, d1, d2, n, multiple;
 445   char s1, s2, sse[3];
 446   float out;
 447   double s, x;
 448
 449   sse[0] = '#';
 450   sse[1] = '#';
 451   sse[2] = '#';
 452
 453   if (pdata->sov_what == 0) {
 454     sse[0] = 'H';
 455     sse[1] = 'E';
 456     sse[2] = 'C';
 457   }
 458   if (pdata->sov_what == 1) {
 459     sse[0] = 'H';
 460     sse[1] = 'H';
 461     sse[2] = 'H';
 462   }
 463   if (pdata->sov_what == 2) {
 464     sse[0] = 'E';
 465     sse[1] = 'E';
 466     sse[2] = 'E';
 467   }
 468   if (pdata->sov_what == 3) {
 469     sse[0] = 'C';
 470     sse[1] = 'C';
 471     sse[2] = 'C';
 472   }
 473   n = 0;
 474   for (i = 0; i < n_aa; i++) {
 475     s1 = sss1[i];
 476     if (s1 == sse[0] || s1 == sse[1] || s1 == sse[2]) {
 477       n++;
 478     }
 479   }
 480   out = 0.0;
 481   s = 0.0;
 482   length1 = 0;
 483   length2 = 0;
 484   i = 0;
 485   while (i < n_aa) {
 486     beg_s1 = i;
 487     s1 = sss1[i];
 488     while (sss1[i] == s1 && i < n_aa) {
 489       i++;
 490     }
 491     end_s1 = i - 1;
 492     length1 = end_s1 - beg_s1 + 1;
 493     multiple = 0;
 494     k = 0;
 495     while (k < n_aa) {
 496       beg_s2 = k;
 497       s2 = sss2[k];
 498       while (sss2[k] == s2 && k < n_aa) {
 499         k++;
 500       }
 501       end_s2 = k - 1;
 502       length2 = end_s2 - beg_s2 + 1;
 503       if (s1 == sse[0] || s1 == sse[1] || s1 == sse[2]) {
 504         if (s1 == s2 && end_s2 >= beg_s1 && beg_s2 <= end_s1) {
 505           if (multiple > 0 && pdata->sov_method == 1) {
 506             n = n + length1;
 507           }
 508           multiple++;
 509           if (beg_s1 > beg_s2) {
 510             j1 = beg_s1;
 511             j2 = beg_s2;
 512           } else {
 513             j1 = beg_s2;
 514             j2 = beg_s1;
 515           }
 516           if (end_s1 < end_s2) {
 517             k1 = end_s1;
 518             k2 = end_s2;
 519           } else {
 520             k1 = end_s2;
 521             k2 = end_s1;
 522           }
 523           minov = k1 - j1 + 1;
 524           maxov = k2 - j2 + 1;
 525           d1 = floor (length1 * pdata->sov_delta_s);
 526           d2 = floor (length2 * pdata->sov_delta_s);
 527           if (d1 > d2)
 528             d = d2;
 529           if (d1 <= d2 || pdata->sov_method == 0)
 530             d = d1;
 531           if (d > minov) {
 532             d = minov;
 533           }
 534           if (d > maxov - minov) {
 535             d = maxov - minov;
 536           }
 537           x = pdata->sov_delta * d;
 538           x = (minov + x) * length1;
 539           if (maxov > 0) {
 540             s = s + x / maxov;
 541           } else {
 542             printf ("\n ERROR! minov = %-4d maxov = %-4d length = %-4d d = %-4d   %4d %4d  %4d %4d", minov, maxov, length1, d, beg_s1 + 1, end_s1 + 1, beg_s2 + 1, end_s2 + 1);
 543           }
 544           if (pdata->sov_out == 2) {
 545             printf ("\n TEST: minov = %-4d maxov = %-4d length = %-4d d = %-4d   %4d %4d  %4d %4d", minov, maxov, length1, d, beg_s1 + 1, end_s1 + 1, beg_s2 + 1, end_s2 + 1);
 546           }
 547         }
 548       }
 549     }
 550   }
 551   if (pdata->sov_out == 2) {
 552     printf ("\n TEST: Number of considered residues = %d", n);
 553   }
 554   if (n > 0) {
 555     out = s / n;
 556   } else {
 557     out = 1.0;
 558   }
 559   return out;
 560 }
 561
 562 /*-----------------------------------------------------------
 563 /
 564 /    Q3 - evaluate SSp by the residues predicted correctly (Q3)
 565 /                Input: secondary structure segments
 566 /
 567 /------------------------------------------------------------*/
 568 float
 569 q3 (int n_aa, char sss1[MAXRES], char sss2[MAXRES], parameters * pdata)
 570 {
 571   int i, n;
 572   float out;
 573   char s, sse[3];
 574
 575   sse[0] = '#';
 576   sse[1] = '#';
 577   sse[2] = '#';
 578
 579   if (pdata->q3_what == 0) {
 580     sse[0] = 'H';
 581     sse[1] = 'E';
 582     sse[2] = 'C';
 583   }
 584   if (pdata->q3_what == 1) {
 585     sse[0] = 'H';
 586     sse[1] = 'H';
 587     sse[2] = 'H';
 588   }
 589   if (pdata->q3_what == 2) {
 590     sse[0] = 'E';
 591     sse[1] = 'E';
 592     sse[2] = 'E';
 593   }
 594   if (pdata->q3_what == 3) {
 595     sse[0] = 'C';
 596     sse[1] = 'C';
 597     sse[2] = 'C';
 598   }
 599
 600   n = 0;
 601   out = 0.0;
 602   for (i = 0; i < n_aa; i++) {
 603     s = sss1[i];
 604     if (s == sse[0] || s == sse[1] || s == sse[2]) {
 605       n++;
 606       if (sss1[i] == sss2[i]) {
 607         out = out + 1.0;
 608       }
 609     }
 610   }
 611   if (n > 0) {
 612     out = out / n;
 613   } else {
 614     out = 1.0;
 615   }
 616
 617   return out;
 618 }