Changeset 12554 in project
- Timestamp:
- 11/18/08 08:31:57 (12 years ago)
- Location:
- release/3/nemo/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
release/3/nemo/trunk/nemo-matlab.scm
r12353 r12554 99 99 (define (s+ . lst) (string-concatenate (map ->string lst))) 100 100 (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))102 101 (define (sl\ p lst) (string-intersperse (map ->string lst) p)) 103 102 (define nl "\n") … … 337 336 (sdoc->string (doc:format width (format-expr/MATLAB 2 x rv))))) 338 337 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) 422 339 (lambda (indent n proc) 423 340 (let ((lst (procedure-data proc)) … … 431 348 (let* ((body1 (canonicalize-expr/MATLAB (rhsexpr body))) 432 349 (lbs (enum-bnds body1 (list)))) 433 ; (if (and table? min-v max-v with)434 ; (match vars435 ; (('v) (pp indent+ (TABLE ,@(if depend `(DEPEND ,depend) `(""))436 ; FROM ,min-v TO ,max-v WITH ,with)))437 ; (else (void))))438 350 (pp indent+ ,(expr->string/MATLAB body1 retval)) 439 351 (pp indent endfunction)) … … 510 422 (fbody (rhsexpr rhs1)) 511 423 (fbody1 (canonicalize-expr/MATLAB fbody))) 512 (cons (list (s+ name "'")fbody1) ax)))))424 (cons (list name fbody1) ax))))) 513 425 (list) nodes))) 514 426 eqs)))) … … 516 428 517 429 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 (else532 (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 system537 (map538 (lambda (e)539 (match e540 (('-> 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 566 430 567 431 (define (state-init n init) … … 570 434 (list (matlab-name n) init1))) 571 435 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 (else587 (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)))594 436 595 437 … … 619 461 620 462 621 (define (poset->state-eq-defs poset sys kinetic)463 (define (poset->state-eq-defs poset sys) 622 464 (fold-right 623 465 (lambda (lst ax) … … 625 467 (match-let (((i . n) x)) 626 468 (let ((en (environment-ref sys n))) 627 (if ( and (not (member n kinetic)) (nemo:quantity? en))469 (if (nemo:quantity? en) 628 470 (cases nemo:quantity en 629 (TSCOMP (name initial open transitions power)471 (TSCOMP (name initial open transitions conserve power) 630 472 (append (state-eqs name initial open transitions power) ax)) 631 473 (else ax)) … … 635 477 636 478 637 (define (poset->kstate-eq-defs poset sys kinetic)638 (fold-right639 (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 en645 (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))651 479 652 480 … … 659 487 (if (nemo:quantity? en) 660 488 (cases nemo:quantity en 661 (TSCOMP (name initial open transitions power)489 (TSCOMP (name initial open transitions conserve power) 662 490 (cons (stcomp-eq name open transitions) ax)) 663 491 (else ax)) … … 675 503 (if (nemo:quantity? en) 676 504 (cases nemo:quantity en 677 (TSCOMP (name initial open transitions power)505 (TSCOMP (name initial open transitions conserve power) 678 506 (if (nemo:rhs? initial) 679 507 (cons* (state-init name initial) … … 686 514 687 515 688 (define (poset->state-init-eq-defs poset sys)689 (fold-right690 (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 en696 (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 706 516 (define (find-locals defs) 707 517 (concatenate (map (lambda (def) (match def (('let bnds _) (map first bnds)) (else (list)))) defs))) … … 712 522 (if (nemo:quantity? en) 713 523 (cases nemo:quantity en 714 (TSCOMP (name initial open transitions power) power)524 (TSCOMP (name initial open transitions conserve power) power) 715 525 (else #f)) #f))) 716 526 … … 746 556 (define (cid x) (second x)) 747 557 (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) ) 750 559 (match-let ((($ nemo:quantity 'DISPATCH dis) (environment-ref sys (nemo-intern 'dispatch)))) 751 560 (let ((imports ((dis 'imports) sys)) … … 767 576 768 577 (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)) 771 580 (perm-ions (fold (lambda (ionch ax) 772 581 (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch))) … … 796 605 (let ((ion (car ep))) 797 606 `(,(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))) 803 608 804 609 (for-each … … 819 624 820 625 (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))) 822 627 (for-each (lambda (fndef) 823 628 (if (not (member (car fndef) builtin-fns)) … … 825 630 defuns)) 826 631 827 (pp indent ,nl ( dy = ,sysname (,(sl\ ", " `(dy t)))))632 (pp indent ,nl (function dy = ,sysname (,(sl\ ", " `(dy t))))) 828 633 829 634 (if (not (null? exports)) (pp indent+ (global ,(sl\ ", " exports)))) … … 871 676 (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)) 872 677 873 #|874 678 (let* ((i-eqs (filter-map 875 679 (lambda (ionch) … … 924 728 ))) 925 729 ionchs)) 730 926 731 (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs)) 732 927 733 (i-eqs (fold (lambda (b ax) 928 734 (match b … … 942 748 (else ax))) 943 749 (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 )))))))) -
release/3/nemo/trunk/nemo.scm
r12469 r12554 200 200 (lambda (options model) 201 201 (nemo:matlab-translator model (lookup-def 'method options) (lookup-def 'table options) 202 -150 150 1 #f 203 (lookup-def 'kinetic options) )) 202 -150 150 1 )) 204 203 (lambda (options model) 205 204 (void)))) … … 665 664 (model->matlab `(;(method . ,matlab-method) 666 665 (table . ,(assoc 't options)) 667 (kinetic 668 ,(let ((opt (assoc 'matlab-kinetic options))) 669 (and opt 670 (if (string? (cdr opt)) 671 (map string->symbol (string-split (cdr opt) ",")) 672 'all))))) 673 model)))) 666 ) 667 model)))) 674 668 )) 675 669 operands)))
Note: See TracChangeset
for help on using the changeset viewer.