Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Poisson-distributed random numbers in Typed Racket
#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))))
#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))
#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))))
#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] );
}
#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?))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment