Changeset 12466 in project


Ignore:
Timestamp:
11/11/08 04:17:50 (13 years ago)
Author:
Ivan Raikov
Message:

Added support for the expeuler method.

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

Legend:

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

    r12369 r12466  
    133133(define (format-op/NMODL indent op args)
    134134  (let ((op1 (doc:text (->string op))))
    135     (if (null? args) op1
    136         (match args
    137                ((x)           (doc:connect op1 x))
    138                ((x y)         (binop/NMODL x op1 y))
    139                ((x y z)       (binop/NMODL x op1 (binop/NMODL y op1 z)))
    140                (lst           (let* ((n   (length lst))
    141                                      (n/2 (inexact->exact (round (/ n 2)))))
    142                                 (binop/NMODL (format-op/NMODL indent op (take lst n/2 )) op1
    143                                          (format-op/NMODL indent op (drop lst n/2 )))))))))
     135    (let ((res
     136           (if (null? args) op1
     137               (match args
     138                      ((x)           (doc:connect op1 x))
     139                      ((x y)         (binop/NMODL x op1 y))
     140                      ((x y z)       (binop/NMODL x op1 (binop/NMODL y op1 z)))
     141                      (lst           (let* ((n   (length lst))
     142                                            (n/2 (inexact->exact (round (/ n 2)))))
     143                                       (binop/NMODL (format-op/NMODL indent op (take lst n/2 )) op1
     144                                                    (format-op/NMODL indent op (drop lst n/2 )))))))))
     145      res)))
    144146
    145147
     
    273275         (fold-right
    274276           (lambda (x ax)
    275              (letblk/NMODL
    276               (match (second x)
    277                      (('if c t e)
    278                       (ifthen/NMODL
    279                        (group/NMODL (format-expr/NMODL indent c))
    280                        (block/NMODL (format-expr/NMODL indent t (first x)))
    281                        (block/NMODL (format-expr/NMODL indent e (first x)))))
    282                      (else
    283                       (format-op/NMODL indent+ " = "
    284                                        (list (format-expr/NMODL indent (first x) )
    285                                              (format-expr/NMODL indent (second x))))))
    286               ax))
     277             (let ((res
     278                    (letblk/NMODL
     279                     (match (second x)
     280                            (('if c t e)
     281                             (ifthen/NMODL
     282                              (group/NMODL (format-expr/NMODL indent c))
     283                              (block/NMODL (format-expr/NMODL indent t (first x)))
     284                              (block/NMODL (format-expr/NMODL indent e (first x)))))
     285                            (else
     286                             (format-op/NMODL indent+ " = "
     287                                              (list (format-expr/NMODL indent (first x) )
     288                                                    (format-expr/NMODL indent (second x))))))
     289                     ax)))
     290               res
     291               ))
    287292           (doc:empty) bindings)
    288293         (let ((body1 (doc:nest indent (format-expr/NMODL indent body))))
     
    454459  (pp indent (,n)))
    455460
    456 
    457 (define (state-eqs n initial open transitions power)
     461(define (expeuler dt name rhs)
     462  (define (isname? x) (equal? x name))
     463  (let ((res
     464         (match rhs
     465                ((or ('- A ('* B (and x (? isname?))))
     466                     ('+ ('neg ('* B (and x (? isname?)))) A))
     467                 (let ((xexp (string->symbol (s+ x 'exp))))
     468                   `(let ((,xexp (exp (* (neg ,B) ,dt))))
     469                      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B))))))
     470
     471                ((or ('- A ('* (and x (? isname?)) . B))
     472                     ('+ ('neg ('* (and x (? isname?)) . B)) A))
     473                 (let ((xexp (string->symbol (s+ x 'exp)))
     474                       (B1   (if (null? (cdr B)) (car B) `(* ,@B))))
     475                   `(let ((,xexp (exp (* (neg ,B1) ,dt))))
     476                      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B1))))))
     477               
     478                (('+ ('neg ('* (and x1 (? isname?)) Alpha))
     479                     ('* ('- 1 (and x2 (? isname?))) Beta))
     480                 (let ((A  Alpha)
     481                       (B  `(+ ,Alpha ,Beta)))
     482                   (let ((xexp (string->symbol (s+ x1 'exp))))
     483                     `(let ((,xexp (exp (* (neg ,B) ,dt))))
     484                        (+ (* ,x1 ,xexp) (* (- 1 ,xexp) (/ ,A ,B)))))))
     485               
     486                (('let bnds body)
     487                 `(let ,bnds ,(expeuler dt name body)))
     488               
     489                (else (nemo:error 'nemo:expeuler ": unable to rewrite equation " rhs
     490                                  "in exponential Euler form")))))
     491
     492    res))
     493 
     494
     495(define (state-eqs n initial open transitions power method)
    458496  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
    459497         (g          (make-digraph n (string-append (->string n) " transitions graph")))
     
    519557                                                      ((and (null? out) (not (null? in)))
    520558                                                       (sum (map third in)))))
    521                                          (fbody  (rhsexpr rhs1))
    522                                          (fbody1 (canonicalize-expr/NMODL fbody)))
    523                                     (cons (list (s+ name "'") fbody1) ax)))))
     559                                         (fbody0 (rhsexpr rhs1)))
     560                                    (case method
     561                                      ((expeuler)  (cons (list name (canonicalize-expr/NMODL (expeuler 'dt name fbody0)))
     562                                                         ax))
     563                                      (else        (cons (list (s+ name "'") (canonicalize-expr/NMODL fbody0))
     564                                                         ax))
     565                                      )))))
    524566                          (list) nodes)))
    525567        eqs))))
     
    630672
    631673
    632 (define (poset->state-eq-defs poset sys kinetic)
     674(define (poset->state-eq-defs poset sys kinetic method)
    633675  (fold-right
    634676   (lambda (lst ax)
     
    639681                               (cases nemo:quantity en
    640682                                      (TSCOMP  (name initial open transitions conserve power)
    641                                                (append (state-eqs name initial open transitions power) ax))
     683                                               (append (state-eqs name initial open transitions power method) ax))
    642684                                      (else  ax))
    643685                               ax))))
     
    938980                             (pp indent+ ,(expr->string/NMODL depth-val depth-name)))))
    939981                       epools))
     982           (case method  ((expeuler)  (pp indent+ dt)))
    940983           (pp indent "}")
    941984           
     
    10821125             (if (not (null? asgns))     (pp indent+ (rates ())))
    10831126             (if has-ode?
    1084                  (if (not method) (pp indent+ (SOLVE states))
    1085                      (pp indent+ (SOLVE states METHOD ,method))))
     1127                 (case method
     1128                   ((#f expeuler)  (pp indent+ (SOLVE states)))
     1129                   (else           (pp indent+ (SOLVE states METHOD ,method)))))
    10861130             (if has-kinetic?   (pp indent+  (SOLVE kstates METHOD sparse)))
    10871131             (if (not (null? stcomps))   (pp indent+ (stcomps ())))
     
    10911135           
    10921136           (if has-ode?
    1093                (begin
    1094                  (pp indent ,nl (DERIVATIVE states "{"))
    1095                  (let* ((state-eq-defs  (reverse (poset->state-eq-defs poset sys kinetic)))
    1096                         (pool-eq-defs
    1097                          (map (lambda (ep)
    1098                                 (let* ((ep-name     (first ep))
    1099                                        (pool-ion    (assoc ep-name pool-ions))
    1100                                        (i-name      (second pool-ion))
    1101                                        (init-name   (nmodl-name (s+ ep-name '-init)))
    1102                                        (temp-name   (nmodl-name (s+ ep-name '-temp-adj)))
    1103                                        (beta-name   (nmodl-name (s+ ep-name '-beta)))
    1104                                        (depth-name  (nmodl-name (s+ ep-name '-depth)))
    1105                                        (rhs         `(let ((F 96485.0))
    1106                                                        (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name))
    1107                                                           (* ,beta-name ,ep-name .
    1108                                                              ,(if temp-name (list temp-name) (list)))))))
    1109                                   `(,(s+ ep-name "'") ,rhs)))
    1110                               epools))
    1111                         (eq-defs (append pool-eq-defs state-eq-defs))
    1112                         (locals  (find-locals (map second eq-defs))))
    1113                    (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
    1114                    (for-each (lambda (def)
    1115                                (let ((n (first def)) (b (second def)))
    1116                                  (pp indent+ ,(expr->string/NMODL b n)))) eq-defs))
    1117                    
     1137               (let* ((state-eq-defs  (reverse (poset->state-eq-defs poset sys kinetic method)))
     1138                      (pool-eq-defs
     1139                       (map (lambda (ep)
     1140                              (let* ((ep-name     (first ep))
     1141                                     (pool-ion    (assoc ep-name pool-ions))
     1142                                     (i-name      (second pool-ion))
     1143                                     (init-name   (nmodl-name (s+ ep-name '-init)))
     1144                                     (temp-name   (nmodl-name (s+ ep-name '-temp-adj)))
     1145                                     (beta-name   (nmodl-name (s+ ep-name '-beta)))
     1146                                     (depth-name  (nmodl-name (s+ ep-name '-depth)))
     1147                                     (rhs         `(let ((F 96485.0))
     1148                                                     (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name))
     1149                                                        (* ,ep-name ,beta-name .
     1150                                                           ,(if temp-name (list temp-name) (list)))))))
     1151                                (case method
     1152                                  ((expeuler)  `(,ep-name ,(canonicalize-expr/NMODL (expeuler 'dt ep-name rhs))))
     1153                                  (else        `(,(s+ ep-name "'") ,(canonicalize-expr/NMODL rhs))))))
     1154                            epools))
     1155                      (eq-defs (append pool-eq-defs state-eq-defs))
     1156                      (locals  (find-locals (map second eq-defs))))
     1157                 (case method
     1158                   ((expeuler) (pp indent ,nl (PROCEDURE states () "{")))
     1159                   (else       (pp indent ,nl (DERIVATIVE states "{"))))
     1160                 (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
     1161                 (for-each (lambda (def)
     1162                             (let ((n (first def)) (b (second def)))
     1163                               (pp indent+ ,(expr->string/NMODL b n))))
     1164                           eq-defs)
    11181165                 (pp indent "}")))
    1119            
    11201166
    11211167           (if has-kinetic?
  • release/3/nemo/trunk/nemo.scm

    r12376 r12466  
    635635                     ((adams runge euler adeuler heun adrunge gear
    636636                             newton simplex simeq seidel sparse derivimplicit cnexp clsoda
    637                              after_cvode cvode_t cvode_t_v #f) method)
     637                             after_cvode cvode_t cvode_t_v expeuler #f) method)
    638638                     (else (error "unknown nmodl-method " method))))))
    639639           (if sxml-fname (with-output-to-file sxml-fname (lambda () (pretty-print (model->ncml options model)))))
Note: See TracChangeset for help on using the changeset viewer.