source: project/release/4/9ML-toolkit/trunk/ivp-octave.scm @ 29952

Last change on this file since 29952 was 29952, checked in by Ivan Raikov, 8 years ago

9ML-toolkit: bringing octave and scheme backends up to date; changed Izhikevich FS example to use heaviside function

File size: 5.4 KB
Line 
1;;
2;;  NineML IVP code generator for Octave.
3;;
4;;
5;; Copyright 2010-2013 Ivan Raikov and the Okinawa Institute of
6;; Science and Technology.
7;;
8;; This program is free software: you can redistribute it and/or
9;; modify it under the terms of the GNU General Public License as
10;; published by the Free Software Foundation, either version 3 of the
11;; License, or (at your option) any later version.
12;;
13;; This program is distributed in the hope that it will be useful, but
14;; WITHOUT ANY WARRANTY; without even the implied warranty of
15;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16;; General Public License for more details.
17;;
18;; A full copy of the GPL license can be found at
19;; <http://www.gnu.org/licenses/>.
20;;
21
22
23(module 9ML-ivp-octave
24
25        (ivp-octave ivp-octave-codegen)
26
27        (import scheme chicken )
28
29(import (only extras sprintf)
30        (only files make-pathname pathname-directory pathname-file absolute-pathname? )
31        (only data-structures conc alist-ref intersperse ->string)
32        (only posix current-directory)
33        (only srfi-1 filter list-index)
34        (only srfi-13 string-concatenate))
35
36(require-extension make datatype signal-diagram 9ML-eval setup-api)
37
38
39(define nl "\n")
40
41       
42(define (octave-value v)
43  (cond
44   ((pair? v)
45    (case (car v)
46      ((realsig)   (octave-value (caddr v)))
47      ((realconst) (octave-value (cadr v)))
48      ((generator) (sprintf "~A ()" (cadr v)))
49      ((random)    (sprintf "~A ()" (cadr v)))
50      ((neg)       (sprintf "-(~A))" (octave-value (cadr v))))
51      ((+ - * / >= <= > <) 
52       (sprintf "(~A ~A ~A)"
53                (octave-value (cadr v)) 
54                (car v) 
55                (octave-value (caddr v))))
56      ((log ln sin cos cosh tanh exp)
57       (sprintf "~A (~A)"
58                (car v) (octave-value (cadr v)) ))
59      (else (error 'octave-value "invalid value" v))))
60   ((boolean? v)  (if v "true" "false"))
61   (else (sprintf "~A" v))))
62
63(define (octave-m ivp-id start ivar hvar dvars events ic imax log-path adir solver-path)
64
65  (let* ((states (cons ivar dvars))
66         (parameters (filter (lambda (x) (not (member (car x) states))) ic)))
67
68  `(,(sprintf "autoload (~S,~S);~%" (->string ivp-id) (make-pathname adir solver-path))
69    ,(sprintf "~A;~%~%" ivp-id)
70   
71    ,(sprintf "~A_initial = struct(~A);~%~%" ivp-id
72              (string-concatenate
73               (intersperse
74                (map (lambda (x)
75                       (let ((n (car x)) (v (octave-value (cdr x))))
76                         (sprintf "\"~A\",~A" n v)))
77                     (cons (cons ivar start) ic))
78                ",")))
79
80    ,(sprintf "i = 1; t = 0; tmax = ~A;~%" imax)
81    ,(sprintf "~A_state = ~A_initial;~%" ivp-id ivp-id)
82    ,(sprintf "data = [[~A]];~%" 
83              (string-concatenate (intersperse (map (lambda (s) (sprintf "~A_state.~A" ivp-id s)) states) ",") ))
84
85    ,(sprintf "tic;~%")
86    ,(sprintf "while (t < tmax) ~%")
87
88    ,(sprintf "  y = ~A(~A_state);~%" ivp-id ivp-id)
89
90    ,(if (null? events) '()
91         `(,(sprintf "  if (~A)~%" (intersperse (map (lambda (x) (sprintf "y.~A" x)) events) " || "))
92           ,(sprintf "    ~A_state.~A = 1e-4;~%" ivp-id hvar)
93           ,(sprintf "    y = ~A(~A_state);~%" ivp-id ivp-id)
94           ,(sprintf "  else~%")
95           ,(sprintf "    if (~A_state.~A < 0.25)~%" ivp-id hvar)
96           ,(sprintf "      ~A_state.~A = ~A_state.~A+1e-2;~%" ivp-id hvar ivp-id hvar)
97           ,(sprintf "    endif~%")
98           ,(sprintf "  endif~%")
99           ))
100
101    ,(map (lambda (s) (sprintf "  ~A_state.~A = y.~A;~%" ivp-id s s)) states) 
102    ,(sprintf "  data = [data; [~A]];~%" 
103              (string-concatenate (intersperse (map (lambda (s) (sprintf "y.~A" s)) states) ",") ))
104   
105    ,(sprintf "  t  = y.~A;~%" ivar )
106    ,(sprintf "  i = i+1;~%")
107    ,(sprintf "endwhile~%")
108    ,(sprintf "toc;~%")
109
110    ,(sprintf "save (\"-ascii\", ~S, \"data\");~%" log-path)
111    )
112  ))
113
114
115
116
117(define (ivp-octave prefix ivp-id hvar ivar dvars pvars events start end ic sd)
118
119  (let* (                       
120         (N (+ 1 (length dvars)))
121         (P (- (length ic) N))
122         (dir        (or (pathname-directory prefix) "."))
123         (adir       (if (absolute-pathname? dir) dir
124                         (make-pathname (current-directory) dir)))
125         (shared-dir         (chicken-home))
126         (solver-path        (make-pathname dir (sprintf "~A_solver.m" ivp-id)))
127         (m-path             (make-pathname dir (sprintf "~A_~A.m" (pathname-file prefix) ivp-id)))
128         (log-path           (make-pathname dir (sprintf "~A_~A.log" (pathname-file prefix) ivp-id)))
129         (octave-path        "octave")
130         )
131   
132    (make (
133           (solver-path (prefix)
134                        (with-output-to-file solver-path
135                          (lambda () 
136                            (codegen/Octave ivp-id sd solver: 'lsode))) )
137           
138           (m-path ()
139                   (with-output-to-file m-path
140                     (lambda () 
141                       (print-fragments
142                        `(,(octave-m ivp-id start ivar hvar dvars events ic end log-path adir solver-path))))) )
143           
144           
145           (log-path (m-path solver-path)
146                     (run (octave -q ,m-path)) )
147
148#|         
149           (clean ()
150                  (run rm ,solver-path ,clib-path ,h-path
151                       ,mlb-path ,clib-so-path ,lib-so-path
152                       ,oct-cc-path ,oct-run-path ,oct-initial-path
153                       ,oct-open-path ,oct-close-path ,m-path ,log-path))
154|#
155           
156           )
157      (list solver-path m-path log-path) )))
158
159
160(define (ivp-octave-codegen prefix ivp-id hvar ivar dvars pvars events ic sd)
161
162  (let* (                       
163         (N (+ 1 (length dvars)))
164         (P (- (length ic) N))
165         (dir          (or (pathname-directory prefix) "."))
166         (solver-path  (make-pathname dir (sprintf "~A_solver.m" ivp-id)))
167         )
168   
169    (make (
170           (solver-path (prefix)
171                        (with-output-to-file solver-path
172                          (lambda () 
173                            (codegen/Octave ivp-id sd solver: 'lsode))) )
174           )
175      (list solver-path) )))
176
177)
178
Note: See TracBrowser for help on using the repository browser.