Changeset 28024 in project for release/4/nemo/trunk/nemo-matlab.scm


Ignore:
Timestamp:
01/03/13 15:10:11 (8 years ago)
Author:
Ivan Raikov
Message:

nemo: overhaul of ionic concentrations interface

File:
1 edited

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/nemo-matlab.scm

    r28020 r28024  
    33;; An extension for generating Matlab/Octave code from NEMO models.
    44;;
    5 ;; Copyright 2008-2012 Ivan Raikov and the Okinawa Institute of Science and Technology
     5;; Copyright 2008-2013 Ivan Raikov and the Okinawa Institute of Science and Technology
    66;;
    77;; This program is free software: you can redistribute it and/or
     
    108108                             (doc:space)))))))))
    109109
     110
    110111(define (format-fncall/MATLAB indent op args)
    111112  (let ((op1 (doc:text (->string op))))
    112113    (doc:cons op1 (group/MATLAB ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
     114
    113115
    114116(define (name-normalize expr)
     
    500502                 
    501503           (map (lambda (pool-ion)
    502                   (let ((n (third pool-ion))
    503                         (b (first pool-ion)))
     504                  (let ((n (pool-ion-in pool-ion))
     505                        (b (pool-ion-inq pool-ion)))
    504506                    (list n b)))
    505                 pool-ions)))
     507                pool-ions)
     508         
     509           (map (lambda (pool-ion)
     510                  (let ((n (pool-ion-out pool-ion))
     511                        (b (pool-ion-outq pool-ion)))
     512                    (list n b)))
     513                pool-ions)
     514           ))
    506515         
    507516         (eq-dag
     
    538547  (if v-eq
    539548      (let ((vi (lookup-def 'v state-index-map)))
    540         (if vi
    541             (pp indent+ ,(expr->string/MATLAB (second v-eq) (s+ "dy(" vi ")"))))))
     549        (if vi (pp indent+ ,(expr->string/MATLAB (second v-eq) (s+ "dy(" vi ")"))))))
    542550
    543551  (case method
     
    546554  (pp indent ,nl (end))
    547555)
     556
    548557
    549558(define (output-steadystate  sysname method globals steady-state-index-map pool-ions
     
    582591                 
    583592                 (map (lambda (pool-ion)
    584                         (let ((n (third pool-ion))
    585                               (b (first pool-ion)))
     593                        (let ((n (pool-ion-in pool-ion))
     594                              (b (pool-ion-inq pool-ion)))
    586595                          (list n b)))
    587596                      pool-ions)))
     
    630639(define (output-init sysname globals state-index-map steady-state-index-map
    631640                     external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    632                      reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
    633 
    634     (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
    635    
    636     (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
    637    
    638     (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
    639    
    640     (if (member 'v globals)
    641         (let ((vi (lookup-def 'v state-index-map)))
    642           (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
    643           (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
    644 
    645     (for-each (lambda (def)
    646                 (let ((n (first def))
    647                       (b (second def)))
    648                   (pp indent+ ,(expr->string/MATLAB b n))))
    649               external-eq-defs)
    650 
    651 
     641                     reaction-eq-defs i-eqs pool-ions perm-ions defaults indent indent+)
     642
     643
     644  (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
     645 
     646  (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
     647 
     648  (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
     649 
     650  (if (member 'v globals)
     651      (let ((vi (lookup-def 'v state-index-map)))
     652        (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
     653        (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
     654 
     655 
    652656    (let* ((init-eqs
    653657            (append
     658             
     659             (map (lambda (def)
     660                    (let ((n (first def))
     661                          (b (second def)))
     662                      (list n b)))
     663                  defaults)
     664             
     665             (map (lambda (def)
     666                    (let ((n (first def))
     667                          (b (second def)))
     668                      (list n b)))
     669                  external-eq-defs)
    654670             
    655671             const-defs
     
    658674             
    659675             (map (lambda (pool-ion)
    660                     (let ((n (third pool-ion))
    661                           (b (first pool-ion)))
     676                    (let ((n (pool-ion-in pool-ion))
     677                          (b (pool-ion-inq pool-ion)))
    662678                      (list n b)))
    663                   pool-ions)))
     679                  pool-ions)
     680
     681             (map (lambda (pool-ion)
     682                    (let ((n (pool-ion-out pool-ion))
     683                          (b (pool-ion-outq pool-ion)))
     684                      (list n b)))
     685                  pool-ions)
     686             ))
     687           
    664688
    665689           (init-dag
     
    667691                   (cons (first def) (enum-freevars (second def) '() '())))
    668692                 init-eqs))
     693
    669694           (init-order
    670695            (reverse
     
    748773             
    749774             (gate-complex-info    (nemo:gate-complex-query sys))
     775             (defaults             (nemo:defaults-query sys))
    750776
    751777             (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
     
    755781                                 (lookup-def 'acc-ions gate-complex-info)))
    756782             (epools        (lookup-def 'pool-ions gate-complex-info))
    757              (pool-ions     (map (lambda (lst) (map matlab-name lst)) epools))
     783             (pool-ions     (pool-ion-name-map matlab-name epools))
    758784
    759785             (i-gates       (lookup-def 'i-gates gate-complex-info))
     
    762788             (mcap          (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
    763789               
     790             (rcomp           (any (match-lambda ((name 'membrane-resistance id) (list name id)) (else #f)) components))
     791             (mres            (and rcomp (car ((dis 'component-exports) sys (cid rcomp)))))
     792               
     793             (rarea           (any (match-lambda ((name 'membrane-area id) (list name id)) (else #f)) components))
     794             (marea           (and rarea (car ((dis 'component-exports) sys (cid rarea)))))
     795               
     796
    764797
    765798             (i-eqs (filter-map
     
    913946                                                        (map third   acc-ions)
    914947                                                        (map fourth  acc-ions)
    915                                                         (map second  pool-ions)
    916                                                         (map third  pool-ions)
     948                                                        (map pool-ion-in  pool-ions)
     949                                                        (map pool-ion-out pool-ions)
    917950                                                        (map first   const-defs)
    918951                                                        (map first   asgn-eq-defs)
     
    921954
    922955                                                        ))))
    923              (v-eq    (if (and mcap (member 'v globals))
    924                           (list 'v (rhsexpr/MATLAB `(/ (neg ,(sum i-names)) ,mcap)))
    925                           (list 'v 0.0)))
     956
     957
     958             (v-eq    (cond ((and mcap mres marea (member 'v globals) (not (null? i-names)))
     959                             (list 'v (rhsexpr/MATLAB `(/ (* (neg ,(sum i-names)) ,marea) (* 2. ,mres ,mcap)))))
     960                            ((and marea (member 'v globals) (not (null? i-names)))
     961                             (list 'v (rhsexpr/MATLAB `(* (neg ,(sum i-names)) ,marea))))
     962                            ((and (member 'v globals) (not (null? i-names)))
     963                             (list 'v (rhsexpr/MATLAB `(neg ,(sum i-names)))))
     964                            (else #f)))
     965
    926966             
    927967             (state-index-map  (let ((acc (fold (lambda (def ax)
     
    9641004                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
    9651005         acc-ions)
    966 
    967            (for-each
    968             (lambda (p)
    969               (let ((pool-ion  (car p)))
    970                 (if (assoc pool-ion perm-ions)
    971                     (nemo:error 'nemo:matlab-translator
    972                                 ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
    973             pool-ions)
    974 
    975            (let ((output  (case mode
    976                             ((single)  (open-output-file (make-output-fname dirname prefix ".m" filename)))
    977                             (else #f))))
     1006       
     1007        (let ((output  (case mode
     1008                         ((single)  (open-output-file (make-output-fname dirname prefix ".m" filename)))
     1009                         (else #f))))
    9781010
    9791011             (if output (with-output-to-port output (lambda () (output-globals))))
     
    10051037                   (output-init sysname globals state-index-map steady-state-index-map
    10061038                                external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    1007                                 reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
     1039                                reaction-eq-defs i-eqs pool-ions perm-ions defaults indent indent+)
    10081040                   (pp indent ,nl)))
    10091041               (if (not output) (close-output-port output1)))
     
    10341066
    10351067
    1036 #|
    1037  
    1038 (define (make-define-dfn fenv)
    1039 
    1040   (define (D vars body)
    1041     (let loop ((vars vars) (exprs (list)))
    1042       (if (null? vars) (if (pair? (cdr exprs))
    1043                            `(+ . ,(reverse exprs))  (car exprs))
    1044           (let ((de (let ((de (differentiate fenv (car vars) body)))
    1045                        (let loop ((e (simplify de)) (de de))
    1046                          (if (equal? de e) de
    1047                              (loop (simplify e) e))))))
    1048             (loop (cdr vars) (cons de exprs))))))
    1049 
    1050   (lambda (indent n proc)
    1051     (let ((lst (procedure-data proc))
    1052           (indent+ (+ 2 indent)))
    1053       (let ((retval   (matlab-name (gensym "retval")))
    1054             (rt       (lookup-def 'rt lst))
    1055             (formals  (lookup-def 'formals lst))
    1056             (vars     (lookup-def 'vars lst))
    1057             (body     (lookup-def 'body lst)))
    1058         (pp indent ,nl (function ,retval = ,(s+ "d_" (matlab-name n) ) (,(slp ", " vars)) ))
    1059         (let* ((dbody (D vars body))
    1060                (dbody (canonicalize-expr/MATLAB (rhsexpr/MATLAB dbody))))
    1061           (pp indent+ ,(expr->string/MATLAB dbody retval))
    1062           (pp indent end))
    1063           )))
    10641068)
    1065 
    1066 (define (output-jac sysname method globals fenv state-index-map
    1067                     rate-eq-defs reaction-eq-defs asgn-eq-defs
    1068                     mcap i-eqs v-eq
    1069                     defuns indent indent+)
    1070 
    1071 
    1072   (define (D var expr)
    1073     (let* ((dexpr  (distribute expr))
    1074            (de (differentiate fenv var dexpr))
    1075            (se (let loop ((e (simplify de)) (de de))
    1076                  (if (equal? de e) de
    1077                      (loop (simplify e) e)))))
    1078       se))
    1079    
    1080   (pp indent ,nl (function jac = ,(s+ sysname "_jac")
    1081                            (,(slp ", " `(t y)))))
    1082 
    1083   (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
    1084  
    1085   (if (member 'v globals)
    1086       (let ((vi (lookup-def 'v state-index-map)))
    1087         (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
    1088  
    1089   (for-each (lambda (def)
    1090               (let ((name (first def)))
    1091                 (pp indent+ ,(expr->string/MATLAB
    1092                               (s+ "y(" (lookup-def name state-index-map) ")") name))))
    1093             rate-eq-defs)
    1094 
    1095   (for-each (lambda (def)
    1096               (let ((n (first def)) (b (second def)))
    1097                 (pp indent+
    1098                     ,(expr->string/MATLAB b n))))
    1099             reaction-eq-defs)
    1100  
    1101   (pp indent+ ,(expr->string/MATLAB `(zeros (length y) (length y)) 'jac))
    1102 
    1103 
    1104   (for-each
    1105    (lambda (def)
    1106      
    1107      (let* ((name  (first def))
    1108             (rhs   (second def))
    1109             (fv    (enum-freevars rhs '() '() ))
    1110             (drvs  (map (lambda (st)
    1111                           (let ((var (first st))
    1112                                 (expr `(let ,(or (and v-eq (list v-eq)) '())
    1113                                          (let ,rate-eq-defs
    1114                                            (let ,reaction-eq-defs
    1115                                              (let ,(filter-map (lambda (x)
    1116                                                                  (or (assoc x i-eqs)
    1117                                                                        (assoc x asgn-eq-defs)))
    1118                                                                fv))
    1119                                              ,rhs)))))
    1120                             (D var expr)))
    1121                         state-index-map))
    1122             (vs (map (lambda (st)
    1123                        (let ((ni (lookup-def name state-index-map))
    1124                              (i (second st)))
    1125                          (s+ "jac(" (slp "," (list ni i)) ")")))
    1126                      state-index-map))
    1127             (ndrvs  (map (lambda (x) (canonicalize-expr/MATLAB (rhsexpr/MATLAB x))) drvs)))
    1128        
    1129        (pp indent+ ,(s+ "% state " name))
    1130        (for-each (lambda (ndrv v)
    1131                    (pp indent+ ,(expr->string/MATLAB ndrv v)))
    1132                  ndrvs vs)
    1133        ))
    1134            
    1135      (if (and mcap (member 'v globals))
    1136          (let ((vi (lookup-def 'v state-index-map)))
    1137            (if vi (cons
    1138                    (list 'v `(/ (neg ,(sum (map first i-eqs))) ,mcap))
    1139                    rate-eq-defs)
    1140                rate-eq-defs))
    1141          rate-eq-defs))
    1142  
    1143   (pp indent ,nl (end))
    1144  
    1145 )
    1146 
    1147 |#
    1148 
    1149 )
Note: See TracChangeset for help on using the changeset viewer.