Changeset 12685 in project


Ignore:
Timestamp:
12/01/08 06:47:14 (13 years ago)
Author:
Ivan Raikov
Message:

Bug fixes.

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

Legend:

Unmodified
Added
Removed
  • release/3/nemo/trunk/examples/AKP06/Na.nemo

    r12367 r12685  
    77             
    88     (const temp_adj = (pow (3 ((celsius - 22) / 10))))
     9
     10     (component (type membrane-capacitance)
     11               (const C = 1)
     12               (output C))
    913     
    1014     (component (type ion-channel) (name Na)
     
    96100                      (conserve (1 = (C1 + C2 + C3 + C4 + C5 + O + B + I1 + I2 + I3 + I4 + I5 + I6)))
    97101                     
    98                      (initial-equilibrium
    99                       (0 = ((I1 * bi1) + (C2 * b01) + (neg (C1 * (fi1 + f01)) )))
    100                       (0 = ((C1 * f01) + (I2 * bi2) + (C3 * b02) + (neg (C2 * (b01 + fi2 + f02)) )))
    101                       (0 = ((C2 * f02) + (I3 * bi3) + (C4 * b03) + (neg (C3 * (b02 + fi3 + f03)) )))
    102                       (0 = ((C3 * f03) + (I4 * bi4) + (C5 * b04) + (neg (C4 * (b03 + fi4 + f04)) )))
    103                       (0 = ((C4 * f04) + (I5 * bi5) + ( O * b0O) + (neg (C5 * (b04 + fi5 + f0O)) )))
    104                       (0 = ((C5 * f0O) + ( B * bip) + (I6 * bin) + (neg ( O * (b0O + fip + fin)) )))
    105                       (0 = (( O * fip) + ( B * bip)))
    106                       (0 = ((C1 * fi1) + (I2 * b11) + (neg (I1 * (bi1 + f11)) )))
    107                       (0 = ((I1 * f11) + (C2 * fi2) + (I3 * b12) + (neg (I2 * (b11 + bi2 + f12)) )))
    108                       (0 = ((I2 * f12) + (C3 * fi3) + (I4 * bi3) + (neg (I3 * (b12 + bi3 + f13)) )))
    109                       (0 = ((I3 * f13) + (C4 * fi4) + (I5 * b14) + (neg (I4 * (b13 + bi4 + f14)) )))
    110                       (0 = ((I4 * f14) + (C5 * fi5) + (I6 * b1n) + (neg (I5 * (b14 + bi5 + f1n)) )))
    111                       (1 = (C1 + C2 + C3 + C4 + C5 + O + B + I1 + I2 + I3 + I4 + I5 + I6 )))
    112                      
    113102                     (open O)   (power 1)))
    114103                   
  • release/3/nemo/trunk/examples/AKP06/Narsg.nemo

    r12367 r12685  
    9696                      (conserve (1 = (C1 + C2 + C3 + C4 + C5 + O + B + I1 + I2 + I3 + I4 + I5 + I6)))
    9797                     
    98                      (initial-equilibrium       
    99                       (0 = ((I1 * bi1) + (C2 * b01) + (neg (C1 * (fi1 + f01)) )))
    100                       (0 = ((C1 * f01) + (I2 * bi2) + (C3 * b02) + (neg (C2 * (b01 + fi2 + f02)) )))
    101                       (0 = ((C2 * f02) + (I3 * bi3) + (C4 * b03) + (neg (C3 * (b02 + fi3 + f03)) )))
    102                       (0 = ((C3 * f03) + (I4 * bi4) + (C5 * b04) + (neg (C4 * (b03 + fi4 + f04)) )))
    103                       (0 = ((C4 * f04) + (I5 * bi5) + ( O * b0O) + (neg (C5 * (b04 + fi5 + f0O)) )))
    104                       (0 = ((C5 * f0O) + ( B * bip) + (I6 * bin) + (neg ( O * (b0O + fip + fin)) )))
    105                       (0 = (( O * fip) + ( B * bip)))
    106                       (0 = ((C1 * fi1) + (I2 * b11) + (neg (I1 * (bi1 + f11)) )))
    107                       (0 = ((I1 * f11) + (C2 * fi2) + (I3 * b12) + (neg (I2 * (b11 + bi2 + f12)) )))
    108                       (0 = ((I2 * f12) + (C3 * fi3) + (I4 * bi3) + (neg (I3 * (b12 + bi3 + f13)) )))
    109                       (0 = ((I3 * f13) + (C4 * fi4) + (I5 * b14) + (neg (I4 * (b13 + bi4 + f14)) )))
    110                       (0 = ((I4 * f14) + (C5 * fi5) + (I6 * b1n) + (neg (I5 * (b14 + bi5 + f1n)) )))
    111                       (1 = (C1 + C2 + C3 + C4 + C5 + O + B + I1 + I2 + I3 + I4 + I5 + I6 )))
    11298                     
    11399                     (open O)   (power 1)))
  • release/3/nemo/trunk/examples/hh-substrate.nemo

    r12567 r12685  
    3535            ;;
    3636            (state-complex (m (transitions (<-> C O (amf (v)) (bmf (v))) )
    37                               (initial (amf (Vrest) / (amf (Vrest) + bmf (Vrest))) )
     37                              ;(initial (amf (Vrest) / (amf (Vrest) + bmf (Vrest))) )
     38                              (conserve (1 = (C + O)))
    3839                              (open O)  (power 3)))
    3940
     
    4243            ;; the value of state complex h is determined by state h1 (open state)
    4344            (state-complex (h (transitions (<-> C O (ahf (v)) (bhf (v))))
    44                               (initial (ahf (Vrest) / (ahf (Vrest) + bhf (Vrest))) )
     45                              ;(initial (ahf (Vrest) / (ahf (Vrest) + bhf (Vrest))) )
     46                              (conserve (1 = (C + O)))
    4547                              (open O) (power 1)))
    4648           
     
    7678            ;; the value of state complex n is determined by state n1 (open state)
    7779            (state-complex (n (transitions (<-> C O (anf (v)) (bnf (v))) )
    78                               (initial (anf (Vrest) / (anf (Vrest) + bnf (Vrest))) )
     80                              ;(initial (anf (Vrest) / (anf (Vrest) + bnf (Vrest))) )
     81                              (conserve (1 = (C + O)))
    7982                              (open O)  (power 4)))
    8083           
  • release/3/nemo/trunk/nemo-core.scm

    r12558 r12685  
    112112  (CONST      (name symbol?) (value number?))
    113113  (TSCOMP     (name symbol?)
    114               (initial      (lambda (x) (or (rhs? x) (and (list? x) (every lineq? x)))))
     114              (initial      (lambda (x) (or (rhs? x) (not x))))
    115115              (open         (lambda (x) (or (symbol? x) (and (list? x) (every symbol? x) ))))
    116116              (transitions  (lambda (x) (and (list? x) (every transition? x))))
     
    919919                             (cond ((and (symbol? id) (list? alst))
    920920                                    (let* ((initial      (lookup-def 'initial alst))
    921                                            (initial-eq   (alist-ref 'initial-equilibrium alst))
    922921                                           (conserve-eq  (alist-ref 'conserve alst))
    923922                                           (power        (lookup-def 'power alst))
     
    940939                                                (or (alist-ref 'transitions alst) (list)))))
    941940
    942                                       (if (not (or initial initial-eq))
     941                                      (if (not (or initial conserve-eq))
    943942                                          (nemo:error 'eval-nemo-system-decls
    944943                                                      "state complex declarations require initial value or "
    945                                                       "initial equilibrium equations"))
    946 
    947                                       (let ((initial-eq
    948                                              (and initial-eq (map (lambda (eq)
    949                                                                     (if (expr? (third eq))
    950                                                                         `(,(first eq) = ,(parse-expr (third eq)))
    951                                                                         (nemo:error 'eval-nemo-system-decls
    952                                                                                     "invalid equation " eq)))
    953                                                                   initial-eq)))
    954                                             (conserve-eq
     944                                                      "conservation equation"))
     945
     946                                      (let ((conserve-eq
    955947                                             (and conserve-eq (map (lambda (eq)
    956948                                                                     (if (expr? (third eq))
     
    960952                                                                   conserve-eq))))
    961953
    962                                         (if (and initial-eq
    963                                                  (or (not (list? initial-eq)) (not (every lineq? initial-eq))))
    964                                             (nemo:error 'eval-nemo-system-decls
    965                                                         "initial equilibrium field in state complex declarations "
    966                                                         "must be a list of linear equations"))
    967                                        
    968954                                        (if (and (list? conserve-eq) (not (every lineq? conserve-eq)))
    969955                                            (nemo:error 'env-extend!
     
    973959                                        (let ((initialv (and initial (eval-const (parse-expr initial)))))
    974960                                          (apply env-extend!
    975                                                  (cons* id '(tscomp) (or initialv initial-eq) `(power ,power)
     961                                                 (cons* id '(tscomp) initialv `(power ,power)
    976962                                                        (alist-update! 'conserve conserve-eq
    977963                                                          (alist-update! 'transitions transitions alst))
  • release/3/nemo/trunk/nemo-matlab.scm

    r12624 r12685  
    4848                           (else #\_))))
    4949            (loop (cons c1 lst) (cdr cs)))))))
    50                
    51 (define (rhsvars rhs)
    52   (enum-freevars rhs (list) (list)))
    5350           
    5451(define (rhsexpr/MATLAB expr)
     
    239236
    240237
    241 (define (define-state indent n)
    242   (pp indent (,n)))
    243 
    244 
    245238(define (state-eqs n initial open transitions power)
    246239  (match-let (((g  node-subs)  (transitions-graph n open transitions matlab-state-name)))
     
    249242            (nodes      ((g 'nodes)))
    250243            (snode      (find (lambda (s) (not (eq? (second s) open))) nodes)))
    251       ;; generate differential equations for each state in the transitions system
     244      ;; generate equations for each state in the transitions system
    252245      (let ((eqs    (fold (lambda (s ax)
    253246                            (if (= (first snode) (first s) ) ax
     
    265258                                         (fbody0  (rhsexpr/MATLAB rhs1))
    266259                                         (fbody1  (canonicalize-expr/MATLAB fbody0)))
    267                                     (cons (list name  fbody1) ax)))))
     260                                    (cons (list name fbody1) ax)))))
    268261                          (list) nodes)))
    269262        eqs))))
     
    274267         (init1 (canonicalize-expr/MATLAB init)))
    275268    (list (matlab-name n) init1)))
    276 
    277269
    278270
     
    355347
    356348
    357 (define (find-locals defs)
    358   (concatenate (map (lambda (def) (match def (('let bnds _) (map first bnds)) (else (list)))) defs)))
     349(define (poset->state-conserve-eq-defs poset sys)
     350  (fold-right
     351   (lambda (lst ax)
     352     (fold  (lambda (x ax)
     353              (match-let (((i . n)  x))
     354                         (let ((en (environment-ref sys n)))
     355                           (if (nemo:quantity? en)
     356                               (cases nemo:quantity en
     357                                      (TSCOMP (name initial open transitions conserve power)
     358                                              (if (and (list? conserve) (every nemo:lineq? conserve))
     359                                                  (cons (state-lineqs (matlab-name name) transitions conserve
     360                                                                      matlab-state-name) ax)
     361                                                  ax))
     362                                      (else  ax))
     363                               ax))))
     364            ax lst))
     365   (list) poset))
    359366
    360367
     
    417424
    418425        (match-let (((state-list asgn-list g) deps*))
    419          (let* ((const-defs       (filter-map
     426         (let* (
     427                (const-defs       (filter-map
    420428                                   (lambda (nv)
    421429                                     (and (not (member (first nv) matlab-builtin-consts))
     
    454462                                          (let ((ion (car ep)))
    455463                                            `(,(matlab-name ion) ,(matlab-name (s+ 'i ion)) ,(matlab-name (s+ ion 'i)))))
    456                                         epools)))
    457 
    458            (for-each
    459             (lambda (a)
    460               (let ((acc-ion   (car a)))
    461                 (if (assoc acc-ion perm-ions)
    462                     (nemo:error 'nemo:matlab-translator
    463                                 ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
    464             acc-ions)
    465 
    466            (for-each
    467             (lambda (p)
    468               (let ((pool-ion  (car p)))
    469                 (if (assoc pool-ion perm-ions)
    470                     (nemo:error 'nemo:matlab-translator
    471                                 ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
    472             pool-ions)
    473 
    474 
    475            (for-each (lambda (ep)
    476                        (let* ((ep-name     (first ep))
    477                               (ep-props    (second ep))
    478                               (init-expr   (lookup-def 'initial ep-props))
    479                               (temp-expr   (lookup-def 'temp-adj ep-props))
    480                               (beta-expr   (lookup-def 'beta ep-props))
    481                               (depth-expr  (lookup-def 'depth ep-props))
    482                               (init-name   (matlab-name (s+ ep-name '-init)))
    483                               (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
    484                               (beta-name   (matlab-name (s+ ep-name '-beta)))
    485                               (depth-name  (matlab-name (s+ ep-name '-depth))))
    486                          (if (or (not beta-expr) (not depth-expr) (not init-expr))
    487                              (nemo:error 'nemo:matlab-translator
    488                                          ": ion pool " ep-name " requires initial value, depth and beta parameters"))
    489                          (let ((temp-val  (and temp-expr (eval-const sys temp-expr)))
    490                                (init-val  (eval-const sys init-expr))
    491                                (beta-val  (eval-const sys beta-expr))
    492                                (depth-val (eval-const sys depth-expr)))
    493                            (pp indent
    494                                ,(expr->string/MATLAB init-val init-name)
    495                                ,(expr->string/MATLAB temp-val temp-name)
    496                                ,(expr->string/MATLAB beta-val beta-name)
    497                                ,(expr->string/MATLAB depth-val depth-name)))))
    498                      epools)
    499          
    500            (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
    501 
    502            (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
    503                                                                 ((lsode)  `(y t))
    504                                                                 ((odepkg) `(t y))
    505                                                                 (else     `(y t)))))))
    506 
    507            (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
    508            
    509            (if (not (null? asgns))
    510                (for-each (lambda (def)
    511                            (let ((n (matlab-name (first def)) )
    512                                  (b (second def)))
    513                              (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs))
    514          
    515            (let* ((i-eqs (filter-map
     464                                        epools))
     465               
     466               (i-eqs (filter-map
    516467                          (lambda (ionch)
    517468                           
     
    589540                 
    590541                  (stcomp-eq-defs   (poset->stcomp-eq-defs poset sys))
    591                  
     542
     543                  (conserve-eq-defs  (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq))))
     544                                          (poset->state-conserve-eq-defs poset sys)))
    592545                  (pool-eq-defs
    593546                   (map (lambda (ep)
     
    605558                            `(,(s+ ep-name)  ,rhs)))
    606559                        epools))
     560
     561                  (rate-eq-defs     (append pool-eq-defs state-eq-defs))
    607562                 
    608                   (rate-eq-defs (append pool-eq-defs state-eq-defs))
     563                  (init-eq-defs     (poset->state-init-defs poset sys))
    609564                 
    610565                  (state-index-map  (let ((acc (fold (lambda (def ax)
     
    614569                                                     (list 1 (list))
    615570                                                     (cons (list 'v) rate-eq-defs))))
     571                                     
     572                                      (second acc)))
     573                 
     574                  (steady-state-index-map  (let ((acc (fold (lambda (def ax)
     575                                                              (let ((st-name (first def)))
     576                                                                (if (not (alist-ref st-name init-eq-defs))
     577                                                                    (list (+ 1 (first ax))
     578                                                                          (cons `(,st-name ,(first ax)) (second ax)))
     579                                                                    ax)))
     580                                                            (list 1 (list))
     581                                                            state-eq-defs)))
    616582                                                             
    617                                       (second acc)))
    618                   )
    619 
    620 
    621              (if (member 'v globals)
    622                  (let ((vi (lookup-def 'v state-index-map)))
    623                    (if vi (pp indent+ (v = ,(s+ "y(" vi ")") #\;)))))
    624 
    625              (for-each (lambda (def)
    626                          (let ((n (first def)) )
    627                            (pp indent+ (,n = ,(s+ "y(" (lookup-def n state-index-map) ");"))))) rate-eq-defs)
    628 
    629              (for-each (lambda (def)
    630                          (let ((n (first def)) (b (second def)))
    631                            (pp indent+ ,(expr->string/MATLAB b n)))) stcomp-eq-defs)
    632 
    633              (pp indent+ (dy = zeros(length(y) #\, 1) #\;))
    634 
    635              (for-each (lambda (def)
    636                          (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) (b (second def)))
    637                            (pp indent+ ,(expr->string/MATLAB b n)))) rate-eq-defs)
    638 
    639              (for-each (lambda (pool-ion)
    640                          (pp indent+ (,(third pool-ion) = ,(first pool-ion)))) pool-ions)
    641 
    642              (for-each (lambda (def)
    643                          (pp indent+ ,(expr->string/MATLAB (second def) (first def)))) i-eqs)
    644 
    645              (if (and mcap (member 'v globals))
    646                  (let ((vi (lookup-def 'v state-index-map)))
    647                    (if vi (pp indent+  ,(expr->string/MATLAB `(/ (neg ,(sum (map first i-eqs))) ,mcap) 
    648                                                              (s+ "dy(" vi ")"))))))
    649              
    650              (pp indent ,nl (endfunction))
    651 
    652              (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
    653 
    654              (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
    655 
    656              (let* ((init-defs         (poset->state-init-defs poset sys)))
    657 
    658                (pp indent+ (y0 = zeros(,(length state-index-map) #\, 1) #\;))
    659 
    660                (if (member 'v globals)
    661                    (let ((vi (lookup-def 'v state-index-map)))
    662                      (pp indent+ (v = Vinit #\;))
    663                      (if vi (pp indent+ (,(s+ "y0(" vi ")") = v #\;)))))
    664 
    665                
    666                (for-each (lambda (def)
    667                            (let ((n (first def)) (b (second def)))
    668                              (pp indent+ ,(expr->string/MATLAB b n)))) const-defs)
    669 
     583                                             (second acc)))
     584               )
     585
     586           (for-each
     587            (lambda (a)
     588              (let ((acc-ion   (car a)))
     589                (if (assoc acc-ion perm-ions)
     590                    (nemo:error 'nemo:matlab-translator
     591                                ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
     592            acc-ions)
     593
     594           (for-each
     595            (lambda (p)
     596              (let ((pool-ion  (car p)))
     597                (if (assoc pool-ion perm-ions)
     598                    (nemo:error 'nemo:matlab-translator
     599                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
     600            pool-ions)
     601
     602
     603           (for-each (lambda (ep)
     604                       (let* ((ep-name     (first ep))
     605                              (ep-props    (second ep))
     606                              (init-expr   (lookup-def 'initial ep-props))
     607                              (temp-expr   (lookup-def 'temp-adj ep-props))
     608                              (beta-expr   (lookup-def 'beta ep-props))
     609                              (depth-expr  (lookup-def 'depth ep-props))
     610                              (init-name   (matlab-name (s+ ep-name '-init)))
     611                              (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
     612                              (beta-name   (matlab-name (s+ ep-name '-beta)))
     613                              (depth-name  (matlab-name (s+ ep-name '-depth))))
     614                         (if (or (not beta-expr) (not depth-expr) (not init-expr))
     615                             (nemo:error 'nemo:matlab-translator
     616                                         ": ion pool " ep-name " requires initial value, depth and beta parameters"))
     617                         (let ((temp-val  (and temp-expr (eval-const sys temp-expr)))
     618                               (init-val  (eval-const sys init-expr))
     619                               (beta-val  (eval-const sys beta-expr))
     620                               (depth-val (eval-const sys depth-expr)))
     621                           (pp indent
     622                               ,(expr->string/MATLAB init-val init-name)
     623                               ,(expr->string/MATLAB temp-val temp-name)
     624                               ,(expr->string/MATLAB beta-val beta-name)
     625                               ,(expr->string/MATLAB depth-val depth-name)))))
     626                     epools)
     627         
     628           (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     629
     630           (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
     631                                                                ((lsode)  `(y t))
     632                                                                ((odepkg) `(t y))
     633                                                                (else     `(y t)))))))
     634
     635           (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     636           
     637           (if (member 'v globals)
     638               (let ((vi (lookup-def 'v state-index-map)))
     639                 (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
     640           
     641           (for-each (lambda (def)
     642                       (let ((n (first def)) )
     643                         (pp indent+ ,(expr->string/MATLAB
     644                                       (s+ "y(" (lookup-def n state-index-map) ")") n))))
     645                     rate-eq-defs)
     646           
     647           (for-each (lambda (def)
     648                       (let ((n (matlab-name (first def)) )
     649                             (b (second def)))
     650                         (pp indent+ ,(expr->string/MATLAB b n))))
     651                     asgn-eq-defs)
     652         
     653         (for-each (lambda (def)
     654                     (let ((n (first def)) (b (second def)))
     655                       (pp indent+ ,(expr->string/MATLAB b n))))
     656                   stcomp-eq-defs)
     657         
     658         (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'dy))
     659         
     660         (for-each (lambda (def)
     661                     (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) (b (second def)))
     662                       (pp indent+ ,(s+ "# state " (first def))
     663                           ,(expr->string/MATLAB b n))))
     664                   rate-eq-defs)
     665         
     666         (for-each (lambda (pool-ion)
     667                     (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
     668                   pool-ions)
     669         
     670         (for-each (lambda (def)
     671                     (pp indent+ ,(expr->string/MATLAB (second def) (first def))))
     672                   i-eqs)
     673         
     674         (if (and mcap (member 'v globals))
     675             (let ((vi (lookup-def 'v state-index-map)))
     676               (if vi (pp indent+  ,(expr->string/MATLAB `(/ (neg ,(sum (map first i-eqs))) ,mcap) 
     677                                                         (s+ "dy(" vi ")"))))))
     678         
     679         (pp indent ,nl (endfunction))
     680
     681         (if (not (null? steady-state-index-map))
     682             (begin
     683               (pp indent ,nl (function y = ,(s+ sysname "_steadystate") (x)))
     684               
     685               (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     686               
     687               (for-each
     688                (lambda (def)
     689                  (let* ((n   (first def))
     690                         (ni  (lookup-def n steady-state-index-map)))
     691                    (pp indent+ ,(expr->string/MATLAB (if ni (s+ "x(" ni ")") n) n))
     692                    ))
     693                state-eq-defs)
     694               
     695               (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'y))
     696               
    670697               (for-each (lambda (def)
    671698                           (let ((n (matlab-name (first def)) )
    672699                                 (b (second def)))
    673700                             (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)
    674 
    675                (for-each (lambda (def)
    676                            (let* ((n  (first def))
    677                                   (ni (lookup-def n state-index-map))
    678                                   (b (second def)))
    679                              (pp indent+ ,(expr->string/MATLAB b (if ni (s+ "y0(" ni ")") n)))))
    680                          init-defs)
    681 
    682 
    683                )
    684              (pp indent ,nl (endfunction))
     701               
     702               (for-each
     703                (lambda (def)
     704                  (let* ((n   (first def))
     705                         (ni  (lookup-def n steady-state-index-map))
     706                         (b   (second def)))
     707                    (pp indent+ ,(s+ "# state " n)
     708                        ,(expr->string/MATLAB b (if ni (s+ "y(" ni ")") n)))
     709                    ))
     710                state-eq-defs)
     711               
     712               (pp indent ,nl (endfunction))))
    685713             
    686              (pp indent ,nl)
    687 
    688              (let* ((with      (inexact->exact (round (/ (abs (- max-v min-v)) step))))
    689                     (define-fn (make-define-fn table? min-v max-v with)))
    690                (for-each (lambda (fndef)
    691                            (if (not (member (car fndef) builtin-fns))
    692                                (apply define-fn (cons indent fndef))))
    693                          defuns))
    694              
    695              
    696              ))))))))
    697 
     714         
     715         (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
     716         
     717         (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     718         
     719         (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
     720         
     721         (if (member 'v globals)
     722             (let ((vi (lookup-def 'v state-index-map)))
     723               (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
     724               (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
     725         
     726         (for-each (lambda (def)
     727                     (let ((n (first def)) (b (second def)))
     728                       (pp indent+ ,(expr->string/MATLAB b n))))
     729                   const-defs)
     730         
     731         (if (not (null? init-eq-defs))
     732             (for-each (lambda (def)
     733                         (let ((n (matlab-name (first def)) )
     734                               (b (second def)))
     735                           (pp indent+ ,(expr->string/MATLAB b n))))
     736                       asgn-eq-defs))
     737         
     738         (for-each (lambda (def)
     739                     (let* ((n  (first def))
     740                            (ni (lookup-def n state-index-map))
     741                            (b (second def)))
     742                       (pp indent+ ,(expr->string/MATLAB b (if ni (s+ "y0(" ni ")") n)))))
     743                   init-eq-defs)
     744         
     745         (if (not (null? steady-state-index-map))
     746             (begin
     747               (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'xs)
     748                   ,(expr->string/MATLAB
     749                     `(fsolve ,(s+ #\" sysname "_steadystate" #\" ) xs)
     750                     "[ys,fval,info]"))
     751               
     752               (for-each
     753                (lambda (def)
     754                  (let* ((n (first def))
     755                         (si (lookup-def n steady-state-index-map))
     756                         (ni (lookup-def n state-index-map)))
     757                    (if si (pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")"))))))
     758                state-eq-defs)
     759               ))
     760         
     761         (pp indent ,nl (endfunction))
     762         
     763         (pp indent ,nl)
     764         
     765         (let* ((with      (inexact->exact (round (/ (abs (- max-v min-v)) step))))
     766                (define-fn (make-define-fn table? min-v max-v with)))
     767           (for-each (lambda (fndef)
     768                       (if (not (member (car fndef) builtin-fns))
     769                           (apply define-fn (cons indent fndef))))
     770                     defuns))
     771         )))))))
  • release/3/nemo/trunk/nemo-nmodl.scm

    r12624 r12685  
    305305               
    306306
    307          
     307#|       
    308308(define (lineq->string/NMODL x val . rest)
    309309  (let-optionals rest ((width 72))
    310310    (s+ "~ " (sdoc->string (doc:format width (format-lineq/NMODL 2 x #f)))
    311311        " = " (number->string val))))
    312  
     312|# 
    313313         
    314314(define (conserve-lineq->string/NMODL x val . rest)
     
    464464
    465465
    466 (define (state-lineqs n transitions lineqs)
    467   (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
    468          (state-list     (let loop ((lst (list)) (tlst transitions))
    469                            (if (null? tlst)  (delete-duplicates lst eq?)
    470                                (match (car tlst)
    471                                       (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
    472                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    473                                       (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
    474                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    475                                       (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
    476                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    477                                       (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
    478                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    479                                       (else
    480                                        (nemo:error 'nemo:state-lineq ": invalid transition equation "
    481                                                    (car tlst) " in state complex " n))
    482                                       (else (loop lst (cdr tlst)))))))
    483          (state-subs     (fold (lambda (s ax) (subst-extend s (nmodl-state-name n s) ax)) subst-empty state-list))
    484          (lineqs1        (map (lambda (lineq) (match lineq ((i '= . expr) `(,i = . ,(subst-convert expr state-subs)))))
    485                               lineqs)))
    486     (list (nmodl-name n) lineqs1)))
    487 
    488466(define (asgn-eq n rhs)
    489467  (let* ((fbody   (rhsexpr rhs))
     
    576554   (list) poset))
    577555
    578 
    579 (define (poset->state-init-eq-defs poset sys)
    580   (fold-right
    581    (lambda (lst ax)
    582      (fold  (lambda (x ax)
    583               (match-let (((i . n)  x))
    584                          (let ((en (environment-ref sys n)))
    585                            (if (nemo:quantity? en)
    586                                (cases nemo:quantity en
    587                                       (TSCOMP (name initial open transitions conserve power)
    588                                               (if (and (list? initial) (every nemo:lineq? initial))
    589                                                   (cons (state-lineqs name transitions initial) ax)
    590                                                   ax))
    591                                       (else  ax))
    592                                ax))))
    593             ax lst))
    594    (list) poset))
    595 
    596 
    597556(define (poset->state-conserve-eq-defs poset sys)
    598557  (fold-right
     
    605564                                      (TSCOMP (name initial open transitions conserve power)
    606565                                              (if (and (list? conserve) (every nemo:lineq? conserve))
    607                                                   (cons (state-lineqs name transitions conserve) ax)
     566                                                  (cons (state-lineqs (nmodl-name name) transitions conserve
     567                                                                      nmodl-state-name) ax)
    608568                                                  ax))
    609569                                      (else  ax))
     
    858818                   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals)))))
    859819#|
     820                 This seems to cause a segmentation fault in nrnoc:
     821
    860822                 (if (and table? min-v max-v table-with)
    861823                     (pp indent+ (TABLE ,(sl\ ", " (map first asgn-eq-defs))
     
    10451007                        (if conserve-eqs
    10461008                            (for-each (lambda (eq)
    1047                                       (let ((val  (first eq))
    1048                                             (expr (third eq)))
    1049                                         (pp indent+ ,(conserve-lineq->string/NMODL expr val))))
     1009                                        (let ((val  (first eq))
     1010                                              (expr (third eq)))
     1011                                          (pp indent+ ,(conserve-lineq->string/NMODL expr val))))
    10501012                                    conserve-eqs))
    10511013                        ))
     
    10551017           
    10561018           (let* ((init-defs         (poset->state-init-defs poset sys))
    1057                   (init-eq-defs      (poset->state-init-eq-defs poset sys))
    10581019                  (locals            (concatenate (find-locals (map second init-defs)))) )
    10591020               (pp indent ,nl (INITIAL "{"))
     
    10631024                           (let ((n (first def)) (b (second def)))
    10641025                             (pp indent+ ,(expr->string/NMODL b n)))) init-defs)
    1065                (cond ((and linear? (not (null? init-eq-defs)) )
    1066                     (pp indent+ (SOLVE initial_equilibrium)))
    1067                    (has-kinetic?
    1068                     (pp indent+ (SOLVE kstates STEADYSTATE sparse))))
     1026
     1027               (if has-kinetic?
     1028                   (pp indent+ (SOLVE kstates STEADYSTATE sparse)))
    10691029               
    10701030               (pp indent "}")
    10711031
    1072                (if (and linear? (not (null? init-eq-defs)) )
    1073                    (begin
    1074                      (pp indent ,nl (LINEAR initial_equilibrium "{"))
    1075                      (for-each
    1076                       (lambda (x)
    1077                         (let ((lineqs  (second x)))
    1078                           (for-each (lambda (eq)
    1079                                       (let ((val  (first eq))
    1080                                             (expr (third eq)))
    1081                                         (pp indent+ ,(lineq->string/NMODL expr val))))
    1082                                     lineqs)))
    1083                       init-eq-defs)
    1084                      (pp indent "}")))
    10851032               ))))
    10861033        ))))
  • release/3/nemo/trunk/nemo-utils.scm

    r12624 r12685  
    3838          if-convert let-enum let-elim let-lift
    3939          s+ sw+ sl\ nl spaces ppf
    40           transitions-graph
     40          transitions-graph state-lineqs
    4141          ))
    4242
     
    8282         (('if c t e)
    8383          `(if ,(k c subst) ,(k t subst) ,(k e subst)))
     84
    8485         (('let bs e)
    85           (k `(let ,(map (lambda (b) `(,(car b) ,(k (cadr b) subst))) bs) ,(k e subst)) subst))
     86          (k `(let ,(map (lambda (b) `(,(car b) ,(k (cadr b) subst))) bs) ,(k e subst))
     87             subst))
     88
    8689         ((f . es)
    8790          (cons (k f subst) (map (lambda (e) (k e subst)) es)))
     91
    8892         ((? symbol? )  (lookup-def t subst t))
     93
    8994         ((? atom? ) t)))
    9095
     
    188193    (let* ((nodes  ((g 'nodes)))
    189194           (snode   (find (lambda (s) (not (eq? (second s) open))) nodes))
    190            (snex   `(- 1 ,(sum (filter-map (lambda (s) (and (not (= (first s) (first snode))) (second s))) nodes))))
     195           (snex   (let ((nodes/s (filter-map (lambda (s) (and (not (= (first s) (first snode))) (second s))) nodes))
     196                         (sumvar  (gensym "sum")))
     197                     `(let ((,sumvar ,(sum nodes/s))) (- 1 ,sumvar))))
    191198           (add-tredge (lambda (s0 s1 rexpr1 rexpr2)
    192                          (let ((i   (car (alist-ref s0 name->id-map)))
    193                                (j   (car (alist-ref s1 name->id-map)))
    194                                (x0  (if (eq? s0 (second snode)) snex s0))
    195                                (x1  (if (eq? s1 (second snode)) snex s1)))
    196                            (add-edge! (list i j `(* ,(subst-convert x0 node-subs)
    197                                                     ,(subst-convert rexpr1 node-subs))))
    198                            (if rexpr2 (add-edge! (list j i `(* ,(subst-convert x1 node-subs)
    199                                                                ,(subst-convert rexpr2 node-subs)))))))))
     199                         (let* ((i   (car (alist-ref s0 name->id-map)))
     200                                (j   (car (alist-ref s1 name->id-map)))
     201                                (x0  (if (eq? s0 (second snode)) snex s0))
     202                                (x1  (if (eq? s1 (second snode)) snex s1))
     203                                (ij-expr  `(* ,(subst-convert x0 node-subs) ,(subst-convert rexpr1 node-subs)))
     204                                (ji-expr  (and rexpr2
     205                                               `(* ,(subst-convert x1 node-subs) ,(subst-convert rexpr2 node-subs)))))
     206                           (add-edge! (list i j ij-expr))
     207                           (if rexpr2 (add-edge! (list j i ji-expr)))))))
    200208      ;; create rate edges in the graph
    201209      (for-each (lambda (e)
     
    210218      (list g node-subs))))
    211219
     220
     221(define (state-lineqs n transitions lineqs state-name)
     222  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
     223         (state-list     (let loop ((lst (list)) (tlst transitions))
     224                           (if (null? tlst)  (delete-duplicates lst eq?)
     225                               (match (car tlst)
     226                                      (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
     227                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     228                                      (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
     229                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     230                                      (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
     231                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     232                                      (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
     233                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     234                                      (else
     235                                       (nemo:error 'nemo:state-lineq ": invalid transition equation "
     236                                                   (car tlst) " in state complex " n))
     237                                      (else (loop lst (cdr tlst)))))))
     238         (state-subs     (fold (lambda (s ax) (subst-extend s (state-name n s) ax)) subst-empty state-list))
     239         (lineqs1        (map (lambda (lineq) (match lineq ((i '= . expr) `(,i = . ,(subst-convert expr state-subs)))))
     240                              lineqs)))
     241    (list n lineqs1)))
Note: See TracChangeset for help on using the changeset viewer.