source: project/ode/trunk/extensions/ode-bpr.scm @ 7357

Last change on this file since 7357 was 7357, checked in by Ivan Raikov, 12 years ago

Miscelanneous egg updates.

File size: 8.3 KB
Line 
1;;
2;;
3;; An extension for binomial probabilistic rate equations in ODE systems.
4;;
5;; Copyright 2007 Ivan Raikov and the Okinawa Institute of Science and Technology
6;;
7;; This program is free software: you can redistribute it and/or
8;; modify it under the terms of the GNU General Public License as
9;; published by the Free Software Foundation, either version 3 of the
10;; License, or (at your option) any later version.
11;;
12;; This program is distributed in the hope that it will be useful, but
13;; WITHOUT ANY WARRANTY; without even the implied warranty of
14;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15;; General Public License for more details.
16;;
17;; A full copy of the GPL license can be found at
18;; <http://www.gnu.org/licenses/>.
19;;
20
21(require-extension srfi-1)
22(require-extension srfi-4)
23(require-extension ode)
24(require-extension digraph)
25(require-extension environments)
26(require-extension random-mtzig)
27
28(define (rate-name s0 s1 . rest)
29  (let-optionals rest ((prefix ""))
30    (string->symbol (string-append prefix (symbol->string s0) "->" (symbol->string s1)))))
31
32(define (vector-name prefix rate-name)
33  (string->symbol (string-append prefix "_" (symbol->string rate-name))))
34
35(define (sum lst)
36  (if (null? lst) lst
37      (match lst
38             ((x)   x)
39             ((x y) `(+ ,x ,y))
40             ((x y . rest) `(+ (+ ,x ,y) ,(sum rest)))
41             ((x . rest) `(+ ,x ,(sum rest))))))
42
43                               
44(define (lookup-field k lst . rest)
45  (let-optionals rest ((default #f))
46   (let ((v (alist-ref k lst)))
47     (if v (first v) default))))
48
49
50(define binom-seed (random-mtzig:init))
51
52; Embed but do not parse
53#>
54  #include <math.h>
55  double randmtzig_randb(int nnr, double ppr, void *state);
56
57  double cbinom_approx (double x, double n, double p, unsigned int *state)
58  {
59     int i;
60     double result;
61     double np = (n-1)*p;
62     double q = (1-p);
63     double xmin, xmax;
64
65     xmin = 0;
66     xmax = n-1;
67
68     if  (!((p < 1.0) && (np >= 10.0) && ((np *  q) >= 10.0)))
69     {
70      result = randmtzig_randb(n, p, state);
71     }
72     else
73     {
74       double mu;
75       double sigma;
76             
77       sigma = np*q;
78       mu = np;
79
80       result = mu + (sqrt(sigma) * x);
81     }
82     if (result < xmin) result = xmin;
83     if (result > xmax) result = xmax;
84
85     return result;
86  }
87
88  double cnormal_sample (double x, double mu, double sigma)
89  {
90     double R;
91
92     R = sqrt(sigma);
93   
94     return (mu + (R * x));
95  }
96
97  void cnormal_sample_vect (unsigned int n, double *x, double *y, double mu, double sigma)
98  {
99     double R; unsigned int i;
100
101     R = sqrt(sigma);
102   
103     for (i = 0; i < n; i++)
104     {
105        y[i] = (mu + (R * x[i]));
106     }
107  }
108<#
109
110(define cnormal_sample (foreign-lambda double "cnormal_sample" 
111                                       double double double))
112
113(define cnormal_sample_vect (foreign-lambda void "cnormal_sample_vect" 
114                                            integer f64vector f64vector double double))
115
116(define cbinom_approx (foreign-lambda double "cbinom_approx" 
117                                      double double double u32vector))
118
119(define (normal-sample x mu sigma)
120  (if (f64vector? x)
121      (let* ((n (f64vector-length x))
122             (y (make-f64vector n 0.0)))
123        (cnormal_sample_vect n x y mu sigma)
124        y)
125      (cnormal_sample x mu sigma)))
126
127 
128(define (ode:bpr-transformer sys . rest)
129  (let ((new-env    (ode:env-copy sys #t))
130        (ode        (match (environment-ref sys (ode-intern 'odecore))
131                           (($ ode:quantity 'ODECORE value)  value)))
132        (dt         (match (environment-ref sys (ode-intern 'timestep))
133                           (($ ode:quantity 'TIMESTEP name value)  name))))
134    (let ((env-extend! ((ode 'env-extend!) new-env))
135          (eqdef!  ((ode 'eqdef!) new-env)))
136      (if (not (environment-includes? sys 'std-normal-sample))
137          (let ((fn   (lambda (S st) (random-mtzig:f64vector-randn! S st))))
138            (env-extend! 'std-normal-sample '(prim) fn `((formals  (integer scheme-object))
139                                                         (rt  f64vector)))))
140      (if (not (environment-includes? sys 'binom-approx))
141          (let ((fn (lambda (x n p) (cbinom_approx x n p (car binom-seed)))))
142            (env-extend! 'binom-approx '(prim) fn `((formals (double double double))
143                                                    (rt double)))))
144      (if (not (environment-includes? sys 'f64vector-ref))
145          (env-extend! 'f64vector-ref '(prim) f64vector-ref `((formals (scheme-object integer))
146                                                              (rt double))))
147    (environment-for-each sys 
148       (lambda (sym en)
149         (let-values (((name seedv alst)
150                       (match en
151                              (('rate . alst)
152                               (values (lookup-field 'name alst)
153                                       (lookup-field 'seed alst)
154                                       alst)))))
155            (if name
156                (let* ((g  (make-digraph name (string-append (symbol->string name) " probability rate graph")))
157                       (seed       (gensym "seed"))
158                       (ssize      (string->symbol (string-append "sample-size_" (symbol->string name))))
159                       (samvec     (vector-name "sample-vector" name))
160                       (add-node!  (g 'add-node!))
161                       (add-edge!  (g 'add-edge!))
162                       (out-edges  (g 'out-edges))
163                       (in-edges   (g 'in-edges))
164                       (node-info  (g 'node-info))
165                       (node-list  (let loop ((lst (list)) (alst alst))
166                                     (if (null? alst)  (delete-duplicates lst eq?)
167                                         (match (car alst) 
168                                                (('-> s0 s1 rate-expr)
169                                                 (loop (cons s0 (cons s1 lst)) (cdr alst)))
170                                                (('-> _)
171                                                 (ode:error 'ode:pr-transformer ": invalid rate equation " 
172                                                            (car alst) " in system " en))
173                                                (else (loop lst (cdr alst)))
174                                                ))))
175                       (node-ids      (list-tabulate (length node-list) identity))
176                       (name->id-map  (zip node-list node-ids)))
177                  ;; insert state nodes in the dependency graph
178                  (for-each (lambda (i n) (add-node! i n)) node-ids node-list)
179                  ;; create rate edges in the graph
180                  (for-each (lambda (e) 
181                              (match e (('-> s0 s1 rate-expr)
182                                        (let ((i  (car (alist-ref s0 name->id-map)))
183                                              (j  (car (alist-ref s1 name->id-map))))
184                                          (add-edge! (list i j (rate-name s0 s1)))))
185                                     (else (void))))
186                            alst)
187                  ;; generate assigned quantities to hold the
188                  ;; computed rates
189                  (let ((qenv (make-environment)))
190                    (for-each (lambda (x) (environment-extend! qenv x (list))) node-list)
191                    (let ((npqs
192                           (cdr (fold
193                                 (lambda (x ax)
194                                   (match x 
195                                     (('-> s0 s1 rate-expr)
196                                      (match (list (environment-ref sys s0)
197                                                   (environment-ref sys s1))
198                                         ((($ ode:quantity 'STATE n0 in0 v0 df0 rhs0 aerr0 rerr0 d0)
199                                           ($ ode:quantity 'STATE n1 in1 v1 df1 rhs1 aerr1 rerr1 d1))
200                                          (let ((i    (first ax))
201                                                (q    (rate-name s0 s1))
202                                                (r    (rate-name s0 s1  "rate_"))
203                                                (pf   (rate-name s0 s1  "pf_"))
204                                                (p    (rate-name s0 s1  "p_"))
205                                                (oq   (environment-ref qenv s0)))
206                                            (((ode 'defun!) new-env) pf '(rate state dt) '(/ (* rate dt) state))
207                                            (((ode 'env-extend!) new-env) r '(asgn)  'none `(* ,s0 ,rate-expr))
208                                            (((ode 'env-extend!) new-env) p '(asgn)  'none `(,pf ,r (- ,s0 . ,oq) ,dt))
209                                            (((ode 'env-extend!) new-env) q '(asgn)  'none 
210                                             `(binom-approx (f64vector-ref ,samvec (fix ,i)) (- ,s0 . ,oq) ,p))
211                                            (environment-set! qenv s0 (cons q oq))
212                                            (let ((ns  (second ax))
213                                                  (ps  (third  ax))
214                                                  (qs  (fourth ax)))
215                                              (list (+ i 1) (cons s0 ns) (cons p ps) (cons q qs)))))
216                                         (else (ode:error 'ode:pr-transformer 
217                                                          ": rate equations require state quantities"))))
218                                     (else ax)
219                                     ))
220                                 (list 0 (list) (list) (list)) alst))))
221                      (let ((S (length ((g 'edges)))))
222                        (((ode 'env-extend!) new-env) ssize '(prim) (length ((g 'edges))))
223                        ;; generate a declaration for an RNG seed
224                        ;; associated with this rate system
225                        (((ode 'env-extend!) new-env) seed '(prim) seedv)
226                        ;; generate a declaration for a sample vector
227                        ;; associated with this rate system
228                        (((ode 'env-extend!) new-env) samvec `(vasgn ,S) 'none `(std-normal-sample ,ssize ,seed))
229                        )))
230                  ;; generate difference equations for each state in the rate equation system
231                  (for-each (lambda (n) 
232                              (let ((out   (out-edges (first n)))
233                                    (in    (in-edges (first n)))
234                                    (name  (second n))) 
235                                (let ((rhs1  `(+ ,name
236                                                 ,(cond ((and (not (null? out)) (not (null? in)))
237                                                         `(+ (neg ,(sum (map third out)))
238                                                             ,(sum (map third in))))
239                                                        ((and (not (null? out)) (null? in))
240                                                         `(neg ,(sum (map third out))))
241                                                        ((and (null? out) (not (null? in)))
242                                                         (sum (map third in)))))))
243                                  (eqdef! name rhs1 'q))))
244                            ((g 'nodes)))
245                  (environment-remove! new-env sym)
246                  )))))
247       new-env)))
248
Note: See TracBrowser for help on using the repository browser.