How many fume cupboards are needed? -- Scheme version

I am going through Principles of Statistics in order to build a more respectable statistical knowledge. When I got to problem 2.6 I though it was computational heavy for such a lazy person such as I. Apparently I am not alone in that thinking.

The result is that I ended building a scheme version of the code
found in the above page. It was a very interesting exercise. You
can see the problem and the code below:

#lang racket
;; In a certain survey of the work of chemical research workers, it was
;; found, on the basis of extensive data, that on average each man
;; required no fume cupboard for 60 per cent of his time, one cupboard
;; for 30 per cent and two cupboards for 10 per cent; three or more were
;; never required. If a group of four chemists worked independently of
;; one another, how many fume cupboards should be availabe in order to
;; provode adequate facilities for at least 95 per cent of the time?
(require "cartesian-product.rkt"
         rackunit
         rackunit/text-ui)

(define probability-of-cupboards
  #hash((0 . 0.6)
        (1 . 0.3)
        (2 . 0.1)))


;; how-many-cupboards-for-% : integer number hash -> number
;; given a minimum % and a table of probabilities, find the number of
;; cupboards that will be adequated for the number of people given.
(define (how-many-cupboards-for-% number-of-people
                                  minimum-%
                                  table-of-probabilities)
  (local [(define possibilities
            (sort (hash-keys table-of-probabilities) <))
          (define (accumulate-trials-probabilities
                   trials
                   accumulated-probabilities)
            (if (empty? trials)
                accumulated-probabilities
                (accumulate-trials-probabilities
                 (rest trials)
                 (update-or-insert-probability
                  accumulated-probabilities
                  (foldl (λ (x y)
                            (+ x y))
                         0 (first trials))
                  (foldl (λ (trial-event probability-of-trial)
                            (* (hash-ref table-of-probabilities
                                         trial-event)
                               probability-of-trial))
                         1.0 (first trials))))))]
         (probability-table->result-with-%-greater-than 
          (accumulate-trials-probabilities
           (cartesian-product (make-list number-of-people possibilities)) (hash))
          minimum-%)))

