File: genetics2.ll

package info (click to toggle)
lifelines 3.0.50-2
  • links: PTS
  • area: main
  • in suites: etch-m68k
  • size: 11,140 kB
  • ctags: 6,517
  • sloc: ansic: 57,468; xml: 8,014; sh: 4,255; makefile: 848; yacc: 601; perl: 170; sed: 16
file content (113 lines) | stat: -rw-r--r-- 2,962 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/*
 * @progname       genetics2.ll
 * @version        1 of 1995-10-05
 * @author         Alexander Ottl (ottl@informatik.uni-muenchen.de)
 * @category       
 * @output         Text
 * @description

This LifeLines report program computes the degree of blood relatedness
between any two people in a database.

Genetic distance d(A,B) is defined recursively by:
  d(A,A) = 1
  d(A,B) = d(B,A)
  d(A,B) = d(F(A),B) / 2 + d(M(A),B) / 2
with F(A) and M(A) being the father and mother of A.

The recursive procedure computedist() follows that definition.
That's the beauty of recursion.

        by Alexander Ottl (ottl@informatik.uni-muenchen.de)
        Version 1 (5 Oct 1995)

*/

global(R0)
global(R1)

proc main()
{
    getindimsg(A, "First person:")
    getindimsg(B, "Second person:")
    call computedist(A, B)
    print("\nExpected degree of genetic overlap: ",
        d(R0), "/", d(R1), "\n")
}

/* BOOL ancestor( INDI, INDI ) */
func ancestor(A, B)
{
    if(not(strcmp(key(A),key(B)))) {
        return(1)
    }
    families(A, Fam, Spo, Num1) {
        children(Fam, Chl, Num2) {
            if(ancestor(Chl, B)) {
                return(1)
            }
        }
    }
    return(0)
}

/* Actually this should be a function returning a rational number.
   I might use a list, but I chose to use two global variables
   R0 and R1 for the numerator and denominator */
/* VOID computedist( INDI, INDI ) */
proc computedist(A, B)
{
    /* Recursion must terminate some time.
       One's distance to himself is 1/1 */
    if(not(strcmp(key(A),key(B)))) {
        set(R0,1)
        set(R1,1)
    }
    /* If there is a direct line from A down to B, we must work our way
       upwards from B. There must of course then be no line
       from B down to A, but no one is his own ancestor, right? */
    elsif(ancestor(A, B)) {
        /* print("Common ancestor: ", name(A), "\n") */
        call computedist(B, A)
    }
    /* Now we try to work our way upwards through the parents */
    else {
        set(R0,0)
        set(R1,1)
        if(F,father(A)) {
            call computedist(F, B)
            /* Result by half */
            set(R1, mul(2, R1))
        }
        if(M,mother(A)) {
            /* Save previous result */
            set(Res0, R0)
            set(Res1, R1)
            call computedist(M, B)
            /* Result by half */
            set(R1, mul(2, R1))
            /* Adding up with previous result */
            set(common, mul(R1, Res1))
            set(R0, add(mul(R0, Res1), mul(R1, Res0)))
            set(R1, common)
            call normalize()
        }
    }
}

/* This is not an all-purpose normalizing function.
   We expect the denominator R1 to be a power of 2 and
   to be greater than the numerator R0. */
/* VOID normalize(VOID) */
proc normalize()
{
    if(R0) {
        while(not(mod(R0,2))) {
            set(R0, div(R0,2))
            set(R1, div(R1,2))
        }
    }
    else {
        set(R1,1)
    }
}