source: tags/initial/TREEGEN/rns.c

Last change on this file was 2, checked in by oldcode, 23 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.5 KB
Line 
1#include "rns.h"
2
3#ifndef __SPREADIN_H
4    #include "spreadin.h"
5#endif
6#ifndef __LIMITS_H
7    #include <limits.h>
8#endif
9
10/* /------------------------\ */
11/* |  Erzeugung der Ur-RNS  | */
12/* \------------------------/ */
13
14int        orgLen;         /* L„nge der Ur-RNS */
15double     orgHelixPart;   /* Anteil Helix-Bereich */
16static int rnsCreated; /* Anzahl bisher erzeugter RNS-Sequenzen */
17
18/* /------------\ */
19/* |  Mutation  | */
20/* \------------/ */
21
22int    timeSteps;        /* Anzahl Zeitschritte */
23Frand  mrpb_Init,        /* Initialisierungsfunktion fr 'mutationRatePerBase' */
24       l2hrpb_Init,      /* Initialisierungsfunktion fr 'loop2helixRatePerBase' */
25       pairPart,         /* Anteil paarender Helix-Bindungen */
26       mutationRate,     /* Mutationsrate */
27       splitRate,        /* Spaltungsrate */
28       helixGcDruck,     /* G-C-Druck       im Helix-Bereich */
29       helixGcRate,      /* Verh„ltnis G:C  im Helix-Bereich */
30       helixAtRate,      /* Verh„ltnis A:T  im Helix-Bereich */
31       loopGcDruck,      /* G-C-Druck       im Loop-Bereich */
32       loopGcRate,       /* Verh„ltnis G:C  im Loop-Bereich */
33       loopAtRate;       /* Verh„ltnis A:T  im Loop-Bereich */
34double transitionRate,   /* Transition-Rate */
35       transversionRate; /* Transversion-Rate */
36
37static double     *mutationRatePerBase,   /* positionsspez. Mutationsrate (wird nur einmal bestimmt und bleibt dann konstant) */
38                  *loop2helixRatePerBase; /* positionsspez. Rate fr Wechsel Loop-Base in Helix-Base und vv. (wird nur einmal bestimmt und bleibt dann konstant) */
39static int         mrpb_anz,              /* Anzahl Positionen */
40                   mrpb_allocated,        /* wirklich Gr”áe des Arrays */
41                   l2hrpb_anz,            /* Anzahl Positionen */
42                   l2hrpb_allocated;      /* wirklich Gr”áe des Arrays */
43static DoubleProb  helixMutationMatrix,   /* Mutationsmatrix fr Helix-Bereiche */
44                   loopMutationMatrix;    /* Mutationsmatrix fr Loop-Bereiche */
45
46/* /----------------------\ */
47/* |  Ausgabefilepointer  | */
48/* \----------------------/ */
49
50FILE *topo, /* Topologie */
51     *seq;  /* Sequenzen */
52
53/* /-------------\ */
54/* |  Sonstiges  | */
55/* \-------------/ */
56
57static int minDepth = INT_MAX, /* minimale Tiefe (Astanzahl) der Blattspitzen */
58           maxDepth = INT_MIN; /* maximale Tiefe der Blattspitzen */
59
60/* -------------------------------------------------------------------------- */
61/*      void dumpDepths(void) */
62/* ------------------------------------------------------ 24.05.95 22.27 ---- */
63void dumpDepths(void)
64{
65    printf("Minimale Baumtiefe = %i\n", minDepth);
66    printf("Maximale Baumtiefe = %i\n", maxDepth);
67}
68/* -------------------------------------------------------------------------- */
69/*      static void dumpRNS(RNS rns) */
70/* ------------------------------------------------------ 26.05.95 11.29 ---- */
71static void dumpRNS(RNS rns)
72{
73    int        b,
74               b1,
75               b2;
76    static int cleared,
77               h_cnt[BASETYPES+1][BASETYPES+1],
78               l_cnt[BASETYPES+1],
79               loop,
80               helix;
81
82    if (!cleared)
83    {
84        for (b1 = 0; b1<(BASETYPES+1); b1++)
85        {
86            for (b2 = 0; b2<(BASETYPES+1); b2++) h_cnt[b1][b2] = 0;
87            l_cnt[b1] = 0;
88        }
89
90        loop  = 0;
91        helix = 0;
92
93        cleared = 1;
94    }
95
96    if (rns)
97    {
98        for (b = 0; b<(rns->bases); b++)
99        {
100            char base = rns->base[b];
101
102            if (isHelical(base))
103            {
104                int bt1 = char2BaseType(base),
105                    bt2 = char2BaseType(rns->base[b+1]);
106
107                h_cnt[bt1][bt2]++;
108                helix++;
109                b++;
110            }
111            else
112            {
113                int bt = char2BaseType(base);
114
115                l_cnt[bt]++;
116                loop++;
117            }
118        }
119    }
120    else
121    {
122        printf("Helix-Basenpaare = %i\n"
123               "Loop-Basen       = %i\n"
124               "Helix:Loop       = %f\n",
125               helix,
126               loop,
127               (double)helix/(double)loop);
128
129        {
130            int gc      = h_cnt[BASE_C][BASE_G]+h_cnt[BASE_G][BASE_C],
131                at      = h_cnt[BASE_A][BASE_T]+h_cnt[BASE_T][BASE_A],
132                paarend = gc+at;
133
134            printf("GC-Paare              = %i\n"
135                   "AT-Paare              = %i\n"
136                   "Paare:Helix-Bindungen = %f\n"
137                   "GC-Paare:Paare        = %f\n",
138                   gc,
139                   at,
140                   (double)paarend/(double)helix,
141                   (double)gc/(double)paarend);
142        }
143
144        printf("\n");
145    }
146}
147/* -------------------------------------------------------------------------- */
148/*      static void initBaseSpecificProbs(int bases) */
149/* ------------------------------------------------------ 24.05.95 12.51 ---- */
150static void initBaseSpecificProbs(int bases)
151{
152    int b;
153
154    mrpb_anz            = bases;
155    mrpb_allocated      = bases;
156    mutationRatePerBase = malloc(bases*sizeof(double));
157
158    l2hrpb_anz            = bases;
159    l2hrpb_allocated      = bases;
160    loop2helixRatePerBase = malloc(bases*sizeof(double));
161
162    if (!mutationRatePerBase || !loop2helixRatePerBase) outOfMemory();
163
164    for (b = 0; b<bases; b++)
165    {
166        mutationRatePerBase[b]   = getFrand(mrpb_Init);
167        loop2helixRatePerBase[b] = getFrand(l2hrpb_Init);
168    }
169}
170/* -------------------------------------------------------------------------- */
171/*      static RNS allocRNS(int len) */
172/* ------------------------------------------------------ 20.05.95 16.04 ---- */
173static RNS allocRNS(int len)
174{
175    RNS rns = malloc(sizeof(*rns));
176
177    if (!rns) outOfMemory();
178
179    rns->bases = orgLen;
180    rns->base  = malloc(sizeof(*(rns->base))*orgLen);
181
182    if (!rns->base) outOfMemory();
183
184    return rns;
185}
186/* -------------------------------------------------------------------------- */
187/*      RNS createOriginRNS(void) */
188/* ------------------------------------------------------ 14.05.95 14:54 ---- */
189/* */
190/*  Erzeugt eine Ur-RNS */
191/* */
192RNS createOriginRNS(void)
193{
194    RNS rns      = allocRNS(orgLen);
195    int helixLen = orgLen*orgHelixPart,
196        l;
197    str base     = rns->base;
198
199    initBaseSpecificProbs(orgLen);
200
201    rns->laufNr = rnsCreated++;
202
203    /* /------------------\ */
204    /* |  Helix erzeugen  |                                                             im Loop-Bereich */
205    /* \------------------/ */
206
207    if (helixLen%1) helixLen--;                            /* muá gerade L„nge haben, da nur Paare! */
208
209    assert(helixLen<=orgLen);
210
211    rns->helix   = helixLen/2;
212    rns->pairing = 0;
213
214    {
215        DoubleProb orgHelixProb;
216        Spreading  s;
217        int        b1,
218                   b2;
219        double     actPairPart     = getFrand(pairPart),
220                   actHelixGcDruck = getFrand(helixGcDruck),
221                   actHelixGcRate  = getFrand(helixGcRate),
222                   actHelixAtRate  = getFrand(helixAtRate),
223                   nonPairProb     = (1.0-actPairPart)/2.0;
224
225        for (b1 = 0; b1<BASETYPES; b1++)
226        {
227            for (b2 = 0; b2<BASETYPES; b2++)
228            {
229                if (isPairing(b1, b2))
230                {
231                    switch (b1)
232                    {
233                        case BASE_A:
234                        case BASE_T:
235                        {
236                            orgHelixProb[b1][b2] = (actPairPart*(1.0-actHelixGcDruck))/2.0;
237                            break;
238                        }
239                        case BASE_C:
240                        case BASE_G:
241                        {
242                            orgHelixProb[b1][b2] = (actPairPart*actHelixGcDruck)/2.0;
243                            break;
244                        }
245                    }
246                }
247                else
248                {
249                    double prob = nonPairProb;
250                    int    b    = b1;
251
252                    while (1)                              /* wird je einmal mit b1 und b2 ausgefhrt */
253                    {
254                        switch (b)
255                        {
256                            case BASE_A:
257                            {
258                                prob = prob*(1.0-actHelixGcDruck)*actHelixAtRate;
259                                break;
260                            }
261                            case BASE_C:
262                            {
263                                prob = prob*actHelixGcDruck*(1.0-actHelixGcRate);
264                                break;
265                            }
266                            case BASE_G:
267                            {
268                                prob = prob*actHelixGcDruck*actHelixGcRate;
269                                break;
270                            }
271                            case BASE_T:
272                            {
273                                prob = prob*(1.0-actHelixGcDruck)*(1.0-actHelixAtRate);
274                                break;
275                            }
276                        }
277
278                        if (b==b2) break;
279                        b = b2;
280                    }
281
282                    orgHelixProb[b1][b2] = prob;
283                }
284            }
285        }
286
287        s = newSpreading((double*)orgHelixProb, BASEQUAD);
288
289        for (l = 0; l<helixLen; l+=2)
290        {
291            int val = spreadRand(s),
292                B1  = val%BASETYPES,
293                B2  = val/BASETYPES;
294
295            base[l]   = helixBaseChar[B1];
296            base[l+1] = helixBaseChar[B2];
297
298            rns->pairing += isPairing(B1, B2);
299        }
300
301        freeSpreading(s);
302    }
303
304    /* /-----------------\ */
305    /* |  Loop erzeugen  | */
306    /* \-----------------/ */
307
308    {
309        SingleProb orgLoopProb;
310        Spreading  s;
311        double     actLoopGcDruck = getFrand(loopGcDruck),
312                   actLoopGcRate  = getFrand(loopGcRate),
313                   actLoopAtRate  = getFrand(loopAtRate);
314
315        orgLoopProb[BASE_A] = (1.0-actLoopGcDruck)*actLoopAtRate;
316        orgLoopProb[BASE_C] = actLoopGcDruck*(1.0-actLoopGcRate);
317        orgLoopProb[BASE_G] = actLoopGcDruck*actLoopGcRate;
318        orgLoopProb[BASE_T] = (1.0-actLoopGcDruck)*(1.0-actLoopAtRate);
319
320        s = newSpreading((double*)orgLoopProb, BASETYPES);
321        for (; l<orgLen; l++) base[l] = loopBaseChar[spreadRand(s)];
322        freeSpreading(s);
323    }
324
325    return rns;
326}
327/* -------------------------------------------------------------------------- */
328/*      void freeRNS(RNS rns) */
329/* ------------------------------------------------------ 20.05.95 19.45 ---- */
330void freeRNS(RNS rns)
331{
332    free(rns->base);
333    free(rns);
334}
335/* -------------------------------------------------------------------------- */
336/*      static RNS dupRNS(RNS rns) */
337/* ------------------------------------------------------ 20.05.95 20.32 ---- */
338static RNS dupRNS(RNS rns)
339{
340    RNS neu = malloc(sizeof(*rns));
341
342    if (!neu) outOfMemory();
343
344    memcpy(neu, rns, sizeof(*rns));
345
346    neu->base = malloc(rns->bases*sizeof(*(neu->base)));
347    memcpy(neu->base, rns->base, rns->bases);
348
349    neu->laufNr = rnsCreated++;
350
351    return neu;
352}
353/* -------------------------------------------------------------------------- */
354/*      static void dumpDoubleProb(double *d, int anz) */
355/* ------------------------------------------------------ 25.05.95 01.31 ---- */
356/*static void dumpDoubleProb(double *d, int anz) */
357/*{ */
358/*    while (anz--) printf("%-10f", *d++); */
359/*    printf("\n\n"); */
360/*} */
361/* -------------------------------------------------------------------------- */
362/*      static void calcMutationMatrix(DoubleProb mutationMatrix, double mu... */
363/* ------------------------------------------------------ 24.05.95 13.58 ---- */
364static void calcMutationMatrix(DoubleProb mutationMatrix, double muteRate, double gcDruck, double gcRate, double atRate, double *pairProb)
365{
366    double k   = transitionRate/transversionRate,
367           fa  = (1.0-gcDruck)*atRate,
368           fc  = gcDruck*(1.0-gcRate),
369           fg  = gcDruck*gcRate,
370           ft  = (1.0-gcDruck)*(1.0-atRate),
371           bfa = transversionRate*fa,
372           bfc = transversionRate*fc,
373           bfg = transversionRate*fg,
374           bft = transversionRate*ft,
375           kag = k/(fa+fg),
376           kct = k/(fc+ft);
377/*           sa  = (kag+3.0)*bfa,                            // Summe der "mutierenden" Positionen jeder Zeile */
378/*           sc  = (kct+3.0)*bfc, */
379/*           sg  = (kag+3.0)*bfg, */
380/*           st  = (kct+3.0)*bft; */
381
382    /* Auf aktuelle Mutationsrate normieren */
383
384/*    bfa = bfa*muteRate/sa; */
385/*    bfc = bfc*muteRate/sc; */
386/*    bfg = bfg*muteRate/sg; */
387/*    bft = bft*muteRate/st; */
388
389/*    printf("bfa=%f\n", bfa); */
390/*    printf("bfc=%f\n", bfc); */
391/*    printf("bfg=%f\n", bfg); */
392/*    printf("bft=%f\n", bft); */
393
394    /* Matrix besetzen */
395
396    mutationMatrix[BASE_A][BASE_A] = 1.0-(kag+3.0)*bfa;
397    mutationMatrix[BASE_C][BASE_A] = bfa;
398    mutationMatrix[BASE_G][BASE_A] = (kag+1.0)*bfa;
399    mutationMatrix[BASE_T][BASE_A] = bfa;
400
401    mutationMatrix[BASE_A][BASE_C] = bfc;
402    mutationMatrix[BASE_C][BASE_C] = 1.0-(kct+3.0)*bfc;
403    mutationMatrix[BASE_G][BASE_C] = bfc;
404    mutationMatrix[BASE_T][BASE_C] = (kct+1.0)*bfc;
405
406    mutationMatrix[BASE_A][BASE_G] = (kag+1.0)*bfg;
407    mutationMatrix[BASE_C][BASE_G] = bfg;
408    mutationMatrix[BASE_G][BASE_G] = 1.0-(kag+3.0)*bfg;
409    mutationMatrix[BASE_T][BASE_G] = bfg;
410
411    mutationMatrix[BASE_A][BASE_T] = bft;
412    mutationMatrix[BASE_C][BASE_T] = (kct+1.0)*bft;
413    mutationMatrix[BASE_G][BASE_T] = bft;
414    mutationMatrix[BASE_T][BASE_T] = 1.0-(kct+3.0)*bft;
415
416/*    { */
417/*        int von,                                         // Matrix ausgeben */
418/*            nach; */
419/* */
420/*        printf("       von %c     von %c     von %c     von %c     \n", */
421/*                helixBaseChar[0], */
422/*                helixBaseChar[1], */
423/*                helixBaseChar[2], */
424/*                helixBaseChar[3] ); */
425/* */
426/*        for (nach = BASE_A; nach<=BASE_T; nach++) */
427/*        { */
428/*            double sum = 0.0; */
429/* */
430/*            printf("nach %c ", helixBaseChar[nach]); */
431/* */
432/*            for (von = BASE_A; von<=BASE_T; von++) */
433/*            { */
434/*                printf("%-10f", mutationMatrix[von][nach]); */
435/*                sum += mutationMatrix[von][nach]; */
436/*            } */
437/* */
438/*            printf("    sum = %-10f\n", sum); */
439/*        } */
440/* */
441/*        printf("\n"); */
442/*    } */
443
444    if (pairProb)                                          /* soll pairProb berechnet werden? */
445    {
446        double mutatesTo[BASETYPES],
447               freq[BASETYPES];                            /* H„ufigkeit der einzelnen Basen */
448        int    von,
449               nach;
450
451        freq[BASE_A] = fa;
452        freq[BASE_C] = fc;
453        freq[BASE_G] = fg;
454        freq[BASE_T] = ft;
455
456        for (nach = 0; nach<BASETYPES; nach++)
457            mutatesTo[nach] = 0.0;
458
459        for (von = 0; von<BASETYPES; von++)
460            for (nach = 0; nach<BASETYPES; nach++)
461                mutatesTo[nach] += mutationMatrix[von][nach]*freq[von];
462
463        *pairProb = 2.0*mutatesTo[BASE_A]*mutatesTo[BASE_T] + 2.0*mutatesTo[BASE_C]*mutatesTo[BASE_G];
464    }
465}
466/* -------------------------------------------------------------------------- */
467/*      static int calcPairTrials(double pairProb, double actPairPart) */
468/* ------------------------------------------------------ 25.05.95 13.31 ---- */
469/* */
470/*  Berechnet die Anzahl Mutations-Wiederholungen, die notwendig sind, um */
471/*  mindestens 'actPairPart' Prozent paarende Bindungen zu erhalten, falls */
472/*  die Wahrscheinlichkeit eine paarende Bindung zu erzeugen gleich */
473/*  'pairProb' ist. */
474/* */
475static int calcPairTrials(double pairProb, double actPairPart)
476{
477    int    trials   = 1;
478    double failProb = 1.0-pairProb,
479           succProb = pairProb;
480
481    while (succProb<actPairPart)
482    {
483        pairProb *= failProb;
484        succProb += pairProb;
485        trials++;
486
487/*        printf("trials=%i succProb=%f actPairPart=%f\n", trials, succProb, actPairPart); */
488    }
489
490    return trials;
491}
492/* -------------------------------------------------------------------------- */
493/*      static void indent(int depth) */
494/* ------------------------------------------------------ 24.05.95 21.08 ---- */
495/*static void indent(int depth) */
496/*{ */
497/*    while (depth--) fputc(' ', topo); */
498/*} */
499/* -------------------------------------------------------------------------- */
500/*      static void mutateRNS(RNS rns, int steps, int depth) */
501/* ------------------------------------------------------ 20.05.95 19.50 ---- */
502/* */
503/*  Mutiert eine RNS bis zur n„chsten Spaltung */
504/* */
505/*  'steps'     Anzahl noch zu durchlaufender Zeitschritte */
506/* */
507static void mutateRNS(RNS rns, int steps, int depth)
508{
509    int    splitInSteps,
510           s;
511    double mutationTime = 0.0;
512
513    /* /---------------------------------------\ */
514    /* |  Schritte bis zur Spaltung berechnen  | */
515    /* \---------------------------------------/ */
516
517    {
518        double actualSplitRate = getFrand(splitRate);
519
520        assert(actualSplitRate!=0);
521
522        splitInSteps = (int)(1.0/actualSplitRate);
523        if (splitInSteps>steps) splitInSteps = steps;
524
525        assert(splitInSteps>=1);
526    }
527
528    /* /----------------------------\ */
529    /* |  Zeitschritte durchlaufen  | */
530    /* \----------------------------/ */
531
532    for (s = 0; s<splitInSteps; s++)
533    {
534        int       b,
535                  pairTrials;                              /* Anzahl Versuche eine paarende Helixbindung herzustellen */
536        double    actMutationRate    = getFrand(mutationRate),
537                  actPairPart        = getFrand(pairPart);
538        Spreading s_helix[BASETYPES],
539                  s_loop[BASETYPES];
540
541        {
542            double pairProb;                               /* Wahrscheinlichkeit, daá ein Paar im helikalen Bereich entsteht */
543
544            calcMutationMatrix(helixMutationMatrix, 1.0, getFrand(helixGcDruck), getFrand(helixGcRate), getFrand(helixAtRate), &pairProb);
545            calcMutationMatrix(loopMutationMatrix, actMutationRate, getFrand(loopGcDruck), getFrand(loopGcRate), getFrand(loopAtRate), NULL);
546
547            pairTrials = calcPairTrials(pairProb, actPairPart);
548/*            printf("pairProb=%f pairTrials=%i\n", pairProb, pairTrials); */
549        }
550
551        for (b = 0; b<BASETYPES; b++)
552        {
553            s_helix[b] = newSpreading(&(helixMutationMatrix[b][0]), BASETYPES);
554            s_loop[b]  = newSpreading(&(loopMutationMatrix[b][0]), BASETYPES);
555        }
556
557        mutationTime += actMutationRate;                   /* Mutationszeit aufaddieren (Einheit ist Mutationsrate*Zeitschritte) */
558
559        /* /----------------------------------\ */
560        /* |  Alle Basen(-paare) durchlaufen  | */
561        /* \----------------------------------/ */
562
563        for (b = 0; b<(rns->bases); )
564        {
565            char base = rns->base[b];
566
567            if (!isDeleted(base))                          /* Deletes ignorieren */
568            {
569                /* /---------------------\ */
570                /* |  Helicale Bereiche  | */
571                /* \---------------------/ */
572
573                if (isHelical(base))
574                {
575                    int  trials = pairTrials,
576                         mut1   = randProb()<mutationRatePerBase[b]*actMutationRate,
577                         mut2   = randProb()<mutationRatePerBase[b+1]*actMutationRate;
578                    char base2  = rns->base[b+1];
579
580                    assert(isHelical(base2));
581
582                    if (mut1 || mut2)
583                    {
584                        int bt1 = char2BaseType(base),
585                            bt2 = char2BaseType(base2);
586
587                        if (isPairing(bt1, bt2))
588                        {
589                            rns->pairing--;
590                        }
591
592                        while (trials--)
593                        {
594                            if (mut1)
595                            {
596                                if (mut2)                  /* beide Basen mutieren */
597                                {
598                                    bt1 = spreadRand(s_helix[bt1]);
599                                    bt2 = spreadRand(s_helix[bt2]);
600                                }
601                                else                       /* nur 1.Base mutieren */
602                                {
603                                    bt1 = spreadRand(s_helix[bt1]);
604                                }
605                            }
606                            else                           /* nur 2.Base mutieren */
607                            {
608                                bt2 = spreadRand(s_helix[bt2]);
609                            }
610
611                            if (isPairing(bt1, bt2))       /* paarend? ja->abbrechen */
612                            {
613                                rns->pairing++;
614                                break;
615                            }
616                        }
617
618                        rns->base[b]   = helixBaseChar[bt1];
619                        rns->base[b+1] = helixBaseChar[bt2];
620                    }
621
622                    b++;
623                }
624
625                /* /-----------------\ */
626                /* |  Loop-Bereiche  | */
627                /* \-----------------/ */
628
629                else
630                {
631                    double mutationProb = actMutationRate*mutationRatePerBase[b];
632
633                    if (randProb()<mutationProb)
634                    {
635                        rns->base[b] = loopBaseChar[spreadRand(s_loop[char2BaseType(base)])];
636                    }
637                }
638            }
639
640            b++;
641        }
642
643        for (b = 0; b<BASETYPES; b++)
644        {
645            freeSpreading(s_helix[b]);
646            freeSpreading(s_loop[b]);
647        }
648    }
649
650    splitRNS(rns, mutationTime, steps-splitInSteps, depth+1);
651}
652/* -------------------------------------------------------------------------- */
653/*      void splitRNS(RNS origin, double age, int steps, int depth) */
654/* ------------------------------------------------------ 20.05.95 20.13 ---- */
655/* */
656/*  Spaltet eine RNS in zwei Species auf */
657/* */
658void splitRNS(RNS origin, double age, int steps, int depth)
659{
660    int x;
661
662    dumpRNS(origin);
663
664    /* /---------------------\ */
665    /* |  Sequenz schreiben  | */
666    /* \---------------------/ */
667
668    fprintf(seq, ">no%i\n", origin->laufNr);
669    for (x = 0; x<(origin->bases); x++) fputc(origin->base[x], seq);
670    fputc('\n', seq);
671
672    if (steps)                                             /* Species splitten! */
673    {
674        double gcDruck_val      = helixGcDruck->val,       /* Frand-Werte merken */
675               pairPart_val     = pairPart->val,
676               mutationRate_val = mutationRate->val,
677               splitRate_val    = splitRate->val;
678
679/*        indent(depth); */
680        fprintf(topo, "(no%i:%f,\n", origin->laufNr, age);
681
682        {
683            RNS left = dupRNS(origin);                     /* linker Sohn */
684
685            mutateRNS(left, steps, depth);
686            freeRNS(left);
687        }
688
689        fputs(",\n", topo);
690
691        helixGcDruck->val = gcDruck_val;                   /* Frand-Werte wiederherstellen */
692        pairPart->val     = pairPart_val;
693        mutationRate->val = mutationRate_val;
694        splitRate->val    = splitRate_val;
695
696        {
697            RNS right = dupRNS(origin);                    /* rechter Sohn */
698
699            mutateRNS(right, steps, depth);
700            freeRNS(right);
701        }
702
703        fputc(')', topo);
704    }
705    else                                                   /* Baumspitze */
706    {
707        if      (depth>maxDepth) maxDepth = depth;
708        else if (depth<minDepth) minDepth = depth;
709
710/*        indent(depth); */
711        fprintf(topo, "no%i:%f", origin->laufNr, age);
712    }
713
714    if (age==0.0) dumpRNS(NULL);
715}
716
717
Note: See TracBrowser for help on using the repository browser.