Changeset 12167 in project


Ignore:
Timestamp:
10/15/08 06:28:44 (12 years ago)
Author:
Ivan Raikov
Message:

Incorporating initial linear equations.

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

Legend:

Unmodified
Added
Removed
  • release/3/nemo/trunk/examples/AKP06/PotIhCa.scm

    r12129 r12167  
    316316                                   (b03 = (* 3.0 Narsg_beta (exp (/ v Narsg_x2))))
    317317                                   (b04 = (* 4.0 Narsg_beta (exp (/ v Narsg_x2))))
    318                                    (b0O = (* Narsg_delta * (exp (/ v Narsg_x4))))
    319                                    (bip = (* Narsg_zeta * (exp (/ v Narsg_x6))))
     318                                   (b0O = (* Narsg_delta (exp (/ v Narsg_x4))))
     319                                   (bip = (* Narsg_zeta (exp (/ v Narsg_x6))))
    320320
    321321                                   (b11 = (* Narsg_beta Narsg_btfac (exp (/ v Narsg_x2))))
     
    354354
    355355                                     (initial-equilibrium
    356                                         (0 = (+ (* I1 bi1) (* C2 b01) (neg (* C1 (+ fi1 + f01)) )))
     356                                        (0 = (+ (* I1 bi1) (* C2 b01) (neg (* C1 (+ fi1 f01)) )))
    357357                                        (0 = (+ (* C1 f01) (* I2 bi2) (* C3 b02) (neg (* C2 (+ b01 fi2 f02)) )))
    358358                                        (0 = (+ (* C2 f02) (* I3 bi3) (* C4 b03) (neg (* C3 (+ b02 fi3 f03)) )))
     
    364364                                        (0 = (+ (* I1 f11) (* C2 fi2) (* I3 b12) (neg (* I2 (+ b11 bi2 f12)) )))
    365365                                        (0 = (+ (* I2 f12) (* C3 fi3) (* I4 bi3) (neg (* I3 (+ b12 bi3 f13)) )))
    366                                         (0 = (+ (* I3 * f13) (* C4 fi4) (* I5 b14) (neg (* I4 (+ b13 bi4 f14)) )))
     366                                        (0 = (+ (* I3 f13) (* C4 fi4) (* I5 b14) (neg (* I4 (+ b13 bi4 f14)) )))
    367367                                        (0 = (+ (* I4 f14) (* C5 fi5) (* I6 b1n) (neg (* I5 (+ b14 bi5 f1n)) )))
    368368                                        (1 = (+ C1 C2 C3 C4 C5 O B I1 I2 I3 I4 I5 I6 )))
  • release/3/nemo/trunk/nemo-core.scm

    r12129 r12167  
    3939(declare (export make-nemo-core nemo:error nemo:warning
    4040                 nemo:env-copy nemo-intern nemo:quantity?
     41                 nemo:rhs? nemo:lineq?
    4142                 eval-nemo-system-decls
    4243                 TSCOMP ASGN CONST PRIM))
     
    8081(define (optional pred?) (lambda (x) (or (not x) (pred? x))))
    8182
    82 (define (rhs? x)  (or (symbol? x) (number? x) (and (list? x) (every rhs? x))))
     83(define (rhs? x) 
     84  (or (symbol? x) (number? x)
     85      (match x (((? symbol?) . rest)  (every rhs? rest)) (else #f))))
    8386
    8487(define (lineq? x)  (match x (((? integer?) '= (? rhs?))  #t) (else #f)))
     88
     89(define nemo:rhs?    rhs?)
     90(define nemo:lineq?  lineq?)
    8591
    8692(define (transition? x)
     
    96102  (ASGN       (name symbol?) (value number?) (rhs rhs?) )
    97103  (CONST      (name symbol?) (value number?))
    98   (TSCOMP     (name symbol?) (initial (optional rhs?)) (initial-eq (optional lineq?))
     104  (TSCOMP     (name symbol?) (initial (lambda (x) (or (rhs? x) (and (list? x) (every lineq? x)))))
    99105              (open (lambda (x) (or (symbol? x) (and (list? x) (every symbol? x) ))))
    100106              (transitions (lambda (x) (and (list? x) (every transition? x))))
     
    897903                            (('state-complex (id . alst) )
    898904                             (cond ((and (symbol? id) (list? alst))
    899                                     (let ((initial (eval-const (lookup-def 'initial alst)))
     905                                    (let ((initial    (lookup-def 'initial alst))
     906                                          (initial-eq (alist-ref 'initial-equilibrium alst))
    900907                                          (power   (eval-const (lookup-def 'power alst))))
    901                                       (apply ((nemo-core 'env-extend!) sys)
    902                                              (cons* id '(tscomp) initial `(power ,power) alst)))
    903                                     (cons id qs))
     908                                      (if (not (or initial initial-eq))
     909                                          (nemo:error 'eval-nemo-system-decls
     910                                                      "state complex declarations require initial value or "
     911                                                      "initial equilibrium equations"))
     912                                      (if (and initial-eq
     913                                               (or (not (list? initial-eq)) (not (every lineq? initial-eq))))
     914                                          (nemo:error 'eval-nemo-system-decls
     915                                                      "initial equilibrium field in state complex declarations "
     916                                                      "must be a list of linear equations"))
     917                                      (let ((initialv (and initial (eval-const initial))))
     918                                        (apply ((nemo-core 'env-extend!) sys)
     919                                               (cons* id '(tscomp) (or initialv initial-eq) `(power ,power) alst))
     920                                        (cons id qs))))
    904921                                   (else (nemo:error 'eval-nemo-system-decls
    905922                                                        "state complex declarations must be of the form: "
  • release/3/nemo/trunk/nemo-nmodl.scm

    r12129 r12167  
    11
    2 ;; TODO: * include option for generating kinetic eqs
    3 ;;       * check that open states are valid
     2;; TODO: * check that open states are valid
    43;;       
    54;;
     
    237236           (expr4 (name-normalize expr3)))
    238237      expr4)))
     238
    239239
    240240(define (format-expr/NMODL indent expr . rest) 
     
    314314  (let-optionals rest ((rv #f) (width 72))
    315315    (sdoc->string (doc:format width (format-expr/NMODL 2 x rv)))))
     316 
     317
     318(define (format-lineq/NMODL indent expr . rest) 
     319  (let-optionals rest ((rv #f))
     320   (let ((indent+ (+ 2 indent)))
     321    (match expr
     322       (('let bindings body)
     323        (letblk/NMODL
     324         (fold-right
     325           (lambda (x ax)
     326             (letblk/NMODL
     327              (match (second x)
     328                     (('if c t e)
     329                      (ifthen/NMODL
     330                       (group/NMODL (format-lineq/NMODL indent c))
     331                       (block/NMODL (format-lineq/NMODL indent t (first x)))
     332                       (block/NMODL (format-lineq/NMODL indent e (first x)))))
     333                     (else
     334                      (format-op/NMODL indent+ " = "
     335                                       (list (format-lineq/NMODL indent (first x) )
     336                                             (format-lineq/NMODL indent (second x))))))
     337              ax))
     338           (doc:empty) bindings)
     339         (let ((body1 (doc:nest indent (format-lineq/NMODL indent body))))
     340           (if rv  (format-op/NMODL indent " = " (list (format-lineq/NMODL indent+ rv ) body1))
     341               body1))))
     342       
     343       (('if . rest) (error 'format-lineq/NMODL "invalid if statement " expr))
     344
     345       ((op . rest) 
     346        (let ((op (case op ((pow)  '^) ((abs) 'fabs) (else op))))
     347          (let ((fe
     348                (if (member op nmodl-ops)
     349                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
     350                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
     351                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
     352                      (case op
     353                        ((/) 
     354                         (format-op/NMODL indent op
     355                                          (map (lambda (x)
     356                                                 (let ((fx (format-lineq/NMODL indent+ x)))
     357                                                   (if (or (symbol? x) (number? x)) fx
     358                                                       (if (or mul? plmin?) (group/NMODL fx) fx)))) rest)))
     359                        ((*) 
     360                         (format-op/NMODL indent op
     361                                          (map (lambda (x)
     362                                                 (let ((fx (format-lineq/NMODL indent+ x)))
     363                                                   (if (or (symbol? x) (number? x)) fx
     364                                                       (if plmin? (group/NMODL fx) fx)))) rest)))
     365                       
     366                        ((^) 
     367                         (format-op/NMODL indent op
     368                                          (map (lambda (x)
     369                                                 (let ((fx (format-lineq/NMODL indent+ x)))
     370                                                   (if (or (symbol? x)  (number? x)) fx
     371                                                       (if (or mdiv? plmin?) (group/NMODL fx) fx)))) rest)))
     372                       
     373                        (else
     374                         (format-op/NMODL indent op
     375                                          (map (lambda (x)
     376                                                 (let ((fx (format-lineq/NMODL indent+ x))) fx)) rest)))))
     377                   
     378                    (case op
     379                      ((neg) (format-op/NMODL indent '* (map (lambda (x) (format-lineq/NMODL indent+ x))
     380                                                             (cons "(-1)" rest))))
     381                      (else  (format-fncall/NMODL indent op (map (lambda (x) (format-lineq/NMODL indent+ x))
     382                                                                 rest)))))))
     383
     384           (if rv (format-op/NMODL indent " = " (list (format-lineq/NMODL indent+ rv ) fe)) fe))))
     385     
     386      (else  (let ((fe (doc:text (->string expr))))
     387               (if rv
     388                   (format-op/NMODL indent " = " (list (format-lineq/NMODL indent+ rv ) fe))
     389                   fe)))))))
     390               
     391
     392         
     393(define (lineq->string/NMODL x val . rest)
     394  (let-optionals rest ((width 72))
     395    (s+ "~ " (sdoc->string (doc:format width (format-lineq/NMODL 2 x #f)))
     396        " = " (number->string val))))
    316397 
    317398
     
    467548    (list (nmodl-name n) init1)))
    468549
     550
     551(define (state-init-eq n transitions init)
     552  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
     553         (state-list     (let loop ((lst (list)) (tlst transitions))
     554                           (if (null? tlst)  (delete-duplicates lst eq?)
     555                               (match (car tlst)
     556                                      (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
     557                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     558                                      (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
     559                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     560                                      (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
     561                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     562                                      (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
     563                                       (loop (cons* s0 s1 lst) (cdr tlst)))
     564                                      (else
     565                                       (nemo:error 'nemo:state-init-eq ": invalid transition equation "
     566                                                   (car tlst) " in state complex " n))
     567                                      (else (loop lst (cdr tlst)))))))
     568         (state-subs     (fold (lambda (s ax) (subst-extend s (nmodl-state-name n s) ax)) subst-empty state-list))
     569         (init1          (map (lambda (lineq) (match lineq ((i '= . expr) `(,i = . ,(subst-convert expr state-subs)))))
     570                              init)))
     571    (list (nmodl-name n) init1)))
     572
    469573(define (asgn-eq n rhs)
    470574  (let* ((fbody   (rhsexpr rhs))
     
    547651                           (if (nemo:quantity? en)
    548652                               (cases nemo:quantity en
    549                                       (TSCOMP  (name initial open transitions power)
    550                                                (cons* (state-init name initial)
    551                                                       (state-init (nmodl-state-name name open) name) ax))
     653                                      (TSCOMP  (name initial open transitions power)
     654                                               (if (nemo:rhs? initial)
     655                                                   (cons* (state-init name initial)
     656                                                          (state-init (nmodl-state-name name open) name) ax)
     657                                                   ax))
    552658                                      (else  ax))
    553659                               ax))))
    554660            ax lst))
    555661   (list) poset))
     662
     663
     664(define (poset->state-init-eq-defs poset sys)
     665  (fold-right
     666   (lambda (lst ax)
     667     (fold  (lambda (x ax)
     668              (match-let (((i . n)  x))
     669                         (let ((en (environment-ref sys n)))
     670                           (if (nemo:quantity? en)
     671                               (cases nemo:quantity en
     672                                      (TSCOMP (name initial open transitions power)
     673                                              (if (and (list? initial) (every nemo:lineq? initial))
     674                                                  (cons (state-init-eq name transitions initial) ax)
     675                                                  ax))
     676                                      (else  ax))
     677                               ax))))
     678            ax lst))
     679   (list) poset))
     680
    556681
    557682(define (find-locals defs)
     
    816941                 
    817942             
    818              (pp indent ,nl (INITIAL "{"))
    819              (let* ((init-defs  (poset->state-init-defs poset sys))
    820                     (locals     (concatenate (find-locals (map second init-defs)))) )
     943             (let* ((init-defs     (poset->state-init-defs poset sys))
     944                    (init-eq-defs  (poset->state-init-eq-defs poset sys))
     945                    (locals        (concatenate (find-locals (map second init-defs)))) )
     946               (pp indent ,nl (INITIAL "{"))
    821947               (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
    822948               (if (not (null? asgns))  (pp indent+ (rates ())))
     
    824950                           (let ((n (first def)) (b (second def)))
    825951                             (pp indent+ ,(expr->string/NMODL b n)))) init-defs)
    826                (for-each (lambda (x)  (pp indent+ (,(third x) = ,(fourth x))))
    827                          perm-ions))
    828              (pp indent "}")
    829              
    830              )))
     952               (for-each (lambda (x)  (pp indent+ (,(third x) = ,(fourth x))))  perm-ions)
     953               (if (not (null? init-eq-defs)) (pp indent+ (SOLVE initial_equilibrium)))
     954               (pp indent "}")
     955               (if (not (null? init-eq-defs))
     956                   (begin
     957                     (pp indent ,nl (LINEAR initial_equilibrium "{"))
     958                     (for-each
     959                      (lambda (x)
     960                        (let ((lineqs  (second x)))
     961                          (for-each (lambda (eq)
     962                                      (let ((val  (first eq))
     963                                            (expr (third eq)))
     964                                        (pp indent+ ,(lineq->string/NMODL expr val))))
     965                                    lineqs)))
     966                      init-eq-defs)
     967                     (pp indent "}")))
     968               ))))
    831969        )))))
Note: See TracChangeset for help on using the changeset viewer.