Changeset 30616 in project


Ignore:
Timestamp:
03/30/14 13:56:10 (7 years ago)
Author:
Ivan Raikov
Message:

nemo: moving symbolic differentiation code away from utils module

File:
1 edited

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/nemo-utils.scm

    r29841 r30616  
    4343  poset->init-defs*
    4444
    45   differentiate simplify distribute
    4645  make-output-fname
    4746  s+ sw+ slp nl spaces ppf
     
    828827
    829828
    830 (define (differentiate fenv x t)
    831   (define subst-convert
    832     (subst-driver
    833      (lambda (x) (and (symbol? x) x))
    834      nemo:binding? identity nemo:bind nemo:subst-term))
    835  
    836   (cond ((number? t)  0.0)
    837         ((symbol? t)  (cond ((string=? (->string x) (->string t)) 1.0)
    838                             (else 0.0)))
    839         (else (match t
    840                 (('neg u)  `(neg ,(differentiate fenv x u)))
    841 
    842                 (('+ u v)  `(+ ,(differentiate fenv x u) 
    843                                ,(differentiate fenv x v)))
    844 
    845                 (('- u v)  `(- ,(differentiate fenv x u)
    846                                ,(differentiate fenv x v)))
    847 
    848                 (('* (and u (? number?)) v)    `(* ,u ,(differentiate fenv x v)))
    849                 (('* v (and u (? number?)))    `(* ,u ,(differentiate fenv x v)))
    850 
    851                 (('* u v)     `(+ (* ,(differentiate fenv x u) ,v)
    852                                   (* ,u ,(differentiate fenv x v))))
    853 
    854                 (('/ u v)     `(/ (- (* ,(differentiate fenv x u) ,v)
    855                                      (* ,u ,(differentiate fenv x v)))
    856                                   (pow ,v 2.0)))
    857 
    858                 (('cube u)     (differentiate fenv x `(pow ,u 3.0)))
    859 
    860                 (('pow u n)    (chain fenv x u `(* ,n (pow ,u (- ,n 1.0)))))
    861                
    862                 (('sqrt u)     (chain fenv x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
    863 
    864                 (('exp u)      (chain fenv x u `(exp ,u)))
    865 
    866                 (('log u)      (chain fenv x u `(/ 1.0 ,u)))
    867 
    868                 (('log10 u)    (chain fenv x u `(* ,LOG10E (/ ,(differentiate fenv x u) ,u))))
    869        
    870                 (('log2 u)     (chain fenv x u `(* ,LOG2E (/ ,(differentiate fenv x u) ,u))))
    871 
    872                 (('log1p u)    (differentiate fenv x `(log (+ 1.0 ,u))))
    873 
    874                 (('ldexp u n)  (differentiate fenv x `(* ,u ,(expt 2 n))))
    875        
    876                 (('sin u)      (chain fenv x u `(cos ,u)))
    877                
    878                 (('cos u)      (chain fenv x u `(neg (sin ,u))))
    879 
    880                 (('tan u)      (differentiate fenv x `(* (sin ,u) (/ 1.0 (cos ,u)))))
    881                
    882                 (('asin u)     (chain fenv x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
    883                                              
    884                 (('acos u)     (chain fenv x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
    885                                              
    886                 (('atan u)     (chain fenv x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
    887                                              
    888                 (('sinh u)     (differentiate fenv x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
    889 
    890                 (('cosh u)     (differentiate fenv x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
    891                
    892                 (('tanh u)     (differentiate fenv x `(/ (sinh ,u) (cosh ,u))))
    893 
    894                 (('let bnds body)  (let* ((body1 (subst-convert body bnds)))
    895                                      (differentiate fenv x body1)))
    896 
    897                 ((op . us)     (let ((fv (enum-freevars t '() '())))
    898                                  (if (member x fv)
    899                                      (cond ((lookup-def op fenv) =>
    900                                             (lambda (fs)
    901                                               (cond ((and (pair? fs) (pair? us))
    902                                                      `(+ . ,(map (lambda (fu u) (chain fenv x u `(,fu ,u)))
    903                                                                  fs us)))
    904                                                     (else (chain fenv x us `(,fs ,us))))))
    905                                            (else #f))
    906                                      0.0)))
    907 
    908                 (else          #f)))))
    909 
    910 (define (chain fenv x t u)
    911   (if (symbol? t) u
    912       `(* ,(differentiate fenv x t) ,u)))
    913 
    914 
    915 (define (simplify t)
    916   (match t
    917          (('neg 0.0)   0.0)
    918 
    919          (('+ 0.0 0.0)  0.0)
    920          (('+ 0.0 t1)  t1)
    921          (('+ t1 0.0)  t1)
    922          (('+ t1 ('neg t2))  `(- ,t1 ,t2))
    923          (('+ (and t1 (? number?)) (and t2 (? number?))) (+ t1 t2))
    924          
    925          (('- 0.0 0.0)  0.0)
    926          (('- 0.0 t1)  `(neg ,t1))
    927          (('- t1 0.0)  t1)
    928          (('neg ('neg t1))  t1)
    929          (('- (and t1 (? number?)) (and t2 (? number?))) (- t1 t2))
    930          
    931          (('* 0.0 0.0)  0.0)
    932          (('* 0.0 t1)  0.0)
    933          (('* t1 0.0)  0.0)
    934          (('* 1.0 t1)  t1)
    935          (('* t1 1.0)  t1)
    936          (('* ('neg t1) ('neg t2))  `(* ,t1 ,t2))
    937          (('* (and t1 (? number?)) (and t2 (? number?))) (* t1 t2))
    938 
    939          (('/ 0.0 t1)  0.0)
    940          
    941          (('pow t1 0.0)  1.0)
    942          (('pow t1 1.0)  t1)
    943          (('pow (and t1 (? number?)) (and t2 (? number?))) (expt t1 t2))
    944          
    945          (('let () body) (simplify body))
    946          
    947          (('let bnds body)
    948           `(let ,(map (match-lambda ((v b) `(v ,(simplify b)))
    949                                     (else #f)) bnds)
    950              ,(simplify body)))
    951          
    952          ((op . ts) 
    953           `(,op . ,(map simplify ts)))
    954          
    955          (else t)))
    956 
    957 (define (distribute t)
    958   (match t
    959 
    960          (((and (or '+ '- '* '/) op) x y)
    961           `(,op ,(distribute x) ,(distribute y)))
    962 
    963          (((and (or '+ '- '* '/) op) x y z)
    964           `(,op ,(distribute x) (,op ,(distribute y) ,(distribute z))))
    965 
    966          (((and (or '+ '- '* '/) op) . lst)
    967           (let* ((n   (length lst))
    968                  (n/2 (inexact->exact (round (/ n 2)))))
    969             `(,op ,(distribute `(,op . ,(take lst n/2)))
    970                   ,(distribute `(,op . ,(drop lst n/2 ))))))
    971 
    972          (('let bnds body)
    973           `(let ,(map (match-lambda ((v b) `(,v ,(distribute b)))
    974                                     (else #f)) bnds)
    975              ,(distribute body)))
    976          
    977          ((op . ts) 
    978           `(,op . ,(map distribute ts)))
    979 
    980          (else t)))
    981 
    982 
    983829)
Note: See TracChangeset for help on using the changeset viewer.