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


Ignore:
Timestamp:
11/18/08 08:31:57 (11 years ago)
Author:
Ivan Raikov
Message:

Many fixes to the matlab code generator.

File:
1 edited

Legend:

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

    r12353 r12554  
    9999(define (s+ . lst)    (string-concatenate (map ->string lst)))
    100100(define (sw+ lst)     (string-intersperse (filter-map (lambda (x) (and x (->string x))) lst) " "))
    101 (define (s\ p . lst)  (string-intersperse (map ->string lst) p))
    102101(define (sl\ p lst)   (string-intersperse (map ->string lst) p))
    103102(define nl "\n")
     
    337336    (sdoc->string (doc:format width (format-expr/MATLAB 2 x rv)))))
    338337 
    339 
    340 (define (format-lineq/MATLAB indent expr . rest) 
    341   (let-optionals rest ((rv #f))
    342    (let ((indent+ (+ 2 indent)))
    343     (match expr
    344        (('let bindings body)
    345         (letblk/MATLAB
    346          (fold-right
    347            (lambda (x ax)
    348              (letblk/MATLAB
    349               (match (second x)
    350                      (('if c t e)
    351                       (ifthen/MATLAB
    352                        (group/MATLAB (format-lineq/MATLAB indent c))
    353                        (block/MATLAB (format-lineq/MATLAB indent t (first x)))
    354                        (block/MATLAB (format-lineq/MATLAB indent e (first x)))))
    355                      (else
    356                       (format-op/MATLAB indent+ " = "
    357                                        (list (format-lineq/MATLAB indent (first x) )
    358                                              (format-lineq/MATLAB indent (second x))))))
    359               ax))
    360            (doc:empty) bindings)
    361          (let ((body1 (doc:nest indent (format-lineq/MATLAB indent body))))
    362            (if rv  (format-op/MATLAB indent " = " (list (format-lineq/MATLAB indent+ rv ) body1))
    363                body1))))
    364        
    365        (('if . rest) (error 'format-lineq/MATLAB "invalid if statement " expr))
    366 
    367        ((op . rest) 
    368         (let ((op (case op ((pow)  '^) ((abs) 'fabs) (else op))))
    369           (let ((fe
    370                 (if (member op matlab-ops)
    371                     (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
    372                           (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
    373                           (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
    374                       (case op
    375                         ((/) 
    376                          (format-op/MATLAB indent op
    377                                           (map (lambda (x)
    378                                                  (let ((fx (format-lineq/MATLAB indent+ x)))
    379                                                    (if (or (symbol? x) (number? x)) fx
    380                                                        (if (or mul? plmin?) (group/MATLAB fx) fx)))) rest)))
    381                         ((*) 
    382                          (format-op/MATLAB indent op
    383                                           (map (lambda (x)
    384                                                  (let ((fx (format-lineq/MATLAB indent+ x)))
    385                                                    (if (or (symbol? x) (number? x)) fx
    386                                                        (if plmin? (group/MATLAB fx) fx)))) rest)))
    387                        
    388                         ((^) 
    389                          (format-op/MATLAB indent op
    390                                           (map (lambda (x)
    391                                                  (let ((fx (format-lineq/MATLAB indent+ x)))
    392                                                    (if (or (symbol? x)  (number? x)) fx
    393                                                        (if (or mdiv? plmin?) (group/MATLAB fx) fx)))) rest)))
    394                        
    395                         (else
    396                          (format-op/MATLAB indent op
    397                                           (map (lambda (x)
    398                                                  (let ((fx (format-lineq/MATLAB indent+ x))) fx)) rest)))))
    399                    
    400                     (case op
    401                       ((neg) (format-op/MATLAB indent '* (map (lambda (x) (format-lineq/MATLAB indent+ x))
    402                                                              (cons "(-1)" rest))))
    403                       (else  (format-fncall/MATLAB indent op (map (lambda (x) (format-lineq/MATLAB indent+ x))
    404                                                                  rest)))))))
    405 
    406            (if rv (format-op/MATLAB indent " = " (list (format-lineq/MATLAB indent+ rv ) fe)) fe))))
    407      
    408       (else  (let ((fe (doc:text (->string expr))))
    409                (if rv
    410                    (format-op/MATLAB indent " = " (list (format-lineq/MATLAB indent+ rv ) fe))
    411                    fe)))))))
    412                
    413 
    414          
    415 (define (lineq->string/MATLAB x val . rest)
    416   (let-optionals rest ((width 72))
    417     (s+ "~ " (sdoc->string (doc:format width (format-lineq/MATLAB 2 x #f)))
    418         " = " (number->string val))))
    419  
    420 
    421 (define (make-define-fn table? min-v max-v with depend)
     338(define (make-define-fn table? min-v max-v with)
    422339  (lambda (indent n proc)
    423340    (let ((lst (procedure-data proc))
     
    431348        (let* ((body1 (canonicalize-expr/MATLAB (rhsexpr body)))
    432349               (lbs   (enum-bnds body1 (list))))
    433 ;         (if (and table? min-v max-v with)
    434 ;             (match vars
    435 ;                    (('v)  (pp indent+ (TABLE ,@(if depend `(DEPEND ,depend) `(""))
    436 ;                                              FROM ,min-v TO ,max-v WITH ,with)))
    437 ;                    (else  (void))))
    438350          (pp indent+ ,(expr->string/MATLAB body1 retval))
    439351          (pp indent endfunction))
     
    510422                                         (fbody  (rhsexpr rhs1))
    511423                                         (fbody1 (canonicalize-expr/MATLAB fbody)))
    512                                     (cons (list (s+ name "'") fbody1) ax)))))
     424                                    (cons (list name fbody1) ax)))))
    513425                          (list) nodes)))
    514426        eqs))))
     
    516428       
    517429
    518 (define (kstate-eqs n initial open transitions power)
    519   (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
    520          (state-list     (let loop ((lst (list)) (tlst transitions))
    521                            (if (null? tlst)  (delete-duplicates lst eq?)
    522                                (match (car tlst)
    523                                       (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
    524                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    525                                       (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
    526                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    527                                       (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
    528                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    529                                       (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
    530                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    531                                       (else
    532                                        (nemo:error 'nemo:matlab-kstate-eqs ": invalid transition equation "
    533                                                    (car tlst) " in state complex " n))
    534                                       (else (loop lst (cdr tlst)))))))
    535          (state-subs     (fold (lambda (s ax) (subst-extend s (matlab-state-name n s) ax)) subst-empty state-list)))
    536     ;; generate kinetic equations for each edge in the transitions system
    537     (map
    538      (lambda (e)
    539        (match e
    540               (('-> s0 s1 rexpr)
    541                (let ((i  (lookup-def s0 state-subs))
    542                      (j  (lookup-def s1 state-subs)))
    543                  `(-> ,i ,j ,(subst-convert rexpr state-subs))))
    544              
    545               ((s0 '-> s1 rexpr)
    546                (let ((i  (lookup-def s0 state-subs))
    547                      (j  (lookup-def s1 state-subs)))
    548                  `(-> ,i ,j ,(subst-convert rexpr state-subs))))
    549              
    550               (('<-> s0 s1 rexpr1 rexpr2)
    551                (let ((i  (lookup-def s0 state-subs))
    552                      (j  (lookup-def s1 state-subs)))
    553                  `(<-> ,i ,j ,(subst-convert rexpr1 state-subs) ,(subst-convert rexpr2 state-subs))))
    554              
    555               ((s0 '<-> s1 rexpr1 rexpr2)
    556                (let ((i  (lookup-def s0 state-subs))
    557                      (j  (lookup-def s1 state-subs)))
    558                  `(<-> ,i ,j ,(subst-convert rexpr1 state-subs) ,(subst-convert rexpr2 state-subs))))
    559              
    560                  
    561               (else (nemo:error 'nemo:matlab-kstate-eqs ": invalid transition equation "
    562                                 e " in state complex " n))))
    563      transitions)))
    564        
    565 
    566430
    567431(define (state-init n init)
     
    570434    (list (matlab-name n) init1)))
    571435
    572 
    573 (define (state-init-eq n transitions init)
    574   (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
    575          (state-list     (let loop ((lst (list)) (tlst transitions))
    576                            (if (null? tlst)  (delete-duplicates lst eq?)
    577                                (match (car tlst)
    578                                       (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
    579                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    580                                       (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
    581                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    582                                       (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
    583                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    584                                       (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
    585                                        (loop (cons* s0 s1 lst) (cdr tlst)))
    586                                       (else
    587                                        (nemo:error 'nemo:state-init-eq ": invalid transition equation "
    588                                                    (car tlst) " in state complex " n))
    589                                       (else (loop lst (cdr tlst)))))))
    590          (state-subs     (fold (lambda (s ax) (subst-extend s (matlab-state-name n s) ax)) subst-empty state-list))
    591          (init1          (map (lambda (lineq) (match lineq ((i '= . expr) `(,i = . ,(subst-convert expr state-subs)))))
    592                               init)))
    593     (list (matlab-name n) init1)))
    594436
    595437
     
    619461
    620462
    621 (define (poset->state-eq-defs poset sys kinetic)
     463(define (poset->state-eq-defs poset sys)
    622464  (fold-right
    623465   (lambda (lst ax)
     
    625467              (match-let (((i . n)  x))
    626468                         (let ((en (environment-ref sys n)))
    627                            (if (and (not (member n kinetic)) (nemo:quantity? en))
     469                           (if (nemo:quantity? en)
    628470                               (cases nemo:quantity en
    629                                       (TSCOMP  (name initial open transitions power)
     471                                      (TSCOMP  (name initial open transitions conserve power)
    630472                                               (append (state-eqs name initial open transitions power) ax))
    631473                                      (else  ax))
     
    635477
    636478
    637 (define (poset->kstate-eq-defs poset sys kinetic)
    638   (fold-right
    639    (lambda (lst ax)
    640      (fold  (lambda (x ax)
    641               (match-let (((i . n)  x))
    642                          (let ((en (environment-ref sys n)))
    643                            (if (and (member n kinetic) (nemo:quantity? en))
    644                                (cases nemo:quantity en
    645                                       (TSCOMP  (name initial open transitions power)
    646                                                (append (kstate-eqs name initial open transitions power) ax))
    647                                       (else  ax))
    648                                ax))))
    649             ax lst))
    650    (list) poset))
    651479
    652480
     
    659487                           (if (nemo:quantity? en)
    660488                               (cases nemo:quantity en
    661                                       (TSCOMP  (name initial open transitions power)
     489                                      (TSCOMP  (name initial open transitions conserve power)
    662490                                               (cons (stcomp-eq name open transitions) ax))
    663491                                      (else  ax))
     
    675503                           (if (nemo:quantity? en)
    676504                               (cases nemo:quantity en
    677                                       (TSCOMP  (name initial open transitions power)
     505                                      (TSCOMP  (name initial open transitions conserve power)
    678506                                               (if (nemo:rhs? initial)
    679507                                                   (cons* (state-init name initial)
     
    686514
    687515
    688 (define (poset->state-init-eq-defs poset sys)
    689   (fold-right
    690    (lambda (lst ax)
    691      (fold  (lambda (x ax)
    692               (match-let (((i . n)  x))
    693                          (let ((en (environment-ref sys n)))
    694                            (if (nemo:quantity? en)
    695                                (cases nemo:quantity en
    696                                       (TSCOMP (name initial open transitions power)
    697                                               (if (and (list? initial) (every nemo:lineq? initial))
    698                                                   (cons (state-init-eq name transitions initial) ax)
    699                                                   ax))
    700                                       (else  ax))
    701                                ax))))
    702             ax lst))
    703    (list) poset))
    704 
    705 
    706516(define (find-locals defs)
    707517  (concatenate (map (lambda (def) (match def (('let bnds _) (map first bnds)) (else (list)))) defs)))
     
    712522    (if (nemo:quantity? en)
    713523        (cases nemo:quantity en
    714                (TSCOMP  (name initial open transitions power)  power)
     524               (TSCOMP  (name initial open transitions conserve power)  power)
    715525               (else  #f))  #f)))
    716526
     
    746556  (define (cid x)  (second x))
    747557  (define (cn x)   (first x))
    748   (let-optionals rest ((method 'cnexp) (table? #f) (min-v -100) (max-v 100) (step 0.5)
    749                        (depend #f)  (kinetic (list)) )
     558  (let-optionals rest ((method 'cnexp) (table? #f) (min-v -100) (max-v 100) (step 0.5) )
    750559  (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
    751560    (let ((imports  ((dis 'imports)  sys))
     
    767576
    768577        (match-let (((state-list asgn-list g) deps*))
    769          (let* ((poset          (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
    770                 (asgn-eq-defs   (poset->asgn-eq-defs poset sys))
     578         (let* ((poset            (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
     579                (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
    771580                (perm-ions (fold (lambda (ionch ax)
    772581                                    (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
     
    796605                                  (let ((ion (car ep)))
    797606                                    `(,(matlab-name ion) ,(matlab-name (s+ 'i ion)) ,(matlab-name (s+ ion 'i)))))
    798                                epools))
    799                (has-kinetic? (or (not (null? (filter (lambda (x) (member (car x) kinetic)) states)))))
    800                (has-ode?     (or (not (null? (filter (lambda (x) (not (member (car x) kinetic))) states)))
    801                                  (not (null? pool-ions)))))
    802 
     607                               epools)))
    803608
    804609           (for-each
     
    819624
    820625           (let* ((with      (inexact->exact (round (/ (abs (- max-v min-v)) step))))
    821                   (define-fn (make-define-fn table? min-v max-v with depend)))
     626                  (define-fn (make-define-fn table? min-v max-v with)))
    822627             (for-each (lambda (fndef)
    823628                         (if (not (member (car fndef) builtin-fns))
     
    825630                       defuns))
    826631
    827            (pp indent ,nl (dy = ,sysname (,(sl\ ", " `(dy t)))))
     632           (pp indent ,nl (function dy = ,sysname (,(sl\ ", " `(dy t)))))
    828633           
    829634           (if (not (null? exports)) (pp indent+ (global ,(sl\ ", " exports))))
     
    871676                             (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs))
    872677
    873 #|
    874678           (let* ((i-eqs (filter-map
    875679                          (lambda (ionch)
     
    924728                                      )))
    925729                          ionchs))
     730
    926731                  (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
     732
    927733                  (i-eqs  (fold (lambda (b ax)
    928734                                  (match b
     
    942748                                         (else ax)))
    943749                                (list) i-bkts))
    944                   (locals (find-locals (map second i-eqs))))
    945              (if (not (null? locals))    (pp indent+ (LOCAL ,(sl\ ", " locals))))
    946              (if (not (null? asgns))     (pp indent+ (rates ())))
    947              (if has-ode?
    948                  (if (not method) (pp indent+ (SOLVE states))
    949                      (pp indent+ (SOLVE states METHOD ,method))))
    950              (if has-kinetic?   (pp indent+  (SOLVE kstates METHOD sparse)))
    951              (if (not (null? stcomps))   (pp indent+ (stcomps ())))
    952              (if (not (null? pool-ions)) (pp indent+ (pools ())))
    953              (for-each (lambda (p) (pp indent+ ,(expr->string/MATLAB (second p) (first p)))) i-eqs)
    954              (pp indent "}"))
    955            
    956            (if has-ode?
    957                (begin
    958                  (pp indent ,nl (DERIVATIVE states "{"))
    959                  (let* ((state-eq-defs  (reverse (poset->state-eq-defs poset sys kinetic)))
    960                         (pool-eq-defs
    961                          (map (lambda (ep)
    962                                 (let* ((ep-name     (first ep))
    963                                        (pool-ion    (assoc ep-name pool-ions))
    964                                        (i-name      (second pool-ion))
    965                                        (init-name   (matlab-name (s+ ep-name '-init)))
    966                                        (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
    967                                        (beta-name   (matlab-name (s+ ep-name '-beta)))
    968                                        (depth-name  (matlab-name (s+ ep-name '-depth)))
    969                                        (rhs         `(let ((F 96485.0))
    970                                                        (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name))
    971                                                           (* ,beta-name ,ep-name .
    972                                                              ,(if temp-name (list temp-name) (list)))))))
    973                                   `(,(s+ ep-name "'") ,rhs)))
    974                               epools))
    975                         (eq-defs (append pool-eq-defs state-eq-defs))
    976                         (locals  (find-locals (map second eq-defs))))
    977                    (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
    978                    (for-each (lambda (def)
    979                                (let ((n (first def)) (b (second def)))
    980                                  (pp indent+ ,(expr->string/MATLAB b n)))) eq-defs))
    981                    
    982                  (pp indent "}")))
    983            
    984 |#
    985 
    986 )))))))
     750
     751                  (state-eq-defs    (reverse (poset->state-eq-defs poset sys)))
     752
     753                  (stcomp-eq-defs   (poset->stcomp-eq-defs poset sys))
     754
     755                  (pool-eq-defs
     756                   (map (lambda (ep)
     757                          (let* ((ep-name     (first ep))
     758                                 (pool-ion    (assoc ep-name pool-ions))
     759                                 (i-name      (second pool-ion))
     760                                 (init-name   (matlab-name (s+ ep-name '-init)))
     761                                 (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
     762                                 (beta-name   (matlab-name (s+ ep-name '-beta)))
     763                                 (depth-name  (matlab-name (s+ ep-name '-depth)))
     764                                 (rhs         `(let ((F 96485.0))
     765                                                 (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name))
     766                                                    (* ,beta-name ,ep-name .
     767                                                       ,(if temp-name (list temp-name) (list)))))))
     768                            `(,(s+ ep-name)  ,rhs)))
     769                        epools))
     770
     771                  (rate-eq-defs (append pool-eq-defs state-eq-defs))
     772                 
     773                  (state-index-map  (let ((acc (fold (lambda (def ax)
     774                                                       (let ((st-name (first def)))
     775                                                         (list (+ 1 (first ax))
     776                                                               (cons `(,st-name ,(first ax)) (second ax)))))
     777                                                     (list 0 (list)) rate-eq-defs)))
     778                                    (second acc)))
     779                  )
     780
     781
     782             (pp indent+ (dy = zeros(length(y),1)))
     783
     784             (for-each (lambda (def)
     785                         (let ((n (first def)) )
     786                           (pp indent+ (,n = ,(s+ "y(" (lookup-def n state-index-map) ")"))))) rate-eq-defs)
     787
     788             (for-each (lambda (def)
     789                         (let ((n (first def)) (b (second def)))
     790                           (pp indent+ ,(expr->string/MATLAB b n)))) stcomp-eq-defs)
     791
     792             (for-each (lambda (def)
     793                         (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) (b (second def)))
     794                           (pp indent+ ,(expr->string/MATLAB b n)))) rate-eq-defs)
     795
     796             (for-each (lambda (pool-ion)
     797                         (pp indent+ (,(third pool-ion) = ,(first pool-ion)))) pool-ions)
     798
     799             (for-each (lambda (def)
     800                         (pp indent+ ,(expr->string/MATLAB (second def) (first def)))) i-eqs)
     801             
     802             ))))))))
Note: See TracChangeset for help on using the changeset viewer.