/* vi: set shiftwidth=4 :
genetics2 - a LifeLines report program to calculate degree of relatedness
by Alexander Ottl (ottl@informatik.uni-muenchen.de)
Version 1 (5 Oct 1995)
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.
*/
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)
}
}