Changeset 15459 in project


Ignore:
Timestamp:
08/14/09 09:39:06 (10 years ago)
Author:
Ivan Raikov
Message:

some fixes and extensions to the differentiation and simplification
routines in nemo-utils

File:
1 edited

Legend:

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

    r15425 r15459  
    221221(define LOG2E  1.44269504088896)
    222222
    223 (define (differentiate x t)
     223(define (differentiate fenv x t)
    224224  (cond ((number? t)  0.0)
    225225        ((symbol? t)  (if (equal? x t) 1.0 0.0))
    226226        (else (match t
    227                 (('neg u)  `(neg ,(differentiate x u)))
    228 
    229                 (('+ u v)  `(+ ,(differentiate x u) ,(differentiate x v)))
    230                 (('- u v)  `(- ,(differentiate x u) ,(differentiate x v)))
    231 
    232                 (('* (and u (? number?)) v)    `(* ,u ,(differentiate x v)))
    233                 (('* v (and u (? number?)))    `(* ,u ,(differentiate x v)))
    234 
    235                 (('* u v)     `(+ (* ,(differentiate x u) ,v)
    236                                   (* ,u ,(differentiate x v))))
    237 
    238                 (('/ u v)     `(/ (- (* ,(differentiate x u) ,v)
    239                                      (* ,u ,(differentiate x v)))
     227                (('neg u)  `(neg ,(differentiate fenv x u)))
     228
     229                (('+ u v)  `(+ ,(differentiate fenv x u) ,(differentiate fenv x v)))
     230                (('- u v)  `(- ,(differentiate fenv x u) ,(differentiate fenv x v)))
     231
     232                (('* (and u (? number?)) v)    `(* ,u ,(differentiate fenv x v)))
     233                (('* v (and u (? number?)))    `(* ,u ,(differentiate fenv x v)))
     234
     235                (('* u v)     `(+ (* ,(differentiate fenv x u) ,v)
     236                                  (* ,u ,(differentiate fenv x v))))
     237
     238                (('/ u v)     `(/ (- (* ,(differentiate fenv x u) ,v)
     239                                     (* ,u ,(differentiate fenv x v)))
    240240                                  (pow ,v 2.0)))
    241241
    242                 (('cube u)     (differentiate x `(pow ,u 3.0)))
    243 
    244                 (('pow u n)    (chain x u `(* ,n (pow ,u (- ,n 1.0)))))
     242                (('cube u)     (differentiate fenv x `(pow ,u 3.0)))
     243
     244                (('pow u n)    (chain fenv x u `(* ,n (pow ,u (- ,n 1.0)))))
    245245               
    246                 (('sqrt u)     (chain x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
    247 
    248                 (('exp u)      (chain x u `(exp ,u)))
    249 
    250                 (('log a u)    (chain x u `(/ 1.0 ,u)))
    251 
    252                 (('log10 u)    (chain x u `(* ,LOG10E (/ ,(differentiate x u) ,u))))
     246                (('sqrt u)     (chain fenv x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
     247
     248                (('exp u)      (chain fenv x u `(exp ,u)))
     249
     250                (('log u)      (chain fenv x u `(/ 1.0 ,u)))
     251
     252                (('log10 u)    (chain fenv x u `(* ,LOG10E (/ ,(differentiate fenv x u) ,u))))
    253253       
    254                 (('log2 u)     (chain x u `(* ,LOG2E (/ ,(differentiate x u) ,u))))
    255 
    256                 (('log1p u)    (differentiate x `(log (+ 1.0 ,u))))
    257 
    258                 (('ldexp u n)  (differentiate x `(* ,u ,(expt 2 n))))
     254                (('log2 u)     (chain fenv x u `(* ,LOG2E (/ ,(differentiate fenv x u) ,u))))
     255
     256                (('log1p u)    (differentiate fenv x `(log (+ 1.0 ,u))))
     257
     258                (('ldexp u n)  (differentiate fenv x `(* ,u ,(expt 2 n))))
    259259       
    260                 (('sin u)      (chain x u `(cos ,u)))
     260                (('sin u)      (chain fenv x u `(cos ,u)))
    261261               
    262                 (('cos u)      (chain x u `(neg (sin ,u))))
    263 
    264                 (('tan u)      (differentiate x `(* (sin ,u) (/ 1.0 (cos ,u)))))
     262                (('cos u)      (chain fenv x u `(neg (sin ,u))))
     263
     264                (('tan u)      (differentiate fenv x `(* (sin ,u) (/ 1.0 (cos ,u)))))
    265265               
    266                 (('asin u)     (chain x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
     266                (('asin u)     (chain fenv x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
    267267                                             
    268                 (('acos u)     (chain x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
     268                (('acos u)     (chain fenv x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
    269269                                             
    270                 (('atan u)     (chain x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
     270                (('atan u)     (chain fenv x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
    271271                                             
    272                 (('sinh u)     (differentiate x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
    273 
    274                 (('cosh u)     (differentiate x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
     272                (('sinh u)     (differentiate fenv x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
     273
     274                (('cosh u)     (differentiate fenv x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
    275275               
    276                 (('tanh u)     (differentiate x `(/ (sinh ,u) (cosh ,u))))
     276                (('tanh u)     (differentiate fenv x `(/ (sinh ,u) (cosh ,u))))
     277
     278                (('let bnds body)  `(let ,(map (match-lambda ((v b) `(,v ,(differentiate fenv x b))))
     279                                               bnds)
     280                                      ,(differentiate fenv x body)))
     281                                                             
     282
     283                ((op u)        (cond ((lookup-def op fenv) =>
     284                                      (match-lambda ((fx) (chain fenv x u `(,fx ,u)))
     285                                                    (else #f)))
     286                                     (else #f)))
     287
     288                ((op . us)     (cond ((lookup-def op fenv) =>
     289                                      (lambda (fs)
     290                                        `(+ . ,(map (lambda (fu u) (chain fenv x u `(,fu ,u)))
     291                                                    fs us))))
     292                                     (else #f)))
    277293
    278294                (else          #f)))))
    279295
    280 (define (chain x t u)
     296(define (chain fenv x t u)
    281297  (if (symbol? t) u
    282       `(* ,(differentiate x t) ,u)))
    283        
     298      `(* ,(differentiate fenv x t) ,u)))
     299
     300
    284301(define (simplify t)
    285302  (match t
     
    287304         (('+ t1 0.0)  t1)
    288305         (('+ t1 ('neg t2))  `(- ,t1 ,t2))
     306         (('+ (and t1 (? number?)) (and t2 (? number?))) (+ t1 t2))
     307         
    289308         
    290309         (('- 0.0 t1)  `(neg ,t1))
    291310         (('- t1 0.0)  t1)
    292311         (('neg ('neg t1))  t1)
    293 
     312         (('- (and t1 (? number?)) (and t2 (? number?))) (- t1 t2))
     313         
     314         
    294315         (('* 0.0 t1)  0.0)
    295316         (('* t1 0.0)  0.0)
    296317         (('* 1.0 t1)  t1)
    297          (('* t1 1.0)  t1)
    298          (('* ('neg t1) ('neg t2))  `(* ,t1 ,t2))
    299 
     318         (('* t1 1.0)  t1)
     319         (('* ('neg t1) ('neg t2))  `(* ,t1 ,t2))
     320         (('* (and t1 (? number?)) (and t2 (? number?))) (* t1 t2))
     321         
     322         (('pow t1 0.0)  1.0)
     323         (('pow t1 1.0)  t1)
     324         (('pow (and t1 (? number?)) (and t2 (? number?))) (expt t1 t2))
     325         
     326         (('let bnds body)
     327          `(let ,(map (match-lambda ((v b) `(v ,(recur b)))
     328                                    (else #f)) bnds)
     329             ,(recur body)))
     330         
     331         ((op . ts) 
     332          `(,op . ,(map simplify ts)))
     333         
    300334         (else t)))
    301335
Note: See TracChangeset for help on using the changeset viewer.