Changeset 29931 in project


Ignore:
Timestamp:
10/16/13 11:10:32 (8 years ago)
Author:
Ivan Raikov
Message:

nemo: added support for steady state solutions in the nest backend using kinsol

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

Legend:

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

    r29882 r29931  
    377377            (pp indent ,nl (function y = ,(s+ sysname "_steadystate") (,(slp ", " '(x v params imports))))))
    378378           
    379        
    380379        (for-each
    381380         (lambda (def)
  • release/4/nemo/trunk/nemo-nest.scm

    r29930 r29931  
    298298
    299299
    300 (define (output-prelude sysname steady-state-index-map indent)
     300(define (output-prelude sysname indent)
    301301
    302302  (pp indent ,#<#EOF
     
    322322
    323323
     324
     325(define fsolve-prelude-header
     326#<<EOF
     327  /**
     328   * Exception to be thrown if a KinSol solver does not return KIN_SUCCESS
     329   * @ingroup KernelExceptions
     330   */
     331  class KINSolverFailure: public KernelException
     332  {
     333  public:
     334    /**
     335    * @note model should be passed from get_name() to ensure that
     336    *             names of copied models are reported correctly.
     337    * @param model name of model causing problem
     338    * @param status exit status of the KINSOL solver
     339    */
     340    KINSolverFailure(const std::string& model,
     341                     const int status)
     342      : KernelException("KINSolverFailure"),
     343        model_(model),
     344        status_(status)
     345      {}
     346    ~KINSolverFailure() throw() {}
     347   
     348    std::string message()
     349    {
     350      std::ostringstream msg;
     351      msg << "In model " << model_ << ", the KINSOL solver "
     352          << "returned with exit status " << status_ << ".\n";
     353      return msg.str();
     354    }
     355
     356    private:
     357      const std::string model_;
     358      const int status_;
     359  };
     360
     361EOF
     362)
    324363
    325364
     
    471510)
    472511
     512(define fsolve-prelude
     513#<<EOF
     514
     515  int fsolve (KINSysFn f, N_Vector fval, void *user_data, std::string name)
     516  {
     517      int status, N;
     518      N_Vector u0, sc;
     519      void *kmem;
     520
     521      N = NV_LENGTH_S(fval);
     522
     523      u0 = N_VNew_Serial(N);
     524      N_VConst_Serial(0.0,u0);
     525      N_VConst_Serial(0.0,fval);
     526
     527      sc = N_VNew_Serial(N);
     528      N_VConst_Serial(1.0,sc);
     529
     530      kmem = KINCreate();
     531
     532      status = KINSetUserData (kmem, user_data);
     533      if (check_flag (&status, "KinSetUserData", 1)) throw KINSolverFailure (name, status);
     534
     535      status = KINInit (kmem, f, u0);
     536      if (check_flag (&status, "KinInit", 1)) throw KINSolverFailure (name, status);
     537
     538      status = KINDense (kmem, N);
     539      if (check_flag (&status, "KinDense", 1)) throw KINSolverFailure (name, status);
     540
     541      status = KINSol (kmem, fval, KIN_NONE, sc, sc);
     542      if (check_flag (&status, "KINSol", 1)) throw KINSolverFailure (name, status);
     543     
     544      N_VDestroy_Serial(sc);
     545      N_VDestroy_Serial(u0);
     546
     547      KINFree (&kmem);
     548
     549      return 0;
     550  }
     551
     552EOF
     553)
     554
     555
    473556
    474557(define (output-event-handle
    475          sysname state-index-map steady-state-index-map
     558         sysname state-index-map
    476559         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    477560         reaction-eq-defs i-eqs pool-ions perm-ions
     
    530613))
    531614
     615
    532616(define (output-header
    533617         sysname method state-index-map steady-state-index-map
     
    556640      ((cvode)
    557641       (pp indent
     642           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
     643           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
    558644           ("#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */")
    559            ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
    560645           ("#include <cvode/cvode_diag.h>        /* prototype for CVDiag */")
    561            ("#include <sundials/sundials_types.h> /* definition of type realtype */")
     646           ("#include <kinsol/kinsol.h>           /* prototypes for KINSOL fcts. */")
     647           ("#include <kinsol/kinsol_dense.h>     /* prototype for KINDense*/")
    562648           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
    563649           ))
    564650      ((ida)
    565651       (pp indent
     652           ("#include <sundials/sundials_types.h> /* definition of type realtype */")
     653           ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
    566654           ("#include <ida/ida.h>                 /* prototypes for IDA fcts., consts. */")
    567655           ("#include <ida/ida_dense.h>")
    568            ("#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */")
     656           ("#include <kinsol/kinsol.h>           /* prototypes for KINSOL fcts. */")
     657           ("#include <kinsol/kinsol_dense.h>     /* prototype for KINDense */")
    569658           ("#define Ith(v,i)    NV_Ith_S(v,i)    /* Ith component in a vector */")
    570659           ))
     
    588677      ((ida)
    589678       (pp indent ,ida-prelude-header)))
     679
     680    (if (not (null? steady-state-index-map))
     681        (case method
     682          ((cvode ida)
     683           (pp indent ,fsolve-prelude-header))))
    590684     
    591685 
     
    619713      )
    620714
     715    (if (not (null? steady-state-index-map))
     716        (case method
     717          ((cvode ida)
     718           (pp indent (extern "\"C\""
     719                              int ,(sprintf "~A_steadystate" sysname)
     720                              "(N_Vector, N_Vector, void*)" #\;) ,nl)
     721           ))
     722          )
     723
    621724    (pp indent
    622725        (,(sprintf "class ~A:" sysname) "public Archiving_Node { "))
     
    697800                         "(double, const double*, double*, void*)" #\;) ,nl))
    698801    )
     802
     803 
     804  (if (not (null? steady-state-index-map))
     805      (case method
     806        ((cvode ida)
     807         (pp indent (friend  int ,(sprintf "~A_steadystate" sysname)
     808                            "(N_Vector, N_Vector, void*)" #\;) ,nl)
     809         ))
     810      )
    699811
    700812  (let ((pscs (lookup-def 'post-synaptic-conductances synapse-info)))
     
    14391551
    14401552(define (output-update
    1441          sysname method state-index-map steady-state-index-map
     1553         sysname method state-index-map
    14421554         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    14431555         reaction-eq-defs transient-event-defs
     
    15041616
    15051617(define (output-nodes-leapfrog
    1506          sysname method state-index-map steady-state-index-map
     1618         sysname method state-index-map
    15071619         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    15081620         reaction-eq-defs i-eqs pool-ions perm-ions
     
    15701682
    15711683(define (output-nodes-cvode
    1572          sysname method state-index-map steady-state-index-map
     1684         sysname method state-index-map
    15731685         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    15741686         reaction-eq-defs i-eqs pool-ions perm-ions
     
    16151727
    16161728    /* Free integrator memory */
    1617 //    CVodeFree(&B_.sys_);
     1729    if (B_.sys_ != NULL)
     1730    {
     1731      CVodeFree(&B_.sys_);
     1732      B_.sys_ = NULL;
     1733    }
    16181734  }
    16191735}
     
    17491865
    17501866(define (output-nodes-ida
    1751          sysname method state-index-map steady-state-index-map
     1867         sysname method state-index-map
    17521868         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    17531869         reaction-eq-defs i-eqs pool-ions perm-ions
     
    17951911
    17961912    /* Free integrator memory */
    1797 //    IDAFree(&B_.sys_);
     1913    if (B_.sys_ != NULL)
     1914    {
     1915      IDAFree(&B_.sys_);
     1916      B_.sys_ = NULL;
     1917    }
     1918
    17981919  }
    17991920}
     
    19322053
    19332054(define (output-nodes-gsl
    1934          sysname method state-index-map steady-state-index-map
     2055         sysname method state-index-map
    19352056         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    19362057         reaction-eq-defs i-eqs pool-ions perm-ions
     
    20722193
    20732194(define (output-buffers
    2074          sysname method state-index-map steady-state-index-map
     2195         sysname method state-index-map
    20752196         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    20762197         reaction-eq-defs i-eqs pool-ions perm-ions
     
    23782499     c-eqs)
    23792500
     2501
     2502    (if (not (null? steady-state-index-map))
     2503
     2504        (case method
     2505
     2506          ((cvode ida)
     2507
     2508           (begin
     2509
     2510             (pp indent+
     2511                 ("N_Vector ssvect;")
     2512                 ,(expr->string/C++ (sprintf "N_VNew_Serial(~A);" (length steady-state-index-map)) 'ssvect)
     2513                 (,(expr->string/C++ `(fsolve ,(sprintf "~A_steadystate" sysname) ssvect "(void *)&p" ,(sprintf "\"~A\"" sysname) )) #\;))
     2514             
     2515             (for-each
     2516              (lambda (def)
     2517                (let* ((n (first def))
     2518                       (si (lookup-def n steady-state-index-map))
     2519                       (ni (lookup-def n state-index-map)))
     2520                  (if si
     2521                        (pp indent+
     2522                            ,(expr->string/C++ (ith 'ssvect si)
     2523                                               (sprintf "y_[~A]" ni))
     2524                            ,(expr->string/C++ (ith 'ssvect si)
     2525                                               n))
     2526                        )
     2527                  ))
     2528              rate-eq-defs)
     2529             
     2530             (pp indent+ ("N_VDestroy_Serial (ssvect);"))
     2531
     2532             ))
     2533          ))
    23802534   
    2381 #|
    2382     (if (not (null? steady-state-index-map))
    2383         (begin
    2384           (gsl-fminimizer sysname N indent+)
    2385           (for-each
    2386            (lambda (def)
    2387              (let* ((n (first def))
    2388                     (si (lookup-def n steady-state-index-map))
    2389                     (ni (lookup-def n state-index-map)))
    2390                (if si (begin
    2391                         (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys, " si ")") n))
    2392                         (pp indent+ ,(expr->string/C++ (s+ "gsl_vector_get(ys," si ")")
    2393                                                        (s+ "y_[" ni "]")))))))
    2394            rate-eq-defs)
    2395       (pp indent+ "gsl_vector_free (ys);")
    2396       ))
    2397 |#
    2398 
    23992535    (for-each
    24002536     (lambda (def)
     
    24872623
    24882624(define (output-recordables
    2489          sysname state-index-map steady-state-index-map
     2625         sysname state-index-map
    24902626         const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    24912627         reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
     
    27622898
    27632899
     2900(define (output-steadystate sysname method state-index-map steady-state-index-map
     2901                             imports external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     2902                             reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints
     2903                             indent indent+)
     2904
     2905 
     2906  (let* ((N (length steady-state-index-map))
     2907
     2908         (c-eqs (lookup-def 'c-eqs constraints))
     2909
     2910         (c-units (map (lambda (x) (let ((n (first x)) (v (second x)))
     2911                                     (list (nest-name n) v)))
     2912                       (lookup-def 'c-units constraints)))
     2913
     2914         (i-eqs
     2915          (map
     2916           (lambda (def) (list (first def)
     2917                               (add-params-to-fncall
     2918                                (canonicalize-expr/C++ (second def)) builtin-fns)))
     2919           i-eqs))
     2920
     2921         (init-eqs
     2922          (append
     2923             
     2924           (map (lambda (def)
     2925                  (let ((n (first def))
     2926                        (b (second def)))
     2927                    (list (nest-name n) (nest-name b))))
     2928                external-eq-defs)
     2929           
     2930           asgn-eq-defs
     2931           init-eq-defs
     2932
     2933           (map (lambda (pool-ion)
     2934                  (let ((n (pool-ion-in pool-ion))
     2935                        (b (pool-ion-inq pool-ion)))
     2936                    (list n b)))
     2937                pool-ions)
     2938
     2939           (map (lambda (pool-ion)
     2940                  (let ((n (pool-ion-out pool-ion))
     2941                        (b (pool-ion-outq pool-ion)))
     2942                    (list n b)))
     2943                pool-ions)
     2944           ))
     2945
     2946         (init-dag
     2947          (map (lambda (def)
     2948                 (cons (first def) (enum-freevars (second def) '() '())))
     2949               init-eqs))
     2950
     2951         (init-order
     2952          (reverse
     2953           (topological-sort init-dag
     2954                             (lambda (x y) (string=? (->string x) (->string y))))))
     2955
     2956         (init-locals  (find-locals (map second (append init-eqs i-eqs reaction-eq-defs))))
     2957
     2958         )
     2959   
     2960    (pp indent ,nl
     2961        (,(sprintf "extern \"C\" int ~A_steadystate" sysname) "(N_Vector u, N_Vector f, void* pnode)" ,nl)
     2962        ( #\{))
     2963
     2964    (pp indent+
     2965        (double ,(slp ", " (delete-duplicates
     2966                            (map ->string
     2967                                 (filter (lambda (x) (not (member x builtin-consts)))
     2968                                         (append
     2969                                          init-locals
     2970                                          init-order
     2971                                          (map first external-eq-defs)
     2972                                          (map pool-ion-in pool-ions)
     2973                                          (map pool-ion-out pool-ions)
     2974                                          (map first i-eqs)
     2975                                          (map first steady-state-index-map)
     2976                                          (map first const-defs)
     2977                                          (map first reaction-eq-defs)
     2978                                          )))
     2979                                 string=?))
     2980                      ";")
     2981        )
     2982
     2983    (pp indent+ ,nl
     2984
     2985        ("// params is a reference to the model parameters ")
     2986        (const struct ,(sprintf "~A::Parameters_" sysname)
     2987               "*params = "
     2988               "(" ,(sprintf "struct ~A::Parameters_ *" sysname) ")" "pnode;")
     2989       
     2990        ,nl
     2991        )
     2992
     2993    (for-each
     2994     (lambda (x)
     2995       (let* ((n  (first x)))
     2996         (pp indent+ ,(expr->string/C++ (sprintf "params->~A" n) n))))
     2997     const-defs)
     2998
     2999    (for-each
     3000     (lambda (def)
     3001       (let* ((n   (first def))
     3002              (ni  (lookup-def n steady-state-index-map)))
     3003         (if ni (pp indent+ ,(expr->string/C++ (ith 'u ni) n)))
     3004         ))
     3005     rate-eq-defs)
     3006
     3007    (for-each (lambda (def)
     3008                (pp indent+ ,(expr->string/C++ 0. (first def))))
     3009              i-eqs)
     3010
     3011    (for-each (lambda (n)
     3012                (let ((b  (lookup-def n init-eqs)))
     3013                  (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
     3014              init-order)
     3015   
     3016    (for-each
     3017     (lambda (def)
     3018       (let* ((n   (first def))
     3019              (ni  (lookup-def n steady-state-index-map))
     3020              (b   (second def))
     3021              (lb  (find-locals (list b))))
     3022         (if ni
     3023             (begin
     3024               (if (not (null? lb))
     3025                   (pp indent+ (double ,(slp ", " (delete-duplicates lb)) #\;)))
     3026               (pp indent+
     3027                   (,(expr->string/C++ b (ith 'f ni))))))
     3028         ))
     3029     rate-eq-defs)
     3030
     3031    (pp indent ,nl
     3032        ( "return 0;" )
     3033        ( #\}))
     3034
     3035    ))
     3036
     3037
    27643038(define (nemo:nest-translator sys . rest)
    27653039
     
    29413215            (lambda ()
    29423216              (pp indent ,nl ("/* " This file was generated by ,(nemo:version-string) on ,(seconds->string (current-seconds)) " */" ,nl))
    2943               (output-prelude sysname steady-state-index-map indent)
     3217              (output-prelude sysname indent)
    29443218              ))
    29453219         
     
    29523226               (lambda ()
    29533227                 (pp indent ,sundials-prelude)))))
     3228
     3229
     3230          (if (not (null? steady-state-index-map))
     3231              (case method
     3232                ((cvode ida)
     3233                 (with-output-to-port cpp-output
     3234                   (lambda ()
     3235                     (pp indent ,fsolve-prelude))))))
    29543236
    29553237          ;; user-defined functions
     
    29953277              ))
    29963278
     3279          ;; steady state function (used by CVODE and IDA)
     3280          (with-output-to-port cpp-output
     3281            (lambda ()
     3282              (output-steadystate sysname method state-index-map steady-state-index-map
     3283                                  imports-sans-v external-eq-defs const-defs asgn-eq-defs init-eq-defs rate-eq-defs
     3284                                  reaction-eq-defs i-eqs pool-ions perm-ions defaults constraints indent indent+)
     3285              (pp indent ,nl)
     3286              ))
     3287
    29973288          ;; Jacobian function
    29983289          (with-output-to-port cpp-output
     
    30053296            (lambda ()
    30063297              (output-recordables
    3007                sysname state-index-map steady-state-index-map
     3298               sysname state-index-map
    30083299               const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    30093300               reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
     
    30273318              ))
    30283319         
     3320         
    30293321          ;; accessors and modifiers for states and parameters
    30303322          (with-output-to-port cpp-output
     
    30443336            (lambda ()
    30453337              (output-buffers
    3046                sysname method state-index-map steady-state-index-map
     3338               sysname method state-index-map
    30473339               const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    30483340               reaction-eq-defs i-eqs pool-ions perm-ions
     
    30653357                  output-nodes-ida)
    30663358                 )
    3067                sysname method state-index-map steady-state-index-map
     3359               sysname method state-index-map
    30683360               const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    30693361               reaction-eq-defs i-eqs pool-ions perm-ions
     
    31113403            (lambda ()
    31123404              (output-update
    3113                sysname method state-index-map steady-state-index-map
     3405               sysname method state-index-map
    31143406               const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    31153407               reaction-eq-defs transient-event-defs i-eqs
     
    31253417            (lambda ()
    31263418              (output-event-handle
    3127                sysname state-index-map steady-state-index-map
     3419               sysname state-index-map
    31283420               const-defs asgn-eq-defs init-eq-defs rate-eq-defs
    31293421               reaction-eq-defs i-eqs pool-ions perm-ions
  • release/4/nemo/trunk/nemo-nmodl.scm

    r29862 r29931  
    11511151
    11521152          (if has-kinetic?
    1153               (pp indent+ (SOLVE kstates STEADYSTATE sparse)))
     1153              (pp indent+
     1154                  (SOLVE kstates STEADYSTATE sparse)
     1155                  ("reactions ()")
     1156                ))
    11541157         
    11551158          (for-each
Note: See TracChangeset for help on using the changeset viewer.