# Mendelian inheritance

I started out with the Rosalind problems a while back. Turns out it's a great way to learn a new language. I'm still trying to learn Ruby, and while I haven't managed to use regular expressions as much as I wanted to (learning those being one of the reasons I started learning Ruby in the first place), I put together a repository of whatever challenges I have solved so far.

I learned how easy it is to benchmark with Ruby, so I put that to the test as well. While the resulting solution is good, it's quite possible there's a better one out there.

### Problem: IPRB

Probability is the mathematical study of randomly occurring phenomena. We will model such a phenomenon with a random variable, which is simply a variable that can take a number of different distinct outcomes depending on the result of an underlying random process.

For example, say that we have a bag containing 3 red balls and 2 blue balls. If we let $X$ represent the random variable corresponding to the color of a drawn ball, then the probability of each of the two outcomes is given by $Pr[X=red]=35$ and $Pr[X=blue]=25$.

Random variables can be combined to yield new random variables. Returning to the ball example, let $Y$ model the color of a second ball drawn from the bag (without replacing the first ball). The probability of $Y$ being red depends on whether the first ball was red or blue. To represent all outcomes of $X$ and $Y$, we therefore use a probability tree diagram. This branching diagram represents all possible individual probabilities for $X$ and $Y$, with outcomes at the endpoints ("leaves") of the tree.

### Given:

Three positive integers $k, m,$ and $n$, representing a population containing $k+m+n$ organisms: $k$ individuals are homozygous dominant for a factor, $m$ are heterozygous, and $n$ are homozygous recessive.

### Return:

The probability that two randomly selected mating organisms will produce an individual possessing a dominant allele (and thus displaying the dominant phenotype). Assume that any two organisms can mate.

### Sample dataset:

2 2 2


### Sample output:

0.78333


My solution solves the sample dataset and the actual dataset given:

### Dataset:

22 26 27


### Solution:

0.71775


### IPRB.rb:

def probability_dominant(k, m, n)
(k * k + k * (2 * m + 2* n - 1) + m * (0.75 * m + n - 0.75)) / ((k + m + n - 1) * (k + m + n))
end

user_input = gets.chomp

split_input = user_input.split().map { |number| number.to_i() }
puts probability_dominant(split_input[0], split_input[1], split_input[2]).round(5)


The math involved appears to be the performance bottleneck.

### IPRB_benchmark.rb:

def probability_dominant_1(k, m, n)
((k * k - k) + 2 * (k * m) + 2 * (k * n) + (0.75 * (m * m - m)) + 2 * (0.5 * m * n))/((k + m + n)*(k + m + n -1))
end

def probability_dominant_2(k, m, n)
(k * k + k * (2 * m + 2* n - 1) + m * (0.75 * m + n - 0.75)) / ((k + m + n - 1) * (k + m + n))
end

def probability_dominant_3(k, m, n)
(k * k + k * m + k + 0.25 * m * m + 0.75 * m)/(k + m + n) - (k * k + k * m - k + 0.25 * m * m - 0.25 * m) / (k + m + n - 1)
end

require 'benchmark'

n = 500000
Benchmark.bm(6) do |x|
x.report("first:")   { for i in 1..n; probability_dominant_1(n, n, n); end }
x.report("second:")   { for i in 1..n; probability_dominant_2(n, n, n); end }
x.report("third:")   { for i in 1..n; probability_dominant_3(n, n, n); end }
end


### Output:

             user     system      total        real
first:   0.594000   0.000000   0.594000 (  0.597825)
second:  0.328000   0.000000   0.328000 (  0.316578)
third:   0.500000   0.000000   0.500000 (  0.517067)


I optimized the math the best I could, but optimizing for a computer is a lot sketchier than optimizing for a human. It's quite hard to keep it fast for the machine and readable for a human at the same time. I foresee the upcoming challenges to be more and more mathematically inclined, so I better tackle this here and now.

Any and all advice about how to math in Ruby is appreciated. Nitpicks are welcome too.

### Profiling and Benchmarking

Rather than benchmarking, which just does total time, you would be better off profiling, which is like recursive benchmarking. I recommend you use the ruby-prof gem.

require 'ruby-prof'

n = 500000
puts "========================DOM 1============================="
RubyProf.start
(1..n).each { |n| probability_dominant_1(n, n, n) }
result = RubyProf.stop
printer = RubyProf::FlatPrinter.new(result)
printer.print(STDOUT)

puts "========================DOM 2============================="
RubyProf.start
(1..n).each { |n| probability_dominant_2(n, n, n) }
result = RubyProf.stop
printer = RubyProf::FlatPrinter.new(result)
printer.print(STDOUT)

puts "========================DOM 3============================="
RubyProf.start
(1..n).each { |n| probability_dominant_3(n, n, n) }
result = RubyProf.stop
printer = RubyProf::FlatPrinter.new(result)
printer.print(STDOUT)


And got the following results:

========================DOM 1=============================
Measure Mode: wall_time
Fiber ID: 22475080
Total: 4.755476
Sort by: self_time

%self      total      self      wait     child     calls  name
55.06      4.501     2.618     0.000     1.883   500000   Object#probability_dominant_1
11.69      0.556     0.556     0.000     0.000  1434466   Fixnum#*
10.14      0.482     0.482     0.000     0.000  1433660   Bignum#+
7.95      0.378     0.378     0.000     0.000  1500000   Float#*
6.60      0.314     0.314     0.000     0.000   934466   Bignum#-
5.34      4.755     0.254     0.000     4.501        1   Range#each
3.13      0.149     0.149     0.000     0.000   500000   Float#/
0.08      0.004     0.004     0.000     0.000    24252   Fixnum#+
0.00      4.755     0.000     0.000     4.755        1   Global#[No method]

* indicates recursively called methods
========================DOM 2=============================
Measure Mode: wall_time
Fiber ID: 22475080
Total: 2.474247
Sort by: self_time

%self      total      self      wait     child     calls  name
60.31      2.258     1.492     0.000     0.766   500000   Object#probability_dominant_2
10.67      0.264     0.264     0.000     0.000   952579   Bignum#+
8.73      2.474     0.216     0.000     2.258        1   Range#each
5.38      0.133     0.133     0.000     0.000   500000   Fixnum#*
5.17      0.128     0.128     0.000     0.000   500000   Float#/
4.73      0.117     0.117     0.000     0.000   500000   Float#+
4.73      0.117     0.117     0.000     0.000   500000   Float#*
0.28      0.007     0.007     0.000     0.000    31037   Fixnum#+
0.00      2.474     0.000     0.000     2.474        1   Global#[No method]

* indicates recursively called methods
========================DOM 3=============================
Measure Mode: wall_time
Fiber ID: 22475080
Total: 4.334434
Sort by: self_time

%self      total      self      wait     child     calls  name
50.85      4.074     2.204     0.000     1.870   500000   Object#probability_dominant_3
17.44      0.756     0.756     0.000     0.000  2364956   Bignum#+
16.66      0.722     0.722     0.000     0.000  3000000   Float#*
6.00      4.334     0.260     0.000     4.074        1   Range#each
5.42      0.235     0.235     0.000     0.000  1000000   Float#/
3.39      0.147     0.147     0.000     0.000   476830   Bignum#-
0.23      0.010     0.010     0.000     0.000    46340   Fixnum#+
0.00      4.334     0.000     0.000     4.334        1   Global#[No method]

* indicates recursively called methods


Foreach function, these tables show the number of calls and time spent on function call, recursively all the way down. Some portions, like the "range#each" can be ignored as part of of our testing overhead.

Looking at the results, I can spot why dominant_2 is the fastest. Integer operations (FixNum) on a processor will always be faster. Floating point (Float) is slow and floating point divisiondividing is slow, and big integers (integer numbers bigger than that hardware can natively support, typically 2^32 or 2^64) will also be slower.

• Dom1 has 2M Float, ~2.4M BigNum, and ~1.4M FixNum operations.
• Dom2 has 1.5M Float, ~1M BigNum, and ~530K FixNum operations.
• Dom3 has 4M Float, ~2.3M BigNum, and ~46K FixNum operations.

Dom2 was the fastest because it has the fewest number of operations by far.

As a weak rule of thumb, when optimizing math, try to avoid BigNum and Float when feasible. Also try to avoid things like converting numbers to strings and then back to numbers.

### Conclusion

Know that more often though, that your slowdown will usually come from lots of branching and flow control, and being inefficient (like calculating the same things over and over rather than caching them), and those are the places to optimize. If this problem were part of a larger, more complicated task, I can practically guarantee that your bottleneck would not be ANY of these functions, but something else entirely.

• I was hoping something like this would be easily available, thank you for pointing me in the right direction. With a bit of trial and error, this should help optimizing the formula even further. I agree with your assessment that it will most likely not be the bottleneck in any complicated process, but as I tried to explain in my question, the upcoming processes are quite likely math-heavy. I better learn how to deal with it and this answer provides valuable information for the process of doing so.
– Mast
Apr 21 '16 at 18:23