Changeset 12685 in project for release/3/nemo/trunk/nemo-matlab.scm


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

Bug fixes.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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         )))))))
Note: See TracChangeset for help on using the changeset viewer.