From: Alex Martelli <aleaxit@...>
Date: 2003-09-15T01:44:09+09:00
Subject: performance and style advice requested

I'm trying to learn some Ruby, so I want to write some Ruby code, starting 
with my usual favourite fun field -- combinatorial arithmetic related to 
the game of contract bridge.  As a warm-up, I begin with the "probability
of pointcount".  Problem definition: a hand is a set of 13 cards out of 52.
16 of those cards (Aces, Kings, Queens, Jacks) are known as "honors" and
are assigned a "pointcount" (4 for each Ace, 3 for each King, 2 for each
Queen, 1 for each Jack).  Consider all possible hands and determine the
exact probability of each possible point-count.

Basically, this boils down to considering each of the 2**16 possible
sets of honors -- each has an easily computed total pointcount -- and
for each of them computing in how many different ways it can occur --
which is C(36, 13-numberofhonors) -- the ways of choosing out of the
36 non-honors in the deck the 13-numberofhonors non-honors in this
hand.  (That's with C(x,y) extended to be 0 for y<0, of course).

So I sit down and code it in the simplest way -- in Python first for
practice, since that's what I'm familiar with, then in Ruby.  Fine
both ways, except that (on my laptop, which is all I have at hand as
I'm on a business trip -- with Python 2.3 and Ruby 1.8 and Win/XP)
the Python version runs in just above 4 seconds while the Ruby one
takes over 9.  Hmmm, probably my lack of Ruby instincts, of course,
but anyway -- just as of course -- this needs to be optimized by
quite a bit (not for its own sake, but because later I'll extend it
to compute much more interesting stuff, such as correlations between
points, controls and shape, etc, etc) so who cares about performance
for the "coded in the simplest way" versions (oh, I should clarify:
simplest _except_ that I have of course memoized C(), computing
combinations aka binomial coefficients, and the factorial it calls).

So I start optimizing and squeezing every cycle I can -- and I get
down to about 0.65 seconds for Python, 2.3 seconds for Ruby (best
case, for each of them).  Eep -- the ratio keeps increasing, so my
lack of sound Ruby instincts must be really hurting.  Oh well, say
I, let's try "ruby -rprofile".  EEK!  After a couple of times where
I was convinced it had just seized up I decided to let it run all
the way -- and it took almost 1000 seconds.  Is a slowdown of over
400 times normal for Ruby profiling, or is something weird in my
installation...?  (It's the one-exe native Ruby 1.8 install for Win
which comes with a lot of goodies such as SciTE, Programming Ruby
in HTML-Help format, etc).

Anyway, the profiler's output isn't illuminating to me at all, so
that's when I decide to turn to the famous Ruby community -- I'm
sure you'll find lots to critique in my Ruby program, particularly
with an eye to performance but not necessarily just that (I _am_
quite aware that I may be guilty of "coding Python in Ruby" -- and
this may produce suboptimal performance _and_ other issues).

So, first, for reference, here's the optimized Python program:


#
# Memoized implementation of factorial
#
def fact(x, _fact_memo = {0:1, 1:1, 2:2, 3:6, 4:24, 5:120, 6:720, 7:5040} ):
    " factorial of x (0 if x<0) "
    try: return _fact_memo[x]
    except KeyError:
        if x < 0: return 0
        return _fact_memo.setdefault(x, x * fact(x-1))

#
# Memoized implementation of C(x, y) [binomial coefficient, aka comb]
#
def comb(x, y, _comb_memo = {}):
    " combinations of x things y at a time (0 if x<0, y<0, or y>x) "
    try: return _comb_memo[x, y]
    except KeyError:
        if x<0 or y<0 or y>x: return 0
        return _comb_memo.setdefault((x, y), fact(x) / (fact(y)*fact(x-y)))

#
# we represent a "honor" by a its points-worth
#
honors = [ p for p in range(1,5) for s in range(4) ]

#
# we represent a "honors set" by a pair (points-worth, total-#-of-cards)
#
def honor_subsets(lis, fromindex=0):
    """ recursive generator of all sets from a list and starting point;
        each set is represented by a honor-set as above, while the list
        must hold points-worths representing "honors" again as above.
    """
    ph = lis[fromindex]
    if fromindex == len(lis)-1:
        yield 0, 0
        yield ph, 1
        return
    for pw, le in honor_subsets(lis, fromindex+1):
        yield pw, le
        yield ph+pw, 1+le

def main():
    " main computation "
    # hist: histogram  points -> number of occurrences
    hist = {}
    # totcom: total number of occurrences
    totcom = 0L
    # range over all possible honorsets
    for hs in honor_subsets(honors):
        # compute number of honors and pointcount for this honorset
        pt, nh = hs
        # compute number of possible occurrences of this honorset
        nc = comb(36, 13-nh)
        # update histogram and total
        hist[pt] = hist.get(pt, 0) + nc
        totcom += nc
    # print total and eyeball-check it
    print totcom, comb(52, 13)

    # sort histogram by decreasing frequency
    aux = [ (nc, pt) for pt, nc in hist.iteritems() ]
    aux.sort()
    aux.reverse()

    # compute frequencies as relative percentages
    divis = 100.0 / totcom

    # give top 10 possibilities
    for nc, pt in aux[:10]:
        print "%2d %d (%.2f)" % (pt, nc, nc * divis)


# run the computation
import time
start = time.time()
main()
stend = time.time()
print stend-start


and here is the corresponding optimized Ruby program:


#
# factorial, memoized via an Array
#
$_fact_memo=[1,1,2,6,24,120,720]
def fact(x)
    # factorial of x (0 if x<0)
    result = $_fact_memo[x]
    if !result
        return 0 if x<0
        result = $_fact_memo[x] = x*fact(x-1)
    end
    return result
end

#
# binomial factor (aka comb), memoized via a Hash
#
$_comb_memo={}
def comb(x, y)
    # combinations of x things y at a time (0 if x<0, y<0, or y>x)
    result = $_comb_memo[[x, y]]
    if !result
        return 0 if x<0 or y<0 or y>x
        result = $_comb_memo[[x, y]] = fact(x)/(fact(y)*fact(x-y))
    end
    return result
end

#
# add to Array a recursive iterator over honor-subsets
#
# a honor is represented by its point-value (1..4); a honor-subset is
# a pair [total-pointcount, number-of-honors]
#
class Array
    def each_subset(from=0)
        if from==size-1:
            #
            # end of array, so, return 2 only possible subsets remaining:
            # empty set, or singleton set with the last element in it
            #
            yield [0, 0]
            yield [self[from], 1]
        else
            #
            # 2+ elements remaining, so, loop iteratively, extracting
            # the current element and getting all possible subsets of
            # all later ones -- for each of those, return 2 possible
            # subsets, i.e., either without or with the current element
            #
            thispoints = self[from]
            each_subset(from+1) do |subpoints, sublength|
                yield subpoints, sublength
                yield thispoints + subpoints, 1 + sublength
            end
        end
    end
end


def main()
    # main computation 
    # hist: histogram (Hash) points -> number of occurrences
    hist = Hash.new(0)
    # totcon: total number of occurrences
    totcom = 0
    # range over all possible honor-sets
    honors = Array(1..4) * 4
    honors.each_subset do |pt, nh|
        # receive number of honors and pointcount for this honorset,
        # then compute number of possible occurrences of this honorset
        nc = comb(36, 13-nh)
        # update histogram and total
        hist[pt] += nc
        totcom += nc
    end
    # print total and eyeball-check it
    print totcom, ' ', comb(52, 13), "\n"

    # sort histogram by decreasing frequency
    aux = hist.sort {|a,b| b[1]<=>a[1]}

    # compute frequencies as relative percentages
    divis = totcom / 100.0

    # give top 10 possibilities
    for pt, nc in aux[0,10]
        printf "%2d %d (%.2f)\n" , pt, nc, nc/divis
    end
end

# run the computation
start = Time.now
main
stend = Time.now
puts stend-start


Finally, FWIW, here's the start of the profiler's output...:


C:\Python23>ruby -rprofile alco2.rb
635013559600 635013559600
10 59723754816 (9.41)
 9 59413313872 (9.36)
11 56799933520 (8.94)
 8 56466608128 (8.89)
 7 50979441968 (8.03)
12 50971682080 (8.03)
13 43906944752 (6.91)
 6 41619399184 (6.55)
14 36153374224 (5.69)
 5 32933031040 (5.19)
948.294
  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
 34.90   317.57    317.57       16 19848.44 909790.50  Array#each_subset
 13.12   437.00    119.43    65537     1.82     6.56  Object#comb
 10.56   533.11     96.11   131073     0.73     2.54  Hash#[]
  8.20   607.70     74.59    65385     1.14     1.82  Array#eql?
  8.11   681.49     73.79    65552     1.13     1.79  Array#hash
  6.17   737.65     56.16   161735     0.35     0.35  Fixnum#+
  4.87   781.99     44.34   130770     0.34     0.34  Numeric#eql?
  4.80   825.63     43.65   131104     0.33     0.33  Kernel.hash
  4.14   863.31     37.68   100426     0.38     0.38  Bignum#+
  2.59   886.84     23.53    65551     0.36     0.36  Hash#[]=
  2.48   909.38     22.54    65613     0.34     0.34  Fixnum#-
  0.02   909.57      0.19       91     2.09    38.57  Object#fact
  0.02   909.75      0.18        1   180.00   180.00  
Profiler__.start_profile
  0.02   909.92      0.17        1   170.00   190.00  Hash#sort
  0.01   910.02      0.10      350     0.29     0.29  Fixnum#<
  0.01   910.10      0.08      193     0.41     0.41  Hash#default
  0.01   910.16      0.06       59     1.02     1.19  Fixnum#*

I've truncated starting from where the "%time" goes to 0.00.


Thanks for any feedback!


Alex