{{ message }}

Instantly share code, notes, and snippets.

# samth/poisson-ffi.rkt

Created Sep 12, 2011
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 #include #include // 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?))))