Changeset 15521 in project


Ignore:
Timestamp:
08/20/09 04:05:07 (10 years ago)
Author:
Ivan Raikov
Message:

Fixes related to symbolic differentiation in nemo

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

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/examples/carelli05_run.m

    r15460 r15521  
    1717TOL = 0.00001;
    1818
     19N = length(y_initial);
     20
    1921x(1) = t;
    2022nstep = ( t(2) - t(1) ) / h;
     
    3133  while ( resmax > TOL & it < 10 )
    3234    r = yp - ( y(:,i) + h * feval ( f, x(i+1), yp ) );
    33     rprime = 1.0 - h * feval ( fy, x(i+1), yp );
     35    rprime = eye(N,N) - h * feval ( fy, x(i+1), yp );
    3436    yp = yp - r / rprime;
    3537    resmax = max ( r );
  • release/4/nemo/trunk/nemo-matlab.scm

    r15460 r15521  
    243243
    244244    res))
    245  
    246 (define (make-define-fn table? min-v max-v with)
    247   (lambda (indent n proc)
    248     (let ((lst (procedure-data proc))
    249           (indent+ (+ 2 indent)))
    250       (let ((retval   (matlab-name (gensym "retval")))
    251             (rt       (lookup-def 'rt lst))
    252             (formals  (lookup-def 'formals lst))
    253             (vars     (lookup-def 'vars lst))
    254             (body     (lookup-def 'body lst)))
    255         (pp indent ,nl (function ,retval = ,(matlab-name n) (,(sl\ ", " vars)) ))
    256         (let* ((body1 (canonicalize-expr/MATLAB (rhsexpr/MATLAB body)))
    257                (lbs   (enum-bnds body1 (list))))
    258           (pp indent+ ,(expr->string/MATLAB body1 retval))
    259           (pp indent endfunction))
    260           ))))
    261245
    262246
     
    433417                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
    434418
     419 
     420(define (make-define-fn table? min-v max-v with)
     421  (lambda (indent n proc)
     422    (let ((lst (procedure-data proc))
     423          (indent+ (+ 2 indent)))
     424      (let ((retval   (matlab-name (gensym "retval")))
     425            (rt       (lookup-def 'rt lst))
     426            (formals  (lookup-def 'formals lst))
     427            (vars     (lookup-def 'vars lst))
     428            (body     (lookup-def 'body lst)))
     429        (pp indent ,nl (function ,retval = ,(matlab-name n) (,(sl\ ", " vars)) ))
     430        (let* ((body1 (canonicalize-expr/MATLAB (rhsexpr/MATLAB body)))
     431               (lbs   (enum-bnds body1 (list))))
     432          (pp indent+ ,(expr->string/MATLAB body1 retval))
     433          (pp indent endfunction))
     434          ))))
     435
     436
     437(define (output-dy sysname method globals state-index-map
     438                   rate-eq-defs reaction-eq-defs asgn-eq-defs
     439                   pool-ions mcap i-eqs indent indent+)
     440
     441  (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
     442                                                       ((lsode)  `(y t))
     443                                                       ((odepkg) `(t y))
     444                                                       (else     `(y t)))))))
     445 
     446  (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     447           
     448  (if (member 'v globals)
     449      (let ((vi (lookup-def 'v state-index-map)))
     450        (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
     451           
     452  (for-each (lambda (def)
     453              (let ((n (first def)) )
     454                (pp indent+ ,(expr->string/MATLAB
     455                              (s+ "y(" (lookup-def n state-index-map) ")") n))))
     456            rate-eq-defs)
     457 
     458  (for-each (lambda (def)
     459              (let ((n (matlab-name (first def)) )
     460                    (b (second def)))
     461                (pp indent+ ,(expr->string/MATLAB b n))))
     462            asgn-eq-defs)
     463 
     464  (for-each (lambda (def)
     465              (let ((n (first def)) (b (second def)))
     466                (pp indent+ ,(expr->string/MATLAB b n))))
     467            reaction-eq-defs)
     468 
     469  (for-each (lambda (pool-ion)
     470              (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
     471            pool-ions)
     472 
     473  (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'dy))
     474 
     475  (for-each (lambda (def)
     476              (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")"))
     477                    (b (second def)))
     478                (pp indent+ ,(s+ "# state " (first def))
     479                    ,(expr->string/MATLAB b n))))
     480            rate-eq-defs)
     481 
     482  (for-each (lambda (def)
     483              (pp indent+ ,(expr->string/MATLAB (second def) (first def))))
     484            i-eqs)
     485 
     486  (if (and mcap (member 'v globals))
     487      (let ((vi (lookup-def 'v state-index-map)))
     488        (if vi (pp indent+  ,(expr->string/MATLAB `(/ (neg ,(sum (map first i-eqs))) ,mcap) 
     489                                                  (s+ "dy(" vi ")"))))))
     490 
     491  (pp indent ,nl (endfunction))
     492)
     493
     494 
     495(define (make-define-dfn fenv)
     496
     497  (define (D vars body)
     498    (let loop ((vars vars) (exprs (list)))
     499      (if (null? vars) (if (pair? (cdr exprs))
     500                           `(+ . ,(reverse exprs))  (car exprs))
     501          (let ((de (let ((de (differentiate fenv (car vars) body)))
     502                       (let loop ((e (simplify de)) (de de))
     503                         (if (equal? de e) de
     504                             (loop (simplify e) e))))))
     505            (loop (cdr vars) (cons de exprs))))))
     506
     507  (lambda (indent n proc)
     508    (let ((lst (procedure-data proc))
     509          (indent+ (+ 2 indent)))
     510      (let ((retval   (matlab-name (gensym "retval")))
     511            (rt       (lookup-def 'rt lst))
     512            (formals  (lookup-def 'formals lst))
     513            (vars     (lookup-def 'vars lst))
     514            (body     (lookup-def 'body lst)))
     515        (pp indent ,nl (function ,retval = ,(s+ "d_" (matlab-name n) ) (,(sl\ ", " vars)) ))
     516        (let* ((dbody (D vars body))
     517               (dbody (canonicalize-expr/MATLAB (rhsexpr/MATLAB dbody))))
     518          (pp indent+ ,(expr->string/MATLAB dbody retval))
     519          (pp indent endfunction))
     520          )))
     521)
     522
     523(define (output-ddy sysname method globals fenv state-index-map
     524                    rate-eq-defs reaction-eq-defs asgn-eq-defs
     525                    defuns indent indent+)
     526
     527
     528  (define (D var expr)
     529    (let ((de (differentiate fenv var expr)))
     530      (let loop ((e (simplify de)) (de de))
     531        (if (equal? de e) de
     532            (loop (simplify e) e)))))
     533   
     534  (pp indent ,nl (function ddy = ,(s+ sysname "_ddy")
     535                           (,(sl\ ", " `(t y)))))
     536 
     537  (if (member 'v globals)
     538      (let ((vi (lookup-def 'v state-index-map)))
     539        (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
     540 
     541  (for-each (lambda (def)
     542              (let ((name (first def)))
     543                (pp indent+ ,(expr->string/MATLAB
     544                              (s+ "y(" (lookup-def name state-index-map) ")") name))))
     545            rate-eq-defs)
     546 
     547  (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'ddy))
     548
     549  (for-each (lambda (def)
     550              (let* ((name  (first def))
     551                     (v     (s+ "ddy(" (lookup-def name state-index-map) ")"))
     552                     (b     (second def))
     553                     (fv    (enum-freevars b '() '() ))
     554                     (db    (D name
     555                               `(let ,(filter-map (lambda (x) (assoc x asgn-eq-defs))
     556                                                  fv)
     557                                  ,b)))
     558                     (cdb    (canonicalize-expr/MATLAB (rhsexpr/MATLAB db))))
     559                (pp indent+ ,(s+ "# state " name)
     560                    ,(expr->string/MATLAB cdb v))))
     561            rate-eq-defs)
     562 
     563  (pp indent ,nl (endfunction))
     564)
     565
     566
     567(define (output-steadystate  sysname method globals steady-state-index-map
     568                             rate-eq-defs reaction-eq-defs asgn-eq-defs
     569                             indent indent+)
     570 
     571  (if (not (null? steady-state-index-map))
     572      (begin
     573        (pp indent ,nl (function y = ,(s+ sysname "_steadystate") (x)))
     574       
     575        (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     576       
     577        (for-each
     578         (lambda (def)
     579           (let* ((n   (first def))
     580                  (ni  (lookup-def n steady-state-index-map)))
     581             (if ni (pp indent+ ,(expr->string/MATLAB (s+ "x(" ni ")") n)))
     582             ))
     583         rate-eq-defs)
     584       
     585        (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'y))
     586       
     587        (for-each (lambda (def)
     588                    (let ((n (matlab-name (first def)) )
     589                          (b (second def)))
     590                      (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)
     591       
     592        (for-each
     593         (lambda (def)
     594           (let* ((n   (first def))
     595                  (ni  (lookup-def n steady-state-index-map))
     596                  (b   (second def)))
     597             (if ni (pp indent+ ,(s+ "# state " n)
     598                        ,(expr->string/MATLAB b (s+ "y(" ni ")"))))
     599             ))
     600         rate-eq-defs)
     601       
     602        (pp indent ,nl (endfunction))))
     603)
     604
     605
     606(define (output-print-state sysname state-index-map indent indent+)
     607  (pp indent ,nl (function  ,(s+ sysname "_print_state") (y)))
     608 
     609  (let ((lst (sort (map (lambda (x) (cons (->string (car x)) (cdr x))) state-index-map)
     610                   (lambda (x y) (string<? (car x) (car y))))))
     611    (for-each (lambda (p)
     612                (let ((n (first p)) (i (second p)))
     613                  (pp indent+ (,n " = " "y(" ,i ")"))))
     614              lst))
     615 
     616  (pp indent ,nl (endfunction))
     617)
     618
     619
     620(define (output-init sysname globals state-index-map steady-state-index-map
     621                     const-defs init-eq-defs rate-eq-defs i-eqs pool-ions
     622                     pool-ion-i  indent indent+)
     623 
     624  (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
     625 
     626  (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
     627 
     628  (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
     629 
     630  (if (member 'v globals)
     631      (let ((vi (lookup-def 'v state-index-map)))
     632        (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
     633        (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
     634 
     635  (for-each (lambda (def)
     636              (let ((n (first def)) (b (second def)))
     637                (pp indent+ ,(expr->string/MATLAB b n))))
     638            const-defs)
     639 
     640  (for-each (lambda (def)
     641              (let* ((n  (first def))
     642                     (ni (lookup-def n state-index-map))
     643                     (b (second def)))
     644                (pp indent+ ,(expr->string/MATLAB b n))
     645                (if ni (pp indent+ ,(expr->string/MATLAB n  (s+ "y0(" ni ")"))))))
     646            init-eq-defs)
     647 
     648  (for-each (lambda (pool-ion)
     649              (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
     650            pool-ions)
     651 
     652  (if (not (null? steady-state-index-map))
     653      (begin
     654        (for-each (lambda (def)
     655                    (if (member (first def) pool-ion-i)
     656                        (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
     657                  i-eqs)
     658       
     659        (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'xs)
     660            ,(expr->string/MATLAB
     661              `(fsolve ,(s+ #\" sysname "_steadystate" #\" ) xs)
     662              "[ys,fval,info]"))
     663       
     664       
     665        (for-each
     666         (lambda (def)
     667           (let* ((n (first def))
     668                  (si (lookup-def n steady-state-index-map))
     669                  (ni (lookup-def n state-index-map)))
     670             (if si (pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")"))))))
     671         rate-eq-defs)
     672        ))
     673 
     674 
     675  (for-each (lambda (def)
     676              (if (not (member (first def) pool-ion-i))
     677                  (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
     678            i-eqs)
     679 
     680  (pp indent ,nl (endfunction))
     681)
    435682
    436683(define (nemo:matlab-translator sys . rest)
     
    565812             (pool-ion-i          (map second pool-ions))
    566813             
    567              
    568              (state-index-map  (let ((acc (fold (lambda (def ax)
    569                                                   (let ((st-name (first def)))
    570                                                     (list (+ 1 (first ax))
    571                                                           (cons `(,st-name ,(first ax)) (second ax)))))
    572                                                 (list 1 (list))
    573                                                 (cons (list 'v) rate-eq-defs))))
    574                                  
    575                                  (second acc)))
    576              
    577              (steady-state-index-map  (let ((acc (fold (lambda (def ax)
    578                                                          (let ((st-name (first def)))
    579                                                            (if (not (alist-ref st-name init-eq-defs))
    580                                                                (list (+ 1 (first ax))
    581                                                                      (cons `(,st-name ,(first ax)) (second ax)))
    582                                                                ax)))
    583                                                        (list 1 (list))
    584                                                        rate-eq-defs)))
    585                                         (second acc)))
    586              
    587814             (globals          (map matlab-name
    588815                                    (delete-duplicates (append
     
    597824                                                        (map first   imports)
    598825                                                        (map first   const-defs)))))
     826             
     827             (state-index-map  (let ((acc (fold (lambda (def ax)
     828                                                  (let ((st-name (first def)))
     829                                                    (list (+ 1 (first ax))
     830                                                          (cons `(,st-name ,(first ax)) (second ax)))))
     831                                                (list 1 (list))
     832                                                (if (member 'v globals)
     833                                                    (cons (list 'v) rate-eq-defs)
     834                                                    rate-eq-defs)
     835                                                )))
     836                                 (second acc)))
     837             
     838             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
     839                                                         (let ((st-name (first def)))
     840                                                           (if (not (alist-ref st-name init-eq-defs))
     841                                                               (list (+ 1 (first ax))
     842                                                                     (cons `(,st-name ,(first ax)) (second ax)))
     843                                                               ax)))
     844                                                       (list 1 (list))
     845                                                       rate-eq-defs)))
     846                                        (second acc)))
     847             
     848             (dfenv
     849              (map (lambda (x) (let ((n (first x)))
     850                                 (list n (matlab-name (s+ "d_" n )))))
     851                   defuns))
    599852
    600853               )
     
    618871           (if (not (null? globals)) (pp indent (global ,(sl\ " " globals))))
    619872
    620            (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
    621                                                                 ((lsode)  `(y t))
    622                                                                 ((odepkg) `(t y))
    623                                                                 (else     `(y t)))))))
    624 
    625            (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
    626            
    627            (if (member 'v globals)
    628                (let ((vi (lookup-def 'v state-index-map)))
    629                  (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
    630 
    631            
    632            (for-each (lambda (def)
    633                        (let ((n (first def)) )
    634                          (pp indent+ ,(expr->string/MATLAB
    635                                        (s+ "y(" (lookup-def n state-index-map) ")") n))))
    636                      rate-eq-defs)
    637 
    638            (for-each (lambda (def)
    639                        (let ((n (matlab-name (first def)) )
    640                              (b (second def)))
    641                          (pp indent+ ,(expr->string/MATLAB b n))))
    642                      asgn-eq-defs)
    643            
    644            (for-each (lambda (def)
    645                        (let ((n (first def)) (b (second def)))
    646                          (pp indent+ ,(expr->string/MATLAB b n))))
    647                      reaction-eq-defs)
    648 
    649            (for-each (lambda (pool-ion)
    650                        (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
    651                      pool-ions)
    652 
    653          (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'dy))
     873           ;; derivative function
     874           (output-dy sysname method globals state-index-map
     875                      rate-eq-defs reaction-eq-defs asgn-eq-defs
     876                      pool-ions mcap i-eqs
     877                      indent indent+)
     878
     879           ;; partial derivative
     880           (output-ddy sysname method globals dfenv state-index-map
     881                       rate-eq-defs reaction-eq-defs asgn-eq-defs
     882                       defuns indent indent+)
     883
     884           ;; steady state solver
     885           (output-steadystate sysname method globals steady-state-index-map
     886                               rate-eq-defs reaction-eq-defs asgn-eq-defs
     887                               indent indent+)
     888
     889           ;; print state vector
     890           (output-print-state sysname state-index-map
     891                               indent indent+)
     892
     893           ;; initial values function
     894           (output-init sysname globals state-index-map steady-state-index-map
     895                        const-defs init-eq-defs rate-eq-defs i-eqs pool-ions
     896                        pool-ion-i indent indent+)
     897
     898           (pp indent ,nl)
    654899         
    655 
    656          (for-each (lambda (def)
    657                      (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")"))
    658                            (b (second def)))
    659                        (pp indent+ ,(s+ "# state " (first def))
    660                            ,(expr->string/MATLAB b n))))
    661                    rate-eq-defs)
     900           ;; user-defined functions
     901           (let* ((with       (inexact->exact (round (/ (abs (- max-v min-v)) step))))
     902                  (define-fn  (make-define-fn table? min-v max-v with)))
     903             (for-each (lambda (fndef)
     904                         (if (not (member (car fndef) builtin-fns))
     905                             (apply define-fn (cons indent fndef))))
     906                       defuns))
    662907         
    663          (for-each (lambda (def)
    664                      (pp indent+ ,(expr->string/MATLAB (second def) (first def))))
    665                    i-eqs)
    666          
    667          (if (and mcap (member 'v globals))
    668              (let ((vi (lookup-def 'v state-index-map)))
    669                (if vi (pp indent+  ,(expr->string/MATLAB `(/ (neg ,(sum (map first i-eqs))) ,mcap) 
    670                                                          (s+ "dy(" vi ")"))))))
    671          
    672          (pp indent ,nl (endfunction))
    673 
    674          (pp indent ,nl (function dy = ,(s+ sysname "_dy")
    675                                   (,(sl\ ", " `(t y)))))
    676 
    677          
    678 
    679          (pp indent ,nl (endfunction))
    680 
    681 
    682          (if (not (null? steady-state-index-map))
    683              (begin
    684                (pp indent ,nl (function y = ,(s+ sysname "_steadystate") (x)))
    685                
    686                (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
    687                
    688                (for-each
    689                 (lambda (def)
    690                   (let* ((n   (first def))
    691                          (ni  (lookup-def n steady-state-index-map)))
    692                     (if ni (pp indent+ ,(expr->string/MATLAB (s+ "x(" ni ")") n)))
    693                     ))
    694                 rate-eq-defs)
    695                
    696                (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'y))
    697                
    698                (for-each (lambda (def)
    699                            (let ((n (matlab-name (first def)) )
    700                                  (b (second def)))
    701                              (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)
    702                
    703                (for-each
    704                 (lambda (def)
    705                   (let* ((n   (first def))
    706                          (ni  (lookup-def n steady-state-index-map))
    707                          (b   (second def)))
    708                     (if ni (pp indent+ ,(s+ "# state " n)
    709                                ,(expr->string/MATLAB b (s+ "y(" ni ")"))))
    710                     ))
    711                 rate-eq-defs)
    712                
    713                (pp indent ,nl (endfunction))))
    714              
    715 
    716          (pp indent ,nl (function  ,(s+ sysname "_print_state") (y)))
    717 
    718          (let ((lst (sort (map (lambda (x) (cons (->string (car x)) (cdr x))) state-index-map)
    719                           (lambda (x y) (string<? (car x) (car y))))))
    720            (for-each (lambda (p)
    721                        (let ((n (first p)) (i (second p)))
    722                          (pp indent+ (,n " = " "y(" ,i ")"))))
    723                      lst))
    724 
    725          (pp indent ,nl (endfunction))
    726 
    727          
    728          (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
    729          
    730          (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
    731          
    732          (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
    733          
    734          (if (member 'v globals)
    735              (let ((vi (lookup-def 'v state-index-map)))
    736                (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
    737                (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
    738          
    739          (for-each (lambda (def)
    740                      (let ((n (first def)) (b (second def)))
    741                        (pp indent+ ,(expr->string/MATLAB b n))))
    742                    const-defs)
    743          
    744          (for-each (lambda (def)
    745                      (let* ((n  (first def))
    746                             (ni (lookup-def n state-index-map))
    747                             (b (second def)))
    748                        (pp indent+ ,(expr->string/MATLAB b n))
    749                        (if ni (pp indent+ ,(expr->string/MATLAB n  (s+ "y0(" ni ")"))))))
    750                    init-eq-defs)
    751          
    752          (for-each (lambda (pool-ion)
    753                      (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
    754                    pool-ions)
    755 
    756          (if (not (null? steady-state-index-map))
    757              (begin
    758                (for-each (lambda (def)
    759                                  (if (member (first def) pool-ion-i)
    760                                      (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
    761                          i-eqs)
    762 
    763                (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'xs)
    764                    ,(expr->string/MATLAB
    765                      `(fsolve ,(s+ #\" sysname "_steadystate" #\" ) xs)
    766                      "[ys,fval,info]"))
    767                
    768 
    769                (for-each
    770                 (lambda (def)
    771                   (let* ((n (first def))
    772                          (si (lookup-def n steady-state-index-map))
    773                          (ni (lookup-def n state-index-map)))
    774                     (if si (pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")"))))))
    775                 rate-eq-defs)
    776                ))
    777 
    778 
    779          (for-each (lambda (def)
    780                      (if (not (member (first def) pool-ion-i))
    781                          (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
    782                    i-eqs)
    783 
    784          (pp indent ,nl (endfunction))
    785          
    786          (pp indent ,nl)
    787          
    788          (let* ((with       (inexact->exact (round (/ (abs (- max-v min-v)) step))))
    789                 (define-fn  (make-define-fn table? min-v max-v with)))
    790            (for-each (lambda (fndef)
    791                        (if (not (member (car fndef) builtin-fns))
    792                            (apply define-fn (cons indent fndef))))
    793                      defuns))
    794          )))))
     908           ;; derivatives of user-defined functions
     909           (let ((define-dfn  (make-define-dfn dfenv)))
     910             (for-each (lambda (fndef)
     911                         (if (not (member (car fndef) builtin-fns))
     912                             (apply define-dfn (cons indent fndef))))
     913                       defuns))
     914
     915
     916           )))))
    795917)
  • release/4/nemo/trunk/nemo-utils.scm

    r15460 r15521  
    112112(define (let-lift expr)
    113113  (let ((bnds (let-enum expr (list))))
    114     (if (null? bnds) expr
     114    (if (null? bnds) (let-elim expr)
    115115        `(let ,(map (lambda (b) (list (car b) (let-elim (cadr b)))) bnds) ,(let-elim expr)))))
    116116
     
    221221(define LOG2E  1.44269504088896)
    222222
    223 (define (differentiate fenv xenv x t)
     223
     224(define (differentiate fenv x t)
     225  (define subst-convert
     226    (subst-driver
     227     (lambda (x) (and (symbol? x) x))
     228     nemo:binding? identity nemo:bind nemo:subst-term))
     229 
    224230  (cond ((number? t)  0.0)
    225231        ((symbol? t)  (cond ((equal? x t) 1.0)
    226                             ((member t xenv) t)
    227232                            (else 0.0)))
    228233        (else (match t
    229                 (('neg u)  `(neg ,(differentiate fenv xenv x u)))
    230 
    231                 (('+ u v)  `(+ ,(differentiate fenv xenv x u) ,(differentiate fenv xenv x v)))
    232                 (('- u v)  `(- ,(differentiate fenv xenv x u) ,(differentiate fenv xenv x v)))
    233 
    234                 (('* (and u (? number?)) v)    `(* ,u ,(differentiate fenv xenv x v)))
    235                 (('* v (and u (? number?)))    `(* ,u ,(differentiate fenv xenv x v)))
    236 
    237                 (('* u v)     `(+ (* ,(differentiate fenv xenv x u) ,v)
    238                                   (* ,u ,(differentiate fenv xenv x v))))
    239 
    240                 (('/ u v)     `(/ (- (* ,(differentiate fenv xenv x u) ,v)
    241                                      (* ,u ,(differentiate fenv xenv x v)))
     234                (('neg u)  `(neg ,(differentiate fenv x u)))
     235
     236                (('+ u v)  `(+ ,(differentiate fenv x u) ,(differentiate fenv x v)))
     237                (('- u v)  `(- ,(differentiate fenv x u) ,(differentiate fenv x v)))
     238
     239                (('* (and u (? number?)) v)    `(* ,u ,(differentiate fenv x v)))
     240                (('* v (and u (? number?)))    `(* ,u ,(differentiate fenv x v)))
     241
     242                (('* u v)     `(+ (* ,(differentiate fenv x u) ,v)
     243                                  (* ,u ,(differentiate fenv x v))))
     244
     245                (('/ u v)     `(/ (- (* ,(differentiate fenv x u) ,v)
     246                                     (* ,u ,(differentiate fenv x v)))
    242247                                  (pow ,v 2.0)))
    243248
    244                 (('cube u)     (differentiate fenv xenv x `(pow ,u 3.0)))
    245 
    246                 (('pow u n)    (chain fenv xenv x u `(* ,n (pow ,u (- ,n 1.0)))))
     249                (('cube u)     (differentiate fenv x `(pow ,u 3.0)))
     250
     251                (('pow u n)    (chain fenv x u `(* ,n (pow ,u (- ,n 1.0)))))
    247252               
    248                 (('sqrt u)     (chain fenv xenv x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
    249 
    250                 (('exp u)      (chain fenv xenv x u `(exp ,u)))
    251 
    252                 (('log u)      (chain fenv xenv x u `(/ 1.0 ,u)))
    253 
    254                 (('log10 u)    (chain fenv xenv x u `(* ,LOG10E (/ ,(differentiate fenv xenv x u) ,u))))
     253                (('sqrt u)     (chain fenv x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
     254
     255                (('exp u)      (chain fenv x u `(exp ,u)))
     256
     257                (('log u)      (chain fenv x u `(/ 1.0 ,u)))
     258
     259                (('log10 u)    (chain fenv x u `(* ,LOG10E (/ ,(differentiate fenv x u) ,u))))
    255260       
    256                 (('log2 u)     (chain fenv xenv x u `(* ,LOG2E (/ ,(differentiate fenv xenv x u) ,u))))
    257 
    258                 (('log1p u)    (differentiate fenv xenv x `(log (+ 1.0 ,u))))
    259 
    260                 (('ldexp u n)  (differentiate fenv xenv x `(* ,u ,(expt 2 n))))
     261                (('log2 u)     (chain fenv x u `(* ,LOG2E (/ ,(differentiate fenv x u) ,u))))
     262
     263                (('log1p u)    (differentiate fenv x `(log (+ 1.0 ,u))))
     264
     265                (('ldexp u n)  (differentiate fenv x `(* ,u ,(expt 2 n))))
    261266       
    262                 (('sin u)      (chain fenv xenv x u `(cos ,u)))
     267                (('sin u)      (chain fenv x u `(cos ,u)))
    263268               
    264                 (('cos u)      (chain fenv xenv x u `(neg (sin ,u))))
    265 
    266                 (('tan u)      (differentiate fenv xenv x `(* (sin ,u) (/ 1.0 (cos ,u)))))
     269                (('cos u)      (chain fenv x u `(neg (sin ,u))))
     270
     271                (('tan u)      (differentiate fenv x `(* (sin ,u) (/ 1.0 (cos ,u)))))
    267272               
    268                 (('asin u)     (chain fenv xenv x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
     273                (('asin u)     (chain fenv x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
    269274                                             
    270                 (('acos u)     (chain fenv xenv x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
     275                (('acos u)     (chain fenv x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
    271276                                             
    272                 (('atan u)     (chain fenv xenv x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
     277                (('atan u)     (chain fenv x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
    273278                                             
    274                 (('sinh u)     (differentiate fenv xenv x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
    275 
    276                 (('cosh u)     (differentiate fenv xenv x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
     279                (('sinh u)     (differentiate fenv x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
     280
     281                (('cosh u)     (differentiate fenv x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
    277282               
    278                 (('tanh u)     (differentiate fenv xenv x `(/ (sinh ,u) (cosh ,u))))
    279 
    280                 (('let bnds body)  (let ((xenv1 (fold (lambda (vb xenv)
    281                                                         (cons (car vb) xenv)) xenv bnds)))
    282                                      `(let ,(map (match-lambda ((v b) `(,v ,(differentiate fenv xenv x b))))
    283                                                  bnds)
    284                                         ,(differentiate fenv xenv1 x body))))
    285 
    286                 ((op u)        (cond ((lookup-def op fenv) =>
    287                                       (match-lambda ((fx) (chain fenv xenv x u `(,fx ,u)))
    288                                                     (else #f)))
    289                                      (else #f)))
    290 
    291                 ((op . us)     (cond ((lookup-def op fenv) =>
    292                                       (lambda (fs)
    293                                         `(+ . ,(map (lambda (fu u) (chain fenv xenv x u `(,fu ,u)))
    294                                                     fs us))))
    295                                      (else #f)))
     283                (('tanh u)     (differentiate fenv x `(/ (sinh ,u) (cosh ,u))))
     284
     285                (('let bnds body)  (let ((body1 (subst-convert body bnds)))
     286                                     (differentiate fenv x body1)))
     287
     288                ((op . us)     (let ((fv (enum-freevars t '() '())))
     289                                 (if (member x fv)
     290                                     (cond ((lookup-def op fenv) =>
     291                                            (lambda (fs)
     292                                              (cond ((and (pair? fs) (pair? us))
     293                                                     `(+ . ,(map (lambda (fu u) (chain fenv x u `(,fu ,u)))
     294                                                                 fs us)))
     295                                                    (else (chain fenv x us `(,fs ,us))))))
     296                                           (else #f))
     297                                     0.0)))
    296298
    297299                (else          #f)))))
    298300
    299 (define (chain fenv xenv x t u)
     301(define (chain fenv x t u)
    300302  (if (symbol? t) u
    301       `(* ,(differentiate fenv xenv x t) ,u)))
     303      `(* ,(differentiate fenv x t) ,u)))
    302304
    303305
    304306(define (simplify t)
    305307  (match t
     308         (('neg 0.0)   0.0)
     309
     310         (('+ 0.0 0.0)  0.0)
    306311         (('+ 0.0 t1)  t1)
    307312         (('+ t1 0.0)  t1)
     
    309314         (('+ (and t1 (? number?)) (and t2 (? number?))) (+ t1 t2))
    310315         
    311          
     316         (('- 0.0 0.0)  0.0)
    312317         (('- 0.0 t1)  `(neg ,t1))
    313318         (('- t1 0.0)  t1)
     
    315320         (('- (and t1 (? number?)) (and t2 (? number?))) (- t1 t2))
    316321         
    317          
     322         (('* 0.0 0.0)  0.0)
    318323         (('* 0.0 t1)  0.0)
    319324         (('* t1 0.0)  0.0)
     
    322327         (('* ('neg t1) ('neg t2))  `(* ,t1 ,t2))
    323328         (('* (and t1 (? number?)) (and t2 (? number?))) (* t1 t2))
     329
     330         (('/ 0.0 t1)  0.0)
    324331         
    325332         (('pow t1 0.0)  1.0)
Note: See TracChangeset for help on using the changeset viewer.