Changeset 13012 in project


Ignore:
Timestamp:
01/15/09 09:14:55 (11 years ago)
Author:
Ivan Raikov
Message:

Save.

Location:
release/3/nemo/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • release/3/nemo/trunk/nemo-macros.scm

    r12710 r13012  
    22;; NEMO macros
    33;;
    4 ;; Copyright 2008 Ivan Raikov and the Okinawa Institute of Science and Technology
     4;; Copyright 2008-2009 Ivan Raikov and the Okinawa Institute of Science and Technology
    55;;
    66;; This program is free software: you can redistribute it and/or
  • release/3/nemo/trunk/nemo-matlab.scm

    r12960 r13012  
    3232(require-extension nemo-core)
    3333(require-extension nemo-utils)
     34(require-extension nemo-ionch)
    3435
    3536(define-extension nemo-matlab)
     
    218219  (let-optionals rest ((rv #f) (width 72))
    219220    (sdoc->string (doc:format width (format-expr/MATLAB 2 x rv)))))
     221
     222
     223(define (expeuler dt name rhs)
     224  (define (isname? x) (equal? x name))
     225  (let ((res
     226         (match rhs
     227                ((or ('- A ('* B (and x (? isname?))))
     228                     ('+ ('neg ('* B (and x (? isname?)))) A))
     229                 (let ((xexp (string->symbol (s+ x 'exp))))
     230                   `(let ((,xexp (exp (* (neg ,B) ,dt))))
     231                      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B))))))
     232
     233                ((or ('- A ('* (and x (? isname?)) . B))
     234                     ('+ ('neg ('* (and x (? isname?)) . B)) A))
     235                 (let ((xexp (string->symbol (s+ x 'exp)))
     236                       (B1   (if (null? (cdr B)) (car B) `(* ,@B))))
     237                   `(let ((,xexp (exp (* (neg ,B1) ,dt))))
     238                      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B1))))))
     239               
     240                (('+ ('neg ('* (and x1 (? isname?)) Alpha))
     241                     ('* ('- 1 (and x2 (? isname?))) Beta))
     242                 (let ((A  Alpha)
     243                       (B  `(+ ,Alpha ,Beta)))
     244                   (let ((xexp (string->symbol (s+ x1 'exp))))
     245                     `(let ((,xexp (exp (* (neg ,B) ,dt))))
     246                        (+ (* ,x1 ,xexp) (* (- 1 ,xexp) (/ ,A ,B)))))))
     247               
     248                (('let bnds body)
     249                 `(let ,bnds ,(expeuler dt name body)))
     250               
     251                (else (nemo:error 'nemo:expeuler ": unable to rewrite equation " rhs
     252                                  "in exponential Euler form")))))
     253
     254    res))
    220255 
    221256(define (make-define-fn table? min-v max-v with)
     
    236271
    237272
    238 (define (state-eqs n initial open transitions power)
     273
     274(define (state-init n init)
     275  (let* ((init  (rhsexpr/MATLAB init))
     276         (init1 (canonicalize-expr/MATLAB init)))
     277    (list (matlab-name n) init1)))
     278
     279
     280(define (asgn-eq n rhs)
     281  (let* ((fbody   (rhsexpr/MATLAB rhs))
     282         (fbody1  (canonicalize-expr/MATLAB fbody)))
     283    (list (matlab-name n) fbody1)))
     284
     285
     286(define (reaction-eq n open transitions)
     287  (list (matlab-name n) (matlab-name (matlab-state-name n open))))
     288
     289
     290(define (reaction-transition-eqs n initial open transitions power method)
    239291  (match-let (((g  node-subs)  (transitions-graph n open transitions matlab-state-name)))
    240292     (let* ((out-edges  (g 'out-edges))
    241293            (in-edges   (g 'in-edges))
    242             (nodes      ((g 'nodes)))
    243             (snode      (find (lambda (s) (not (eq? (second s) open))) nodes)))
    244       ;; generate equations for each state in the transitions system
     294            (nodes   ((g 'nodes)))
     295            (snode   (find (lambda (s) (not (eq? (second s) open))) nodes)))
     296       ;; generate differential equations for each state in the transitions system
    245297      (let ((eqs    (fold (lambda (s ax)
    246298                            (if (= (first snode) (first s) ) ax
     
    256308                                                      ((and (null? out) (not (null? in)))
    257309                                                       (sum (map third in)))))
    258                                          (fbody0  (rhsexpr/MATLAB rhs1))
    259                                          (fbody1  (canonicalize-expr/MATLAB fbody0)))
    260                                     (cons (list name fbody1) ax)))))
     310                                         (fbody0 (rhsexpr/MATLAB rhs1)))
     311                                    (case method
     312                                      ((expeuler)  (cons (list name (canonicalize-expr/MATLAB
     313                                                                     (expeuler 'dt name fbody0)))
     314                                                         ax))
     315                                      (else        (cons (list name (canonicalize-expr/MATLAB fbody0))
     316                                                         ax))
     317                                      )))))
    261318                          (list) nodes)))
    262319        eqs))))
    263 
    264 
    265 (define (state-init n init)
    266   (let* ((init  (rhsexpr/MATLAB init))
    267          (init1 (canonicalize-expr/MATLAB init)))
    268     (list (matlab-name n) init1)))
    269 
    270 
    271 (define (asgn-eq n rhs)
    272   (let* ((fbody   (rhsexpr/MATLAB rhs))
    273          (fbody1  (canonicalize-expr/MATLAB fbody)))
    274     (list (matlab-name n) fbody1)))
    275 
    276 
    277 (define (reaction-eq n open transitions)
    278   (list (matlab-name n) (matlab-name (matlab-state-name n open))))
    279 
     320           
     321       
    280322
    281323(define (poset->asgn-eq-defs poset sys)
     
    288330                               (cases nemo:quantity en
    289331                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
     332                                      (else  ax))
     333                               ax))))
     334            ax lst))
     335   (list) poset))
     336
     337(define (poset->rate-eq-defs poset sys method)
     338  (fold-right
     339   (lambda (lst ax)
     340     (fold  (lambda (x ax)
     341              (match-let (((i . n)  x))
     342                         (let ((en (environment-ref sys n)))
     343                           (if (nemo:quantity? en)
     344                               (cases nemo:quantity en
     345
     346                                      (REACTION  (name initial open transitions conserve power)
     347                                                 (append (reaction-transition-eqs name initial open transitions
     348                                                                                  power method) ax))
     349                                     
     350                                      (RATE (name initial rhs)
     351                                            (let ((fbody0  (rhsexpr/MATLAB rhs))
     352                                                  (dy      name ))
     353                                              (case method
     354                                                ((expeuler) 
     355                                                 (cons (list dy (canonicalize-expr/MATLAB (expeuler 'dt name fbody0)))
     356                                                       ax))
     357                                                (else
     358                                                 (cons (list dy (canonicalize-expr/MATLAB fbody0)) ax)))))
     359
    290360                                      (else  ax))
    291361                               ax))))
     
    323393                                                            (state-init (matlab-state-name name open) name) ax)
    324394                                                     ax))
     395
     396                                      (RATE  (name initial rhs)
     397                                             (if (nemo:rhs? initial)
     398                                                 (cons (state-init name initial) ax)
     399                                                 ax))
     400                                     
    325401                                      (else  ax))
    326402                               ax))))
     
    365441                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
    366442                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
    367 
    368 
    369 (define (collect-epools sys)
    370    (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
    371       (let recur ((comp-name (nemo-intern 'toplevel)) (ax (list)))
    372         (let* ((comp-symbols   ((dis 'component-symbols)   sys comp-name))
    373                (subcomps       ((dis 'component-subcomps)  sys comp-name)))
    374           (fold recur
    375                 (fold (lambda (sym ax)
    376                         (let ((en (environment-ref sys sym)))
    377                           (match en
    378                                  ((or (('decaying 'pool)  ('name (? symbol? ion)) . alst)
    379                                       (('decaying-pool)   ('name (? symbol? ion)) . alst))
    380                                   (cons (list ion alst) ax))
    381                                  (else ax)))) ax comp-symbols)
    382                 (map third subcomps))))))
    383443
    384444
     
    401461             (defuns      ((dis 'defuns)  sys))
    402462             (components  ((dis 'components) sys))
    403              (ionchs      (filter-map (match-lambda ((name 'ion-channel id) (list name id)) (else #f)) components))
    404              (capcomp     (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
    405              (epools      (collect-epools sys)))
    406 
    407         (match-let (((state-list asgn-list g) deps*))
    408          (let* (
    409                 (const-defs       (filter-map
    410                                    (lambda (nv)
    411                                      (and (not (member (first nv) matlab-builtin-consts))
    412                                           (let ((v1 (canonicalize-expr/MATLAB (second nv))))
    413                                             (list (matlab-name (first nv)) v1))))
    414                                    consts))
    415                 (poset            (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
    416                 (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
    417                 (mcap             (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
    418                 (perm-ions        (fold (lambda (ionch ax)
    419                                           (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
    420                                                  (perm      (lookup-def 'permeating-substance subcomps)))
    421                                             (if perm
    422                                                 (case (cn perm)
    423                                                   ((non-specific)   
    424                                                    (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
    425                                                           (i    (matlab-name 'i))
    426                                                           (e    (matlab-name 'e)))
    427                                                      (cons `(,(cn perm) ,i ,e ,erev) ax)))
    428                                                   (else (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
    429                                                                (i    (matlab-name (s+ 'i (cn perm))))
    430                                                                (e    (matlab-name (s+ 'e (cn perm)))))
    431                                                           (cons `(,(cn perm) ,i ,e ,erev) ax))))
    432                                                 ax)))
    433                                         (list) ionchs))
    434                (acc-ions           (fold (lambda (ionch ax)
    435                                            (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
    436                                                   (acc   (lookup-def 'accumulating-substance subcomps))
    437                                                   (i     (and acc (matlab-name (s+ 'i (cn acc)))))
    438                                                   (in    (and acc (matlab-name (s+ (cn acc) 'i))))
    439                                                   (out   (and acc (matlab-name (s+ (cn acc) 'o)))))
    440                                              (if acc  (cons `(,(cn acc) ,i ,in ,out) ax) ax)))
    441                                          (list) ionchs))
    442                (pool-ions          (map (lambda (ep)
    443                                           (let ((ion (car ep)))
    444                                             `(,(matlab-name ion) ,(matlab-name (s+ 'i ion)) ,(matlab-name (s+ ion 'i)))))
    445                                         epools))
    446 
    447                (pool-ion-i        (map (lambda (ep) (let ((ion (car ep))) (matlab-name (s+ 'i ion))))
    448                                        epools))
    449 
    450                
    451                (i-eqs (filter-map
    452                           (lambda (ionch)
    453                            
    454                             (let* ((label     (first ionch))
    455                                    (n         (second ionch))
    456                                    (subcomps  ((dis 'component-subcomps) sys n))
    457                                    (acc       (lookup-def 'accumulating-substance subcomps))
    458                                    (perm      (lookup-def 'permeating-substance subcomps))
    459                                    (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
    460                                    (pore      (lookup-def 'pore subcomps))
    461                                    (gate      (lookup-def 'gate subcomps))
    462                                    (sts       (and gate ((dis 'component-exports) sys (cid gate)))))
    463                              
    464                               (cond ((and perm pore gate)
    465                                      (case (cn perm)
    466                                        ((non-specific)
    467                                         (let* ((i     (matlab-name 'i))
    468                                                (e     (car permqs))
    469                                                (gmax  (car ((dis 'component-exports) sys (cid pore))))
    470                                                (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
    471                                                (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
    472                                                (gion  `(* ,gmax ,@sptms)))
    473                                           (list i e gion)))
    474                                        (else
    475                                         (let* ((i     (matlab-name (s+ 'i (cn perm))))
    476                                                (e     (car permqs))
    477                                                (gmax  (car ((dis 'component-exports) sys (cid pore))))
    478                                                (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
    479                                                (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
    480                                                (gion  `(* ,gmax ,@sptms)))
    481                                           (list i e gion)))))
     463             
     464             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
     465             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
     466
     467             (const-defs       (filter-map
     468                                (lambda (nv)
     469                                  (and (not (member (first nv) matlab-builtin-consts))
     470                                       (let ((v1 (canonicalize-expr/MATLAB (second nv))))
     471                                         (list (matlab-name (first nv)) v1))))
     472                                consts))
     473             
     474             (ionch-info    (nemo:ionch-query sys))
     475             (ionchs        (lookup-def 'ion-channels ionch-info))
     476             (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(matlab-name i) ,(matlab-name e) ,erev)))
     477                                 (lookup-def 'perm-ions ionch-info)))
     478             (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map matlab-name (list i in out)))))
     479                                 (lookup-def 'acc-ions ionch-info)))
     480             (epools        (lookup-def 'pool-ions ionch-info))
     481             (pool-ions     (map (lambda (lst) (map matlab-name lst)) epools))
     482
     483             (i-gates       (lookup-def 'i-gates ionch-info))
     484
     485             (capcomp       (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
     486             (mcap          (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
     487               
     488
     489             (i-eqs0 (filter-map
     490                     (lambda (ionch)
     491                       
     492                       (let* ((label     (first ionch))
     493                              (n         (second ionch))
     494                              (subcomps  ((dis 'component-subcomps) sys n))
     495                              (acc       (lookup-def 'accumulating-substance subcomps))
     496                              (perm      (lookup-def 'permeating-substance subcomps))
     497                              (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
     498                              (pore      (lookup-def 'pore subcomps))
     499                              (gate      (lookup-def 'gate subcomps))
     500                              (sts       (and gate ((dis 'component-exports) sys (cid gate)))))
     501                         
     502                         (cond ((and perm pore gate)
     503                                (case (cn perm)
     504                                  ((non-specific)
     505                                   (let* ((i     (matlab-name 'i))
     506                                          (e     (car permqs))
     507                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
     508                                          (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
     509                                          (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
     510                                          (gion  `(* ,gmax ,@sptms)))
     511                                     (list i e gion)))
     512                                  (else
     513                                   (let* ((i     (matlab-name (s+ 'i (cn perm))))
     514                                          (e     (car permqs))
     515                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
     516                                          (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
     517                                          (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
     518                                          (gion  `(* ,gmax ,@sptms)))
     519                                     (list i e gion)))))
     520                               
     521                               ((and perm pore)
     522                                (case (cn perm)
     523                                  ((non-specific)
     524                                   (let* ((i     (matlab-name 'i))
     525                                          (e     (car permqs))
     526                                          (gmax  (car ((dis 'component-exports) sys (cid pore)))))
     527                                     (list i e gmax)))
     528                                  (else
     529                                   (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
     530                               
     531                               ((and acc pore gate)
     532                                (let* ((i     (matlab-name (s+ 'i (cn acc))))
     533                                       (gmax  (car ((dis 'component-exports) sys (cid pore))))
     534                                       (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
     535                                       (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
     536                                       (gion  `(* ,gmax ,@sptms)))
     537                                  (list i #f gion)))
     538                               (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
     539                               )))
     540                     ionchs))
     541               
     542             (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs0))
     543             
     544             (i-eqs  (fold (lambda (b ax)
     545                             (match b
     546                                    ((and ps ((i e gion) . rst)) 
     547                                     (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b))))
     548                                                                   ps))
     549                                                       (sum (map third ps))))
     550                                            (sum0  (rhsexpr/MATLAB sum))
     551                                            (sum1  (canonicalize-expr/MATLAB sum0)))
     552                                       (cons (list i sum1) ax)))
    482553                                   
    483                                     ((and perm pore)
    484                                      (case (cn perm)
    485                                        ((non-specific)
    486                                         (let* ((i     (matlab-name 'i))
    487                                                (e     (car permqs))
    488                                              (gmax  (car ((dis 'component-exports) sys (cid pore)))))
    489                                           (list i e gmax)))
    490                                        (else
    491                                         (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
     554                                    ((i e gion)
     555                                     (let* ((expr0  (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
     556                                            (expr1  (canonicalize-expr/MATLAB expr0)))
     557                                       (cons (list i expr1) ax)))
    492558                                   
    493                                     ((and acc pore gate)
    494                                      (let* ((i     (matlab-name (s+ 'i (cn acc))))
    495                                             (gmax  (car ((dis 'component-exports) sys (cid pore))))
    496                                             (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
    497                                             (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
    498                                             (gion  `(* ,gmax ,@sptms)))
    499                                        (list i #f gion)))
    500                                     (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
    501                                     )))
    502                           ionchs))
    503                
    504                   (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
    505                  
    506                   (i-eqs  (fold (lambda (b ax)
    507                                   (match b
    508                                          ((and ps ((i e gion) . rst)) 
    509                                           (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b))))
    510                                                                         ps))
    511                                                             (sum (map third ps))))
    512                                                  (sum0  (rhsexpr/MATLAB sum))
    513                                                  (sum1  (canonicalize-expr/MATLAB sum0)))
    514                                             (cons (list i sum1) ax)))
    515                                          
    516                                          ((i e gion)
    517                                           (let* ((expr0  (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
    518                                                  (expr1  (canonicalize-expr/MATLAB expr0)))
    519                                             (cons (list i expr1) ax)))
    520                                          
    521                                        (else ax)))
    522                                 (list) i-bkts))
    523 
    524                   (state-eq-defs      (reverse (poset->state-eq-defs poset sys)))
    525                  
    526                   (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
    527 
    528                   (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq))))
    529                                            (poset->state-conserve-eq-defs poset sys)))
    530 
    531                   (pool-consts       (concatenate
    532                                       (map (lambda (ep)
    533                                              (let* ((ep-name     (first ep))
    534                                                     (pool-ion    (assoc ep-name pool-ions))
    535                                                     (i-name      (second pool-ion))
    536                                                     (init-name   (matlab-name (s+ ep-name '-init)))
    537                                                     (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
    538                                                     (beta-name   (matlab-name (s+ ep-name '-beta)))
    539                                                     (depth-name  (matlab-name (s+ ep-name '-depth))))
    540                                                (list init-name temp-name beta-name depth-name)))
    541                                            epools)))
    542 
    543                   (pool-eq-defs
    544                    (map (lambda (ep)
    545                           (let* ((ep-name     (first ep))
    546                                  (pool-ion    (assoc ep-name pool-ions))
    547                                  (i-name      (second pool-ion))
    548                                  (init-name   (matlab-name (s+ ep-name '-init)))
    549                                  (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
    550                                  (beta-name   (matlab-name (s+ ep-name '-beta)))
    551                                  (depth-name  (matlab-name (s+ ep-name '-depth)))
    552                                  (rhs         `(let ((F 96485.0))
    553                                                  (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name))
    554                                                     (* ,beta-name ,ep-name .
    555                                                        ,(if temp-name (list temp-name) (list)))))))
    556                             `(,(s+ ep-name)  ,rhs)))
    557                         epools))
    558 
    559                   (rate-eq-defs     (append pool-eq-defs state-eq-defs))
    560                  
    561                   (init-eq-defs     (poset->state-init-defs poset sys))
    562 
    563                   (state-index-map  (let ((acc (fold (lambda (def ax)
    564                                                        (let ((st-name (first def)))
    565                                                          (list (+ 1 (first ax))
    566                                                                (cons `(,st-name ,(first ax)) (second ax)))))
    567                                                      (list 1 (list))
    568                                                      (cons (list 'v) rate-eq-defs))))
    569                                      
    570                                       (second acc)))
    571                  
    572                   (steady-state-index-map  (let ((acc (fold (lambda (def ax)
    573                                                               (let ((st-name (first def)))
    574                                                                 (if (not (alist-ref st-name init-eq-defs))
    575                                                                     (list (+ 1 (first ax))
    576                                                                           (cons `(,st-name ,(first ax)) (second ax)))
    577                                                                     ax)))
    578                                                             (list 1 (list))
    579                                                             rate-eq-defs)))
    580                                              (second acc)))
    581 
    582                   (globals          (map matlab-name
    583                                        (delete-duplicates (append
    584                                                            exports
    585                                                            pool-consts
    586                                                            (map second  perm-ions)
    587                                                            (map third   perm-ions)
    588                                                            (map second  acc-ions)
    589                                                            (map third   acc-ions)
    590                                                            (map fourth  acc-ions)
    591                                                            (map second  pool-ions)
    592                                                            (map third   pool-ions)
    593                                                            (map first   imports)
    594                                                            (map first   const-defs)))))
     559                                    (else ax)))
     560                           (list) i-bkts))
     561
     562             (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
     563             
     564             
     565             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys method)))
     566             
     567             (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
     568             
     569             (init-eq-defs       (poset->state-init-defs poset sys))
     570             
     571             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq))))
     572                                      (poset->state-conserve-eq-defs poset sys)))
     573             
     574             (pool-ion-i          (map second pool-ions))
     575             
     576             
     577             (state-index-map  (let ((acc (fold (lambda (def ax)
     578                                                  (let ((st-name (first def)))
     579                                                    (list (+ 1 (first ax))
     580                                                          (cons `(,st-name ,(first ax)) (second ax)))))
     581                                                (list 1 (list))
     582                                                (cons (list 'v) rate-eq-defs))))
     583                                 
     584                                 (second acc)))
     585             
     586             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
     587                                                         (let ((st-name (first def)))
     588                                                           (if (not (alist-ref st-name init-eq-defs))
     589                                                               (list (+ 1 (first ax))
     590                                                                     (cons `(,st-name ,(first ax)) (second ax)))
     591                                                               ax)))
     592                                                       (list 1 (list))
     593                                                       rate-eq-defs)))
     594                                        (second acc)))
     595             
     596             (globals          (map matlab-name
     597                                    (delete-duplicates (append
     598                                                        exports
     599                                                        (map second  perm-ions)
     600                                                        (map third   perm-ions)
     601                                                        (map second  acc-ions)
     602                                                        (map third   acc-ions)
     603                                                        (map fourth  acc-ions)
     604                                                        (map second  pool-ions)
     605                                                        (map third   pool-ions)
     606                                                        (map first   imports)
     607                                                        (map first   const-defs)))))
    595608
    596609               )
     
    612625            pool-ions)
    613626
    614            (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     627           (if (not (null? globals)) (pp indent (global ,(sl\ " " globals))))
    615628
    616629           (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
     
    701714               (pp indent ,nl (endfunction))))
    702715             
     716
     717         (pp indent ,nl (function  ,(s+ sysname "_print_state") (y)))
     718
     719         (let ((lst (sort (map (lambda (x) (cons (->string (car x)) (cdr x))) state-index-map)
     720                          (lambda (x y) (string<? (car x) (car y))))))
     721           (for-each (lambda (p)
     722                       (let ((n (first p)) (i (second p)))
     723                         (pp indent+ (,n " = " "y(" ,i ")"))))
     724                     lst))
     725
     726         (pp indent ,nl (endfunction))
     727
    703728         
    704729         (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
     
    708733         (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
    709734         
    710          (for-each (lambda (ep)
    711                      (let* ((ep-name     (first ep))
    712                             (ep-props    (second ep))
    713                             (init-expr   (lookup-def 'initial ep-props))
    714                             (temp-expr   (lookup-def 'temp-adj ep-props))
    715                             (beta-expr   (lookup-def 'beta ep-props))
    716                             (depth-expr  (lookup-def 'depth ep-props))
    717                             (init-name   (matlab-name (s+ ep-name '-init)))
    718                             (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
    719                             (beta-name   (matlab-name (s+ ep-name '-beta)))
    720                             (depth-name  (matlab-name (s+ ep-name '-depth))))
    721                        (if (or (not beta-expr) (not depth-expr) (not init-expr))
    722                            (nemo:error 'nemo:matlab-translator
    723                                        ": ion pool " ep-name " requires initial value, depth and beta parameters"))
    724                        (let ((temp-val  (and temp-expr (eval-const sys temp-expr)))
    725                              (init-val  (eval-const sys init-expr))
    726                              (beta-val  (eval-const sys beta-expr))
    727                              (depth-val (eval-const sys depth-expr)))
    728                          (pp indent+
    729                              ,(expr->string/MATLAB init-val init-name)
    730                              ,(expr->string/MATLAB temp-val temp-name)
    731                              ,(expr->string/MATLAB beta-val beta-name)
    732                              ,(expr->string/MATLAB depth-val depth-name)))))
    733                    epools)
    734          
    735 
    736735         (if (member 'v globals)
    737736             (let ((vi (lookup-def 'v state-index-map)))
     
    796795                           (apply define-fn (cons indent fndef))))
    797796                     defuns))
    798          )))))))
     797         )))))
  • release/3/nemo/trunk/nemo-nmodl.scm

    r13004 r13012  
    5757  (enum-freevars rhs (list) (list)))
    5858
    59 (define (rhsexpr expr)
     59(define (rhsexpr/NMODL expr)
    6060  (match expr
    61          (('if . es)  `(if . ,(map (lambda (x) (rhsexpr x)) es)))
     61         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/NMODL x)) es)))
    6262         (('pow x y)  (if (and (integer? y)  (positive? y))
    6363                          (if (> y 1)  (let ((tmp (gensym "x")))
     
    6565                              x)
    6666                            expr))
    67          ((s . es)    (if (symbol? s)  (cons s (map (lambda (x) (rhsexpr x)) es)) expr))
     67         ((s . es)    (if (symbol? s)  (cons s (map (lambda (x) (rhsexpr/NMODL x)) es)) expr))
    6868         (id          (if (symbol? id) (nmodl-name id) id))))
    6969
     
    322322            (body     (lookup-def 'body lst)))
    323323        (pp indent ,nl (FUNCTION ,n (,(sl\ ", " vars)) "{" ))
    324         (let* ((body1 (canonicalize-expr/NMODL (rhsexpr body)))
     324        (let* ((body1 (canonicalize-expr/NMODL (rhsexpr/NMODL body)))
    325325               (lbs   (enum-bnds body1 (list))))
    326326          (if (not (null? lbs)) (pp indent+ (LOCAL ,(sl\ ", " lbs))))
     
    368368 
    369369
    370 (define (reaction-eqs n initial open transitions power method)
     370(define (reaction-transition-eqs n initial open transitions power method)
    371371  (match-let (((g  node-subs)  (transitions-graph n open transitions nmodl-state-name)))
    372372     (let* ((out-edges  (g 'out-edges))
     
    388388                                                      ((and (null? out) (not (null? in)))
    389389                                                       (sum (map third in)))))
    390                                          (fbody0 (rhsexpr rhs1)))
     390                                         (fbody0 (rhsexpr/NMODL rhs1)))
    391391                                    (case method
    392392                                      ((expeuler)  (cons (list name (canonicalize-expr/NMODL (expeuler 'dt name fbody0)))
     
    452452
    453453(define (state-init n init)
    454   (let* ((init  (rhsexpr init))
     454  (let* ((init  (rhsexpr/NMODL init))
    455455         (init1 (canonicalize-expr/NMODL init)))
    456456    (list (nmodl-name n) init1)))
     
    458458
    459459(define (asgn-eq n rhs)
    460   (let* ((fbody   (rhsexpr rhs))
     460  (let* ((fbody   (rhsexpr/NMODL rhs))
    461461         (fbody1  (canonicalize-expr/NMODL fbody)))
    462462    (list (nmodl-name n) fbody1)))
     
    498498
    499499
    500 (define (poset->state-eq-defs poset sys kinetic method)
     500(define (poset->rate-eq-defs poset sys kinetic method)
    501501  (fold-right
    502502   (lambda (lst ax)
     
    507507                               (cases nemo:quantity en
    508508                                      (REACTION  (name initial open transitions conserve power)
    509                                                  (append (reaction-eqs name initial open transitions power method) ax))
     509                                                 (append (reaction-transition-eqs name initial open transitions
     510                                                                                  power method) ax))
    510511                                     
    511512                                      (RATE (name initial rhs)
    512                                             (let ((fbody0  (rhsexpr rhs))
     513                                            (let ((fbody0  (rhsexpr/NMODL rhs))
    513514                                                  (dy      name ))
    514515                                              (case method
     
    525526
    526527
    527 (define (poset->kstate-eq-defs poset sys kinetic)
     528(define (poset->kinetic-eq-defs poset sys kinetic)
    528529  (fold-right
    529530   (lambda (lst ax)
     
    631632             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
    632633             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
    633              (asgn-eq-defs  (poset->asgn-eq-defs poset sys))
    634634             (ionch-info    (nemo:ionch-query sys))
    635635             (ionchs        (lookup-def 'ion-channels ionch-info))
     
    645645             (has-kinetic?  (or (not (null? (filter (lambda (x) (member (car x) kinetic)) states)))))
    646646             (has-ode?      (or (not (null? (filter (lambda (x) (not (member (car x) kinetic))) states)))
    647                                 (not (null? pool-ions)))))
     647                                (not (null? pool-ions))))
     648
     649             (asgn-eq-defs       (poset->asgn-eq-defs poset sys))
     650             (reaction-eq-defs   (poset->reaction-eq-defs poset sys kinetic))
     651             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys kinetic method)))
     652             (kstate-eq-defs     (poset->kinetic-eq-defs poset sys kinetic))
     653             (conserve-eq-defs   (poset->state-conserve-eq-defs poset sys))
     654             (state-init-defs    (poset->state-init-defs poset sys))
     655
     656             )
    648657
    649658
     
    762771               (begin
    763772                 (pp indent ,nl (PROCEDURE reactions () "{"))
    764                  (let* ((eq-defs   (poset->reaction-eq-defs poset sys kinetic))
    765                         (locals    (find-locals (map second eq-defs))) )
     773                 (let ((locals    (find-locals (map second reaction-eq-defs))) )
    766774                   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
    767775                   (for-each (lambda (def)
    768776                               (let ((n (nmodl-name (first def))) (b (second def)))
    769                                  (pp indent+ ,(expr->string/NMODL b n)))) eq-defs))
     777                                 (pp indent+ ,(expr->string/NMODL b n)))) reaction-eq-defs))
    770778                 
    771779                 (pp indent "}")))
     
    852860                                            (if (null? ps)
    853861                                                (let* ((sum0  (sum summands))
    854                                                        (sum1  (rhsexpr sum0))
     862                                                       (sum1  (rhsexpr/NMODL sum0))
    855863                                                       (sum2  (canonicalize-expr/NMODL sum1)))
    856864                                                  (cons (list i sum2) ax))
     
    859867                                           
    860868                                         ((i e gion)
    861                                           (let* ((expr0  (rhsexpr (if e `(* ,gion (- v ,e)) gion)))
     869                                          (let* ((expr0  (rhsexpr/NMODL (if e `(* ,gion (- v ,e)) gion)))
    862870                                                 (expr1  (canonicalize-expr/NMODL expr0)))
    863871                                            (cons (list i expr1) ax)))
     
    879887           
    880888           (if has-ode?
    881                (let* ((eq-defs  (reverse (poset->state-eq-defs poset sys kinetic method)))
    882                       (locals   (find-locals (map second eq-defs))))
     889               (let ((locals   (find-locals (map second rate-eq-defs))))
    883890                 (case method
    884891                   ((expeuler) (pp indent ,nl (PROCEDURE states () "{")))
     
    892899                                     (b (second def)))
    893900                                 (pp indent+ ,(expr->string/NMODL b n))))
    894                              eq-defs))
     901                             rate-eq-defs))
    895902                 (pp indent "}")))
    896903
     
    898905               (begin
    899906                 (pp indent ,nl (KINETIC kstates "{"))
    900                  (let* ((keq-defs          (poset->kstate-eq-defs poset sys kinetic))
    901                         (locals            (concatenate (find-locals (map third (map second keq-defs)))))
    902                         (conserve-eq-defs  (poset->state-conserve-eq-defs poset sys)))
     907                 (let ((locals            (concatenate (find-locals (map third (map second kstate-eq-defs))))))
    903908                   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
    904909                   (for-each
     
    907912                             (eqs           (second def))
    908913                             (conserve-eqs  (lookup-def n conserve-eq-defs)))
    909                        
    910914                        (for-each
    911915                         (lambda (eq)
     
    926930                                    conserve-eqs))
    927931                        ))
    928                     keq-defs))
     932                    kstate-eq-defs))
    929933                 (pp indent "}")))
    930934           
    931935           
    932            (let* ((init-defs         (poset->state-init-defs poset sys))
    933                   (locals            (concatenate (find-locals (map second init-defs)))) )
     936           (let ((locals (concatenate (find-locals (map second state-init-defs)))) )
    934937               (pp indent ,nl (INITIAL "{"))
    935938               (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
     
    937940               (for-each (lambda (def)
    938941                           (let ((n (first def)) (b (second def)))
    939                              (pp indent+ ,(expr->string/NMODL b n)))) init-defs)
     942                             (pp indent+ ,(expr->string/NMODL b n)))) state-init-defs)
    940943
    941944               (if has-kinetic?
     
    944947               (pp indent "}")
    945948
     949               (pp indent ,nl (PROCEDURE print_state () "{"))
     950
     951               (let ((lst (sort (map (compose ->string first) rate-eq-defs) string<?)))
     952                 (for-each (lambda (x)
     953                             (pp indent+ (printf (,(s+ #\" x " = %g\\n"  #\") ", " ,x ))))
     954                           lst))
     955
     956               (pp indent "}")
     957
    946958               ))))))
  • release/3/nemo/trunk/nemo-utils.scm

    r12784 r13012  
    33;; Utility procedures for NEMO code generators.
    44;;
    5 ;; Copyright 2008 Ivan Raikov and the Okinawa Institute of Science and Technology
     5;; Copyright 2008-2009 Ivan Raikov and the Okinawa Institute of Science and Technology
    66;;
    77;; This program is free software: you can redistribute it and/or
Note: See TracChangeset for help on using the changeset viewer.