;; probability-table->result-with-%-greater-than : hash number -> number or false
;; takes a probability table with the accumulated results, adds up then
;; in sequence until it surpasses the threshold. False if there is no it never
;; surpasses the threshold. 
(define (probability-table->result-with-%-greater-than table minimum-%)
  (define (accumulate-result list-of-possibilities
                             acc-probability
                             (last-probability #f))
    (cond ((empty? list-of-possibilities)
           (if (< acc-probability minimum-%) #f last-probability))
          ((> acc-probability minimum-%) last-probability)
          (else (accumulate-result (rest list-of-possibilities)
                                   (+ (hash-ref table
                                                (first list-of-possibilities))
                                      acc-probability)
                                   (first list-of-possibilities)))))
  (accumulate-result (sort (hash-keys table) <) 0))

;; update-or-insert-probability : hash integer number -> hash
(define (update-or-insert-probability table
                                      cupboard-number
                                      probability)
  (hash-update table
               cupboard-number
               (λ (old-probability)
                  (+ old-probability probability))
               0))

(define-test-suite cupboards
  (check-equal? (how-many-cupboards-for-% 4 0.95 probability-of-cupboards) 4)

  (check-equal? (probability-table->result-with-%-greater-than 
                 #hash((0 . 0.1296)
                       (1 . 0.2592)
                       (2 . 0.2808)
                       (3 . 0.1944)
                       (4 . 0.094)
                       (5 . 0.0324)) 0.94)
                4)
  (check-equal? (probability-table->result-with-%-greater-than 
                 #hash() 0.0)
                #f)
  (check-equal? (probability-table->result-with-%-greater-than 
                 #hash((0 . 0.4)
                       (1 . 0.2)) 0.7)
                #f))

(run-tests cupboards)

Serendipity, languages, YACC and robots!

While I was playing around with a particular problem this week, a
small discovery led me to a major design fix for the solution I had
developed. Here I will tell you this story because it is not only
about coding, but about inspiration, and its many sources.

The problem was the following:

You will have to trace the steps of a robot to determine its final
commands in the following language:

L - turn 90 degrees left
R - turn 90 degrees right
M - move forward
T - transport to a given location

The commands will be given as such:


10 10 # board size
2 5 N # initial location and the direction the robot is facing
LLRRMMMRLRMMM # series of moves
T 1 3 # transport to position x=1 y=3
LLRRMMRMMRM # another series of moves


Summing up, it is a small textual logo.


My approach to the problem was basically to draft the structures I
would need, in a semi object oriented approach (old-habits). Much like
the following(the first version is in portuguese):

;; tabuleiro é uma estrutura contendo um par de inteiros
(define-struct tabuleiro (x y) #:transparent)

(define-struct robo (x y direção) #:transparent)

;; constrói-tabuleiro : string -> tabuleiro
(define (constrói-tabuleiro dados)
  (let ([dados-do-tabuleiro (string-tokenize dados)])
    (make-tabuleiro (string->number (first dados-do-tabuleiro))
                    (string->number (second dados-do-tabuleiro)))))

;; constrói-tabuleiro : string -> tabuleiro
(define (constrói-robo dados)
  (let ([dados-do-robo (string-tokenize dados)])
    (make-robo (string->number (first dados-do-robo))
               (string->number (second dados-do-robo))
               (string-ref (third dados-do-robo) 0))))

All was going pretty smooth. I had defined some structures, later
the functions that mapped to the commands in the mini logo:

(define DIREÇÕES-POSSIVEIS
  (list '(#\N (#\W #\E))
        '(#\W (#\S #\N))
        '(#\S (#\E #\W))
        '(#\E (#\N #\S))))

(define RESULTADO-DO-PASSO-A-FRENTE
  (list (list #\N (list (λ (n) n) add1))
        (list #\W (list sub1 (λ (n) n)))
        (list #\S (list (λ (n) n) sub1))
        (list #\E (list add1 (λ (n) n)))))

;; muda-direção : robo char -> robo
(define (muda-direção o-robo virar-para)
  (let ([possíveis-direções
         (second (findf (λ (uma-possibilidade)
                          (equal? (first uma-possibilidade)
                                  (robo-direção o-robo)))
                        DIREÇÕES-POSSIVEIS))])
    (struct-copy robo
                 o-robo
                 (direção (if (equal? virar-para #\L)
                              (first possíveis-direções)
                              (second possíveis-direções))))))


;; muda-posição : robo number number -> robo
(define (muda-posição o-robo x y)
  (struct-copy robo
               o-robo
               (x x)
               (y y)))

;; passo-a-frente : robo -> robo
(define (passo-a-frente o-robo)
  (let ([possíveis-funções
         (second (findf (λ (um-resultado-de-uma-posição)
                          (equal? (first um-resultado-de-uma-posição)
                                  (robo-direção o-robo)))
                        RESULTADO-DO-PASSO-A-FRENTE))])
    (muda-posição o-robo
                  ((first possíveis-funções) (robo-x o-robo))
                  ((second possíveis-funções) (robo-y
  o-robo)))))

That also went ok. But the problem was parsing the language. I did
my own ad-hoc parser, but as I was building it, the felling of
suspicion was growing and my faith in my approach diminishing. Take
a look at the final product:

  
;; ler-movimentos-do-robo : input-port robo tabuleiro -> robo
;; le e modifica a posição do robo a partir de uma lista de movimentos
(define (ler-movimentos-do-robo arquivo
                                (o-robo #f)
                                (o-tabuleiro #f))
  (let* ([linha-de-comandos (read-line arquivo 'any)])
    (if (eof-object? linha-de-comandos)
        o-robo
        (let ([elementos-da-linha (string-tokenize linha-de-comandos)])
          (case (length elementos-da-linha)
            [(3)
             (if (equal? (first elementos-da-linha) "T")
                 (ler-movimentos-do-robo arquivo
                                         (muda-posição
                                          o-robo
                                          (string->number (second
                                                           elementos-da-linha))
                                          (string->number (third
                                                           elementos-da-linha)))
                                         o-tabuleiro)
                 (ler-movimentos-do-robo arquivo
                                         (constrói-robo linha-de-comandos)
                                         o-tabuleiro))]
            [(2) (ler-movimentos-do-robo arquivo
                                         o-robo
                                         (constrói-tabuleiro linha-de-comandos))]
            [else
             (ler-movimentos-do-robo arquivo
                                     ; cria um novo robo pra cada comando da linha
                                     (string-fold
                                      (lambda (comando velho-robo)
                                        (if (equal? comando #\M)
                                            (passo-a-frente velho-robo)
                                            (muda-direção velho-robo comando)))
                                      o-robo
                                      linha-de-comandos)
                                     o-tabuleiro)])))))

Big, butt-ugly, hard to read and maintain. Something was wrong
and I wasn't really sure what. Then serendipity struck me in the
form of
this message. Lex
and YACC, of course! It took me a couple of hours to learn Racket's
syntax for it, but I ended up with a lexer and a parser for the
mini logo, in a version that was much more friendly and elegant.

(define-tokens value-tokens (NUMBER
                             DIRECTION
                             FUNCTION
                             COMMANDS))

(define-empty-tokens op-tokens (EOF NEWLINE))

(define-lex-abbrevs
  (teleport     #\T)
  (turn-left    #\L)
  (turn-right   #\R)
  (move-forward #\M)
  (directions (:or #\N
                   #\S
                   #\E
                   #\W))
  (digit        (:/ #\0 #\9)))

(define mech-lexer
  (lexer
   [(eof) 'EOF]
   [#\newline (token-NEWLINE)]
   ;; recursively call the lexer on the remaining input after a tab or space.    
   ;; Returning the result of that operation.  
   ;; This effectively skips all whitespace.
   [#\space (mech-lexer input-port)]
   [directions (token-DIRECTION (string-ref lexeme 0))]
   [(:+ (:or turn-right turn-left move-forward))
    (token-COMMANDS
     (map (λ (command)
             (get-consequences COMMAND-LOOKUP-TABLE command))
          (string->list lexeme)))]
   [teleport (token-FUNCTION (string-ref lexeme 0))]
   [(:+ digit) (token-NUMBER (string->number lexeme))]))

;; mech-parser : (X -> token) board mech -> mech
(define (mech-parser gen
                     (board #f)
                     (mech #f))
  ((parser
    (start start)
    (end EOF NEWLINE)
    (tokens value-tokens op-tokens)
    (error (λ (tok-ok? tok-name tok-value)
              (error 'lexer
                     (format
                      (if (false? tok-ok?)
                          "the token ~a with value ~a was invalid."
                          "unknow error with token ~a with value ~a")
                      tok-name
                      tok-value))))
    
    (grammar
     (start [() mech] ;; end our parsing, return the mech structure
            [(exp) $1])
     
     (exp [(NUMBER NUMBER)
           (mech-parser gen
                        (make-board $1 $2)
                        mech)]
          
          [(NUMBER NUMBER DIRECTION)
           (mech-parser gen
                        board
                        (make-mech $1 $2 $3))]
          
          [(FUNCTION NUMBER NUMBER)
           (with-handlers ((string? (λ (message)
                                       (display message)
                                       (mech-parser gen
                                                    board
                                                    mech))))
                          (mech-parser gen
                                       board
                                       (jump board mech $2 $3)))]
          
          [(COMMANDS)
           (mech-parser gen
                        board
                        (foldl
                         (λ (command old-mech)
                            (with-handlers ((string? (λ (message)
                                                        (display message)
                                                        old-mech)))
                                           (command board old-mech)))
                         mech $1))])))
   gen))

(define (process-mech ip)
  (mech-parser (λ () (mech-lexer ip))))

If you want to compare the full versions, take a look at the
github repository. The
version with the cleaner Lexer/Parser are called 'mech' instead of
'robo'.
Until next time!

Bitvector genealogy update

In the past few days I have been fiddling with several possibilities
for the bitvector genealogy problem. I was aiming to solve the problem
with the root bitvector that my first solution had.

I have found different sources on the web and different ways of
dealing with the problem, but none matched the accuracy of the first
solution, although it did surpass it in speed.

Here you will find the latest racket code with a several different takes on
the solution. The main things to look for are the find-genealogy-*
functions.

For now I am done with this until I can come up with a different
strategy for approaching the problem.

Bitvector Genealogy

Introduction

Recently I have come across something that reminded me of a challenge
I tried some time ago, ITAs bitvector genealogy
challenge(1). Back in the day when I was first
learning scheme, I have tried to use this to build up my skills with
the help of a great tutorial(3).
Needless to say that I failed at that, both for my lack of scheme and
common lisp competence but also for my lack of aptitude as a
programmer.

The something was a thread in Hacker News where ITA was mentioned. I
became curious to see how much I had progressed since my last attempt,
so I grabbed DrRacket, Emacs and jumped into the fray one more time.

Again, most of the original hard work was done by Patrick May. What I
am showing here is a take on his work, translated to scheme
and modified to use different data structures and different points of
views on the same fundamental algorithm.

Ok, this is the description of the problem:

The BitVectors are an ancient and immortal race of 10,000, each with
a 10,000 bit genome. The race evolved from a single individual by
the following process: 9,999 times a BitVector chosen at random from
amongst the population was cloned using an error-prone process that
replicates each bit of the genome with 80% fidelity (any particular
bit is flipped from parent to child 20% of the time, independent of
other bits).

Write a program to guess the reproductive history of BitVectors from
their genetic material. The randomly-ordered file
bitvectors-genes.data.gz contains a 10,000 bit line for each
individual. Your program's output should be, for each input line,
the 0-based line number of that individual's parent, or -1 if it is
the progenitor. Balance performance against probability of mistakes
as you see fit.
(1)

The next step is really the crux of the problem. This is where a bit
of mathematics comes into play

The next step in developing a solution is to determine whether or not
a parent-child relationship exists between two bit vectors. According
to the problem definition, every bit vector other than the original is
created by randomly mutating the elements of an existing bit
vector. The probability of a particular bit changing is 20%. That
means that, on average, the children of a bit vector should be 80%
identical to their parent and the grandchildren should be 64%
similar. Of course, the random element means that some will have more
than 80% similarity and others less. Some children may be more
similar to each other than either is to their parent. The best result
any automated solution will be able to achieve is the most probable
relationship tree.

Assuming that a fair random number generator was used, the set
containing the number of bits different between a parent and each of
its children should follow a binomial distribution. The probability of
two vectors having a parent child relationship is related to how close
the difference between the two is to the expected value (20%, in this
case). In a binomial distribution (which approximates a continuous
normal distribution), 99.73% of the bit difference values will be
within three standard deviations of the expected value. The standard
deviation is computed by:

Body

First, I will define some globals that will be useful throughout the
source code, and the utilities to load the bitvectors data into a
vector of numbers:

(define ROOT-SYMBOL -1)

(define MUTATION-RATE 0.2)

;; http://mathworld.wolfram.com/BinomialDistribution.html
(define DEVIATION 3) ;; 3 = 99.73 percent


(define DNA-LENGTH 0)

(define (get-dna)
  DNA-LENGTH)

(define (set-dna! number)
  (set! DNA-LENGTH number))

;; filename->bitvector : input-port -> (vectorof binary-number)
;; read everything from a given file and returns a 
;; vector containing binary numbers. 
(define (filename->bitvector file-path (radix 2))
  (local [(define (read-accumulative bitvector)
            (let ([line-read (read-line file-path 'any)])
              (if (eof-object? line-read)
                  bitvector
                  (let ([possible-number (string->number line-read radix)])
                    (read-accumulative
                     (if (number? possible-number)
                         (vector-append bitvector (vector possible-number))
                         bitvector))))))]
    (cond
      [(equal? "" file-path)]
      [else
       (read-accumulative (vector))])))

The next part is really the crux of the problem. This is where a bit
of mathematics comes into play (2)
and the most complex part of the problem in my view:

The next step in developing a solution is to determine whether or not
a parent-child relationship exists between two bit vectors. According
to the problem definition, every bit vector other than the original is
created by randomly mutating the elements of an existing bit
vector. The probability of a particular bit changing is 20%. That
means that, on average, the children of a bit vector should be 80%
identical to their parent and the grandchildren should be 64%
similar. Of course, the random element means that some will have more
than 80% similarity and others less. Some children may be more
similar to each other than either is to their parent. The best result
any automated solution will be able to achieve is the most probable
relationship tree.

Assuming that a fair random number generator was used, the set
containing the number of bits different between a parent and each of
its children should follow a binomial distribution. The probability of
two vectors having a parent child relationship is related to how close
the difference between the two is to the expected value (20%, in this
case). In a binomial distribution (which approximates a continuous
normal distribution), 99.73% of the bit difference values will be
within three standard deviations of the expected value.
(3)

Taking that, we can establish the existence of a relationship (with
some degree of confidence) with the following code:

;;  binary-difference : binary-number binary-number -> integer
;; returns the number of bits that are different between 2 binary numbers.
(define (binary-difference bits1 bits2)
  (let ([differences 0])
    (for ((i (in-string (number->string (bitwise-xor bits1 bits2) 2))))
      (when (equal? i #\1)
        (set! differences (add1 differences))))
    differences))

;; binomial-standard-deviation : number -> number
;; Compute the expected binomial standard deviation for N events of probability P.
(define (binomial-standard-deviation n p)
  (sqrt (* n p (- 1 p))))

;; has-kinship? : integer integer -> boolean
;; verify if 2 individuals share a kinship. 
;; The relationship is assumed if the bit difference between the two is
;; less than DEVIATIONS standard deviations from the expected average for
;; the MUTATION-RATE.
(define (has-kinship? individual-1
                      individual-2
                      (mutation-rate MUTATION-RATE)
                      (deviations DEVIATION)
                      (dna-length (get-dna)))
  (let* ([expected-difference (* dna-length mutation-rate)]
         [standard-deviation
          (binomial-standard-deviation dna-length
                                       mutation-rate)]
         [max-difference (+ expected-difference
                            (* deviations standard-deviation))])
    (< (binary-difference individual-1
                          individual-2)
       max-difference)))

The code above shows in some ways where I deviated from Patrick's
example. I used a number instead of a bit array (which does not exist in
Racket, the language I'm using, as far as I know). In the following
code, this difference deepens, with the construction of a hash map
instead of a vector to get the relationship matrix:

;; create-relationship-hash : (vectorof number) -> hash
;; creates a hash mapping an individual to it's list of relationships.
;; Both parents and children are mapped to an individual.
(define (make-relationship-hash bitvector)
  (let ([relationship-hash (make-hash)]
        [dna-range (in-range 0 (vector-length bitvector))])
    (for* ((individual-1-index dna-range)
           (individual-2-index dna-range))
      (begin (unless (hash-has-key? relationship-hash individual-1-index)
               (hash-set! relationship-hash individual-1-index empty))
             (let ([individual-1 (vector-ref bitvector
                                             individual-1-index)]
                   [individual-2 (vector-ref bitvector
                                             individual-2-index)])
               (when (and (has-kinship? individual-1
                                        individual-2)
                          (not (= individual-1 individual-2)))
                 (hash-update! relationship-hash
                               individual-1-index
                               (λ (relationships)
                                 (cons individual-2-index relationships)))))))
    relationship-hash))

This almost closes the issue, but I must figure out the directions of
the relationships(who is who's father and son). To establish that I
should have to:

• Find all the bit vectors with only a single parent-child relationship
  
• Assume that those vectors are the child nodes
  
• Record that information
  
• Remove the child node from the parent's list of relationships
  
• Repeat
  
  

(3)

That is accomplished in the longest and most complex code in the solution:

;; make-kinship-hash : hash -> hash
;; creates a hash mapping an individual only to it's parent.
(define (make-kinship-hash relationship-hash-initial)
  (local [(define (make-kinships kinship-hash
                                 relationship-hash)
            (local [(define (process-possible-kinship individual-index
                                                      relationships-indexes)
                      (let ([number-of-relationships
                             (length relationships-indexes)])
                        (unless (> number-of-relationships 1)
                          (begin
                            (hash-remove! relationship-hash
                                          individual-index)
                            (if (= number-of-relationships 1)
                                ; leaf node
                                (begin
                                  (hash-set! kinship-hash
                                             individual-index
                                             (first relationships-indexes))
                                  (hash-update! relationship-hash
                                                (first relationships-indexes)
                                                (λ (relationships)
                                                  (remove individual-index
                                                          relationships))))
                                ; root node
                                (begin (hash-remove! relationship-hash
                                                     individual-index)
                                       (hash-set! kinship-hash
                                                  individual-index
                                                  ROOT-SYMBOL)))))))]
              (if (false? (hash-iterate-first relationship-hash))
                  kinship-hash
                  (begin
                    (hash-for-each relationship-hash
                                   process-possible-kinship)
                    (make-kinships kinship-hash
                                   relationship-hash)))))]
    (make-kinships (make-hash) relationship-hash-initial)))

And that is the end of the solution! I am still facing the same issues
as Patrick, specially in regard to the root individual, but
reimplementing his solution and trying to achieve a better result by
the use of different data structures led me to understand the problem
more deeply, and it was loads of fun!

PS: There is more code in
this
repository, specially regarding testing and comparison.

Bibliography

1

ITA.

Bitvector genealogy.

URL
http://www.itasoftware.com/careers/work-at-ita/hiring-puzzles.html.

2

MathWorld.

Binomial distribution.

URL http://mathworld.wolfram.com/BinomialDistribution.html.

3

Patrick May.

Bitvector genealogy.

URL http://www.softwarematters.org/bitvector-genealogy.html.