Poisson-distributed random numbers in Typed Racket

poisson-ffi.rkt

#lang racket/base

;; translated from an algorithm in C from
;; Hyberts SG, Takeuchi K, Wagner G (2010) Poisson-gap sampling and
;; forward maximum entropy reconstruction for enhancing the resolution
;; and sensitivity of protein NMR Data. J Am Chem Soc 132:2145–2147

;; generate random Poisson-distributed numbers as given
;; by Donald E. Knuth (1969). Seminumerical Algorithms.
;; The Art of Computer Programming, Volume 2. Addison Wesley


(require ffi/unsafe racket/cmdline)
(define srand48 (get-ffi-obj 'srand48 #f (_fun _long -> _void)))
(define drand48 (get-ffi-obj 'drand48 #f (_fun -> _double)))


(define-syntax-rule (while condition body ...)
  (let L () (when condition body ... (L))))
(define-syntax-rule (do-while condition body ...)
  (let L () (when (let () body ... condition) (L))))

;(: poisson : (Float -> Integer))
(define (poisson lmbd)
  (define L (exp (- lmbd)))
  (define k 0)
  (define p 1.)
  (do-while (p . >= . L)
    (define u (drand48))
    (set! p (* p u))
    (set! k (+ k 1)))
  (- k 1))

;(: run (Exact-Positive-Integer Fixnum Fixnum -> (Listof Integer)))
(define (run s p z)
  (define ld (/ (exact->inexact z) p))
  (define adj (* 2.0 (- ld 1)))
  (define i 0) (define j 0)
  (define k 0) (define n 0)
  (define v (make-vector z))
  (srand48 s)
  (do-while (not (= n p))
    (set! i 0)
    (set! n 0)
    (while (< i z)
      (vector-set! v n i)
      (set! i (+ i 1))
      (set! k (poisson (* adj (sin (* (/ (+ i 0.5) (+ z 1.)) 1.5707963268)))))
      (set! i (+ i k))
      (set! n (+ n 1)))
    (when (> n p) (set! adj (* adj 1.02)))
    (when (< n p) (set! adj (/ adj 1.02))))
  (for/list ([j (in-range p)]) (vector-ref v j)))


(command-line
 #:args
 (s p z)
 (void (run (string->number s) (string->number z) (string->number p))))

poisson-fun.rkt

#lang typed/racket
(require racket/cmdline)

;; translated from an algorithm in C from
;; Hyberts SG, Takeuchi K, Wagner G (2010) Poisson-gap sampling and
;; forward maximum entropy reconstruction for enhancing the resolution
;; and sensitivity of protein NMR Data. J Am Chem Soc 132:2145–2147

;; generate random Poisson-distributed numbers as given
;; by Donald E. Knuth (1969). Seminumerical Algorithms.
;; The Art of Computer Programming, Volume 2. Addison Wesley

;#|
;(require ffi/unsafe)
;(define srand48 (get-ffi-obj 'srand48 #f (_fun _long -> _void)))
;(define drand48 (get-ffi-obj 'drand48 #f (_fun -> _double)))
;|#

(: poisson (Float -> Integer))
(define (poisson lmbd)
  (define L (exp (- lmbd)))
  (let: loop : Integer ([p : Float 1.] [k : Integer 0])    
    (define p* (* p #;(drand48) (random)))
    (if (p* . >= . L) (loop p* (+ k 1)) k)))

#;
(define (main s* p* z*)
  (define s (string->number s*))
  (define p (string->number p*))
  (define z (string->number z*))
  (run s p z))

(: run (Positive-Fixnum Fixnum Index -> (Listof Integer)))
(define (run s p z)
  (define ld (/ (exact->inexact z) p))
  ;(: k Fixnum)
  ;(: i Fixnum) (: n Fixnum)
  (: v (Vectorof Integer))
  (define v (make-vector z))
  (random-seed s) #;(srand48 s)
  (let outer ([n 0] [adj (* 2.0 (- ld 1))] [k 0])
    (: n* Integer) (: k* Integer)
    (define-values (n* k*)
      (let inner ([i 0] [n 0] [k k])
        (cond [(< i z)
               (vector-set! v n i)
               (define i* (+ 1 i)) 
               (define k* (poisson (* adj (sin (* (/ (+ i* 0.5) (+ z 1.)) 1.5707963268)))))
               (define i** (+ i* k*))
               (inner i** (+ n 1) k*)]
              [else (values n k)])))
    (unless (= n* p)
      (define adj* (cond [(> n* p) (* adj 1.02)]
                         [(< n* p) (/ adj 1.02)]
                         [else adj]))
      (outer n* adj* k*))) 
  (for/list: : (Listof Integer) ([e v] [_ (in-range p)]) e))

(define-predicate pos-fix? Positive-Fixnum)

(command-line
 #:args
 (#{s : String} #{p : String} #{z : String})
 (void (run (assert (string->number s) pos-fix?) (assert (string->number p) pos-fix?) (assert (string->number z) index?))))

#;
(unless (equal? '(0 1 2 4 5 6 7 10 11 13 14 16 18 20 22 23 27 32 34 37 41 45 52 58 65 68 77 81 90 96)
                (run 1 30 100))
  (error 'fail))

poisson-untyped.rkt

#lang racket/base

;; translated from an algorithm in C from
;; Hyberts SG, Takeuchi K, Wagner G (2010) Poisson-gap sampling and
;; forward maximum entropy reconstruction for enhancing the resolution
;; and sensitivity of protein NMR Data. J Am Chem Soc 132:2145–2147

;; generate random Poisson-distributed numbers as given
;; by Donald E. Knuth (1969). Seminumerical Algorithms.
;; The Art of Computer Programming, Volume 2. Addison Wesley


(require racket/cmdline)


(define-syntax-rule (while condition body ...)
  (let L () (when condition body ... (L))))
(define-syntax-rule (do-while condition body ...)
  (let L () (when (let () body ... condition) (L))))

;(: poisson : (Float -> Integer))
(define (poisson lmbd)
  (define L (exp (- lmbd)))
  (define k 0)
  (define p 1.)
  (do-while (p . >= . L)
    (define u (random))
    (set! p (* p u))
    (set! k (+ k 1)))
  (- k 1))

;(: run (Exact-Positive-Integer Fixnum Fixnum -> (Listof Integer)))
(define (run s p z)
  (define ld (/ (exact->inexact z) p))
  (define adj (* 2.0 (- ld 1)))
  (define i 0) (define j 0)
  (define k 0) (define n 0)
  (define v (make-vector z))
  (random-seed s)
  (do-while (not (= n p))
    (set! i 0)
    (set! n 0)
    (while (< i z)
      (vector-set! v n i)
      (set! i (+ i 1))
      (set! k (poisson (* adj (sin (* (/ (+ i 0.5) (+ z 1.)) 1.5707963268)))))
      (set! i (+ i k))
      (set! n (+ n 1)))
    (when (> n p) (set! adj (* adj 1.02)))
    (when (< n p) (set! adj (/ adj 1.02))))
  (for/list ([j (in-range p)]) (vector-ref v j)))


(command-line
 #:args
 (s p z)
 (void (run (string->number s) (string->number z) (string->number p))))

poisson.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
// generate random Poisson-distributed numbers as given
// by Donald E. Knuth (1969). Seminumerical Algorithms.
// The Art of Computer Programming, Volume 2. Addison Wesley
int poisson ( double lmbd )
{
  double L = exp( -lmbd );
  int k = 0;
  double p = 1;
  do {
    double u = drand48();
    p *= u;
    k += 1;
  } while ( p >= L );
  return( k-1 );
}
int main ( int argc, char** argv )
{
  float s = atof( argv[1] ); // seed value e.g. 1.0
  int p = atoi( argv[2] ); // sampled # of indices e.g 64
  int z = atoi( argv[3] ); // total # of indices e.g. 256
  int i; // Fourier grid index, e.g. 1 through 256
  int k; // generated gap size
  int n; // temporary # indices
  int *v; // temporary storage vector
  int j; // wrking variable
  float ld = (float) z / (float) p; // establish 1/fraction
  float adj = 2.0*(ld-1); // initial guess of adjustment
  srand48( s );
  v = ( int* ) malloc( z*sizeof( z ) );
  do {
    i = 0; n = 0;
    while ( i < z ) {
      v[n] = i;
      i += 1;
      k =poisson(adj*sin((float)(i+0.5)/(float)(z+1)*1.5707963268) );
      i += k;
      n += 1;
    }
    if ( n > p ) adj *= 1.02; // too many pts created
    if ( n < p ) adj /= 1.02; // too few pts created
  } while ( n != p ); // if not at first, try, try again
  //for ( j = 0 ; j < p ; j++ ) printf( "%d\n", v[j] );
}

poisson.rkt

#lang typed/racket

;; translated from an algorithm in C from
;; Hyberts SG, Takeuchi K, Wagner G (2010) Poisson-gap sampling and
;; forward maximum entropy reconstruction for enhancing the resolution
;; and sensitivity of protein NMR Data. J Am Chem Soc 132:2145–2147

;; generate random Poisson-distributed numbers as given
;; by Donald E. Knuth (1969). Seminumerical Algorithms.
;; The Art of Computer Programming, Volume 2. Addison Wesley

(define-syntax-rule (while condition body ...)
  (let L () (when condition body ... (L))))
(define-syntax-rule (do-while condition body ...)
  (let L () (when (let () body ... condition) (L))))

(: poisson : (Float -> Integer))
(define (poisson lmbd)
  (define L (exp (- lmbd)))
  (define k 0)
  (define p 1.)
  (do-while (p . >= . L)
    (define u (random) #;(drand48))
    (set! p (* p u))
    (set! k (+ k 1)))
  (- k 1))

(: run (Exact-Positive-Integer Fixnum Fixnum -> (Listof Integer)))
(define (run s p z)
  (define ld (/ (exact->inexact z) p))
  (define adj (* 2.0 (- ld 1)))
  (define i 0) (define j 0)
  (define k 0) (define n 0)
  (define v (make-vector z))
  (random-seed s) #;(srand48 s)
  (do-while (not (= n p))
    (set! i 0)
    (set! n 0)
    (while (< i z)
      (vector-set! v n i)
      (set! i (+ i 1))
      (set! k (poisson (* adj (sin (* (/ (+ i 0.5) (+ z 1.)) 1.5707963268)))))
      (set! i (+ i k))
      (set! n (+ n 1)))
    (when (> n p) (set! adj (* adj 1.02)))
    (when (< n p) (set! adj (/ adj 1.02))))
  (for/list ([j (in-range p)]) (vector-ref v j)))

(define-predicate pos-fix? Positive-Fixnum)

(command-line
 #:args
 (#{s : String} #{p : String} #{z : String})
 (void (run (assert (string->number s) pos-fix?) (assert (string->number p) pos-fix?) (assert (string->number z) index?))))