source: project/ode/tags/2.7a/extensions/ode-rate.scm @ 7359

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

Created releases with new author URL.

File size: 4.2 KB
Line 
1;;
2;;
3;; An extension for 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 ode)
23(require-extension digraph)
24(require-extension environments)
25
26(define (rate-name s0 s1)
27  (string->symbol (string-append (symbol->string s0) "->" (symbol->string s1))))
28
29(define (sum lst)
30  (if (null? lst) lst
31      (match lst
32             ((x)   x)
33             ((x y) `(+ ,x ,y))
34             ((x y . rest) `(+ (+ ,x ,y) ,(sum rest)))
35             ((x . rest) `(+ ,x ,(sum rest))))))
36
37
38(define (lookup-field k lst . rest)
39  (let-optionals rest ((default #f))
40   (let ((v (alist-ref k lst)))
41     (if v (first v) default))))
42
43 
44(define (ode:rate-transformer sys)
45  (let ((ode        (match (environment-ref sys (ode-intern 'odecore))
46                           (($ ode:quantity 'ODECORE value)  value)))
47        (new-env    (ode:env-copy sys #t)))
48    (let ((env-extend! ((ode 'env-extend!) new-env))
49          (eqdef! ((ode 'eqdef!) new-env)))
50      (environment-for-each sys 
51       (lambda (sym en)
52         (match en
53                (('rate . alst)
54                 (let* ((name       (lookup-field 'name alst))
55                        (g          (make-digraph name (string-append (symbol->string name) " rate graph")))
56                        (add-node!  (g 'add-node!))
57                        (add-edge!  (g 'add-edge!))
58                        (out-edges  (g 'out-edges))
59                        (in-edges   (g 'in-edges))
60                        (node-info  (g 'node-info))
61                        (node-list  (let loop ((lst (list)) (alst alst))
62                                      (if (null? alst)  (delete-duplicates lst eq?)
63                                          (match (car alst) 
64                                                 (('-> s0 s1 rate-expr)
65                                                  (loop (cons s0 (cons s1 lst)) (cdr alst)))
66                                                 (('-> _)
67                                                  (ode:error 'ode:rate-transformer ": invalid rate equation " 
68                                                             (car alst) " in system " en))
69                                                 (else (loop lst (cdr alst)))))))
70                        (node-ids      (list-tabulate (length node-list) identity))
71                        (name->id-map  (zip node-list node-ids)))
72                   ;; insert state nodes in the dependency graph
73                   (for-each (lambda (i n) (add-node! i n)) node-ids node-list)
74                   ;; create rate edges in the graph
75                   (for-each (lambda (e) 
76                               (match e (('-> s0 s1 rate-expr)
77                                         (let ((i  (car (alist-ref s0 name->id-map)))
78                                               (j  (car (alist-ref s1 name->id-map))))
79                                           (add-edge! (list i j (rate-name s0 s1)))))
80                                      (else (void))))
81                             alst)
82                   ;; generate assigned quantities to hold the computed rates
83                   (for-each
84                    (lambda (x)
85                      (match x 
86                             (('-> s0 s1 rate-expr)
87                              (match (list (environment-ref sys s0)
88                                           (environment-ref sys s1))
89                                     ((($ ode:quantity 'STATE n0 in0 v0 df0 rhs0 aerr0 rerr0 d0)
90                                       ($ ode:quantity 'STATE n1 in1 v1 df1 rhs1 aerr1 rerr1 d1))
91                                      (let ((r (rate-name s0 s1)))
92                                        (env-extend! r '(asgn) 'none  `(* ,s0 ,rate-expr))))
93                                     (else (ode:error 'ode:kinetic-transformer 
94                                                      ": rate equations require state quantities"))))
95                             (else (void))))
96                    alst)
97                   ;; generate differential equations for each state in the rate equation system
98                   (for-each (lambda (n) 
99                               (let ((out   (out-edges (first n)))
100                                     (in    (in-edges (first n)))
101                                     (name  (second n))) 
102                                 (let ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
103                                                     `(+ (neg ,(sum (map third out)))
104                                                         ,(sum (map third in))))
105                                                    ((and (not (null? out)) (null? in))
106                                                     `(neg ,(sum (map third out))))
107                                                    ((and (null? out) (not (null? in)))
108                                                     (sum (map third in))))))
109                                   (match (environment-ref sys name)
110                                          (($ ode:quantity 'STATE n in v df rhs aerr rerr d)
111                                           (eqdef! name rhs1))
112                                          (else (ode:error 'ode:kinetic-transformer 
113                                                           ": rate equations require state quantities"))))))
114                             ((g 'nodes)))
115                   ))
116
117                (else (void))))))
118      new-env))
Note: See TracBrowser for help on using the repository browser.