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

Last change on this file since 25514 was 25514, checked in by Ivan Raikov, 9 years ago

9ML-toolkit: support for lsode solver in octave

File size: 14.1 KB
Line 
1;;
2;;  An IVP solver for NineML.
3;;
4;;
5;; Copyright 2010-2011 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(require-extension setup-api srfi-13 datatype static-modules miniML miniMLsyntax miniMLvalue miniMLeval)
23(require-extension getopt-long ssax sxml-transforms sxpath sxpath-lolevel object-graph signal-diagram)
24(require-extension 9ML-parse 9ML-repr )
25(require-extension 9ML-ivp-octave 9ML-ivp-chicken 9ML-ivp-mlton 9ML-ivp-octave-mlton )
26       
27
28
29(include "SXML.scm")
30(include "SXML-to-XML.scm")
31
32(define-values (env-binding? env-empty env-add-signature env-add-module env-add-type env-add-spec env-add-value
33                env-find-value env-find-type env-find-module env-find)
34  (make-mod-env core-syntax))
35
36(define-values (scope-typedecl scope-modtype scope-signature scope-modterm scope-moddef)
37  (make-mod-scoping core-syntax core-scoping))
38
39(define-values (check-modtype check-signature type-modterm type-moddef type-definition)
40  (make-mod-typing core-syntax core-typing))
41
42
43
44(include "NineMLcore.scm")
45(include "NineMLsignal.scm")
46(include "NineMLdiagram.scm")
47(include "NineMLinterval.scm")
48(include "NineMLgraph.scm")
49(include "NineMLivp.scm")
50
51
52(define init-scope      (make-parameter st-empty))
53(define init-type-env   (make-parameter env-empty))
54(define init-eval-env   (make-parameter env-empty))
55
56
57(define (enter-typedecl id decl)
58  (init-scope (st-enter-type id (init-scope)))
59  (init-type-env   (env-add-type id decl (init-type-env))))
60
61(define (enter-valtype name ty)
62  (let ((id (ident-create name)))
63    (init-scope (st-enter-value id (init-scope)))
64    (init-type-env   (env-add-value id ty (init-type-env)))))
65
66(define (enter-val name val)
67  (let ((id (or (and (ident? name) name) (ident-create name))))
68    (init-eval-env (ident-add id val (init-eval-env)))))
69
70(core-initialize enter-typedecl enter-valtype)
71(eval-cbv-initialize enter-val)
72
73
74(define (enter-module id mty)
75  (init-scope (st-enter-module id (init-scope)))
76  (init-type-env (env-add-module id mty (init-type-env))))
77
78
79(define lookup-def 
80  (lambda (k lst . rest)
81    (let-optionals rest ((default #f))
82      (alist-ref k lst eq? default))))
83
84
85(define opt-defaults
86  `(
87    (platform . chicken)
88    ))
89
90(define (defopt x)
91  (lookup-def x opt-defaults))
92
93                 
94(define opt-grammar
95  `(
96
97    (print-type-env  "prints the type environment of each operand"
98                     (single-char #\t)
99                     (value (optional COMPONENT-LIST)
100                            (default all)
101                            (transformer 
102                             ,(lambda (x) 
103                                (if (string=? x "all") x
104                                    (list (string-split x ",")))))))
105
106    (print-eval-env  "prints the evaluation environment of each operand"
107                     (single-char #\e)
108                     (value (optional COMPONENT-LIST)
109                            (default all)
110                            (transformer 
111                             ,(lambda (x) 
112                                (if (string=? x "all") x
113                                    (list (string-split x ",")))))))
114
115    (print-source-defs  "prints the source definitions of each operand"
116                        (single-char #\s))
117
118    (html-report        "prints out an HTML report of the unified environments of each operand")
119
120    (output-sxml        "sets output format to SXML")
121
122    (output-xml         "sets output format to XML")
123
124    (data            "save data from simulations in files ${OPERAND}_NAME.log"
125                     (single-char #\d))
126
127    (plot            "save PNG plots of simulation data in files ${OPERAND}_NAME.png"
128                     (single-char #\p))
129
130    (plot-eps        "save EPS plots of simulation data in files ${OPERAND}_NAME.eps")
131
132    (plot-index      "specify index of variable to plot"
133                     (value (required INDEX)
134                            (predicate  ,(lambda (x) (integer? (string->number x))))
135                            (transformer ,string->number)
136                             ))
137
138    (xml            "reads canonical NineML XML representation of each operand"
139                    (single-char #\x))
140
141    (platform        "simulation platform (one of chicken, mlton, octave, octave/mlton)"
142                     (value (required PLATFORM)
143                            (predicate 
144                             ,(lambda (x) 
145                                (let ((s (string->symbol (string-downcase x))))
146                                  (case s
147                                    ((chicken mlton octave octave/mlton) s)
148                                    (else (error 'ivp "unrecognized platform" x))))))
149                            (transformer ,string->symbol)
150                             ))
151
152    (verbose          "print commands as they are executed"
153                      (single-char #\v))
154
155
156    (help         (single-char #\h))           
157
158    ))
159
160
161;; Process arguments and collate options and arguments into OPTIONS
162;; alist, and operands (filenames) into OPERANDS.  You can handle
163;; options as they are processed, or afterwards.
164
165(define opts    (getopt-long (command-line-arguments) opt-grammar))
166(define opt     (make-option-dispatch opts opt-grammar))
167
168(define simulation-platform (make-parameter #f))
169(define ivp-verbose (make-parameter 0))
170
171(define (d fstr . args)
172  (let ([port (current-error-port)])
173    (if (positive? (ivp-verbose)) 
174        (begin (apply fprintf port fstr args)
175               (flush-output port) ) )))
176
177(define (run:execute explist)
178  (define (smooth lst)
179    (let ((slst (map ->string lst)))
180      (string-intersperse (cons (car slst) (cdr slst)) " ")))
181  (for-each (lambda (cmd) (system (->string cmd)))
182            (map smooth explist)))
183
184(define (run:execute* explist)
185  (define (smooth lst)
186    (let ((slst (map ->string lst)))
187      (string-intersperse (cons (car slst) (cdr slst)) " ")))
188  (for-each (lambda (cmd) (system* "~a" cmd))
189            (map smooth explist)))
190
191(define-syntax run
192  (syntax-rules ()
193    ((_ exp ...)
194     (begin
195       (d "running ~A ...~%" (list `exp ...))
196       (run:execute* (list `exp ...))))))
197
198(define-syntax run-
199  (syntax-rules ()
200    ((_ exp ...)
201     (begin
202       (d "running ~A ...~%" (list `exp ...))
203       (run:execute (list `exp ...))))))
204
205
206;; Use args:usage to generate a formatted list of options (from OPTS),
207;; suitable for embedding into help text.
208(define (ivp:usage)
209  (print "Usage: " (car (argv)) " [options...] file1... ")
210  (newline)
211  (print "The following options are recognized: ")
212  (newline)
213  (print (parameterize ((indent 5)) (usage opt-grammar)))
214  (exit 1))
215
216
217(define nl "\n")
218
219
220
221(define (generate-ivp-table prefix ivp-id sxml-tuple #!key (platform 'chicken))
222
223  (let ((sexpr 
224         (let ((sexpr (sxml-value->sexpr sxml-tuple)))
225           (case (car sexpr)
226             ((ivp)
227              (and (pair? (cdr sexpr))
228                   (case (cadr sexpr)
229                     ((construct)  (cddr sexpr))
230                     (else #f))))
231             (else #f)))))
232      (and sexpr
233           (let* ((ivar    (cadr sexpr))
234                  (hvar    (caddr sexpr))
235                  (start   (cadddr sexpr))
236                  (end     (cadddr (cdr sexpr)))
237                  (diagram+initial (sexpr->diagram+initial hvar (car sexpr))))
238             (let ((sd       (construct (car diagram+initial)))
239                   (ic       (cadr diagram+initial)))
240               (if (not (alist-ref ivar ic))
241                   (error 'generate-ivp-table "IVP independent variable is not present in given system" ivar))
242               (if (not (alist-ref hvar ic))
243                   (error 'generate-ivp-table "IVP step variable is not present in given system" hvar))
244               (let* ((dfe (dataflow sd '()))
245                      (dvars (lset-difference eq?
246                                              (lset-intersection eq? (alist-ref 'in dfe) (alist-ref 'out dfe))
247                                              (list ivar)))
248                      (pvars (lset-difference eq? (alist-ref 'in dfe) (cons ivar dvars)))
249                      (events (events sd))
250                      )
251
252
253                 (case platform
254                   
255                   ((octave)
256                    (ivp-octave prefix ivp-id hvar ivar dvars pvars events start end ic sd)
257                    (list ivp-id ivar dvars) )
258
259                   ((octave/mlton octave-mlton)
260                    (ivp-octave-mlton prefix ivp-id hvar ivar dvars pvars start end ic sd)
261                    (list ivp-id ivar dvars) )
262                   
263                   ((mlton)
264                    (ivp-mlton  prefix ivp-id ivar dvars pvars start end ic sd)
265                    (list ivp-id ivar dvars) )
266
267                   ((chicken)
268                    (ivp-chicken  prefix ivp-id ivar dvars pvars start end ic sd)
269                    (list ivp-id ivar dvars) )
270
271
272                   (else (error 'generate-ivp-table "unknown platform"  platform))
273                   
274                   )))
275             )))
276  )
277
278
279(define (generate-ivp-plot prefix ivp-info #!key (format 'png) (index #f))
280  (let ((ivp-id (car ivp-info))
281        (dvars (caddr ivp-info)))
282    (let* ((n (+ 1 (length dvars)))
283           (range (if (and index (integer? index)) (number->string index) (sprintf "2:~A" n)))
284           (linewidth 3)
285           (dir (or (pathname-directory prefix) "."))
286           (log-path (make-pathname dir (sprintf "~A_~A.log" (pathname-file prefix) ivp-id)))
287           (png-path (make-pathname dir (sprintf "~A_~A.png" (pathname-file prefix) ivp-id)))
288           (eps-path (make-pathname dir (sprintf "~A_~A.eps" (pathname-file prefix) ivp-id)))
289           )
290      (case format
291             ((png)
292              (run (octave -q --eval 
293                           #\'
294                           log = load (,(sprintf "\"~A\"" log-path)) #\;
295                           h = plot ("log(:,1)" #\, ,(sprintf "log(:,~A)" range)) #\;
296                           ,@(if index 
297                                 (concatenate (list-tabulate 1 (lambda (i) `(set (h(,(+ 1 i)) #\, "\"linewidth\"" #\, ,linewidth) #\;))))
298                                 (concatenate (list-tabulate (- n 1) (lambda (i) `(set (h(,(+ 1 i)) #\, "\"linewidth\"" #\, ,linewidth) #\;)))))
299                           print (,(sprintf "\"~A\"" png-path) #\, "\"-dpng\"") #\;
300                           #\')))
301             ((eps)
302              (run (octave -q --eval 
303                           #\'
304                           log = load (,(sprintf "\"~A\"" log-path)) #\;
305                           h = plot ("log(:,1)" #\, ,(sprintf "log(:,~A)" range)) #\;
306                           ,@(if index 
307                                 (concatenate (list-tabulate 1 (lambda (i) `(set (h(,(+ 1 i)) #\, "\"linewidth\"" #\, ,linewidth) #\;))))
308                                 (concatenate (list-tabulate (- n 1) (lambda (i) `(set (h(,(+ 1 i)) #\, "\"linewidth\"" #\, ,linewidth) #\;)))))
309                           print (,(sprintf "\"~A\"" eps-path) #\, "\"-depsc2\"") #\;
310                           #\')))
311             (else (error 'generate-ivp-plot "unrecognized format" format)))
312      )))
313
314
315
316(define (make-ivp-hook #!key (ivp #f) (diagram #f) (ivp-plot #f) (plot-format 'png) (plot-index #f))
317  (lambda (prefix name label value)
318    (cond
319     ((and diagram
320           (or (and (string? label) (string=? label "diagram")) 
321               (and (pair? label) (string=? (car label) "diagram")))) ;; value is a diagram
322      (let* ((diagram-id (gensym 'diagram))
323             (diagram-link `(img (@ (src ,(sprintf "~A.png" diagram-id))) (alt "NineML diagram"))))
324        (generate-diagram prefix diagram-id value)
325        diagram-link
326        ))
327     
328     ((and ivp (or (and (string? label) (string=? label "ivp"))
329                   (and (pair? label) (string=? (car label) "ivp")))) ;; value is an IVP
330      (let* ((ivp-id (gensym 'ivp))
331             (ivp-plot-link `(img (@ (src ,(sprintf "~A_~A.png" (pathname-file prefix) ivp-id)) (alt "NineML IVP plot"))))
332             (ivp-info (generate-ivp-table prefix ivp-id value platform: (simulation-platform))))
333        (if ivp-plot (generate-ivp-plot prefix ivp-info format: plot-format index: plot-index))
334        (and ivp-plot ivp-plot-link)
335        ))
336     
337     (else #f))))
338           
339
340
341
342(define (parse-xml fpath)
343  (with-input-from-file fpath
344    (lambda () (cons '*TOP* (ssax:xml->sxml (current-input-port) `((nml . ,nineml-xmlns) (nml . "CoModL")))))
345    ))
346
347
348(define (interpreter operand #!key (xml #f))
349
350  (let ((defs (if xml (parse-al-sxml (parse-xml operand))
351                  (call-with-input-file operand (lambda (in) (parse 'NineML in))))))
352
353
354    (let* ((scoped-defs      (scope-moddef (init-scope) defs))
355           (mty              (type-moddef (init-type-env) '() scoped-defs))
356           (type-env         (map (lambda (x) (cases modspec x
357                                                     (Value_sig (id vty) (cons id x))
358                                                     (Type_sig (id decl) (cons id x))
359                                                     (Module_sig (id mty) (cons id x))
360                                                     )) mty))
361           (eval-env         (mod-eval-cbv (init-eval-env) scoped-defs))
362           (unified-env      (list scoped-defs 
363                                   (filter (lambda (x) (not (assoc (car x) (init-type-env)))) type-env) 
364                                   (filter (lambda (x) (not (assoc (car x) (init-eval-env)))) eval-env) ))
365
366           )
367      unified-env
368      )))
369
370
371(define (main options operands)
372
373  (if (options 'help) (ivp:usage))
374
375  (let ((find-module (lambda (x) (env-find-module x (init-type-env)))))
376    (for-each (lambda (init name) (init name enter-module find-module init-eval-env))
377              (list Signal:module-initialize   
378                    Diagram:module-initialize 
379                    Interval:module-initialize 
380                    Graph:module-initialize
381                    IVP:module-initialize )
382              (list "Signal" "Diagram" "Interval" "Graph" "IVP" )) )
383
384  (if (null? operands)
385      (ivp:usage)
386      (let ((output-type (cond ((options 'output-xml)  'xml)
387                               ((options 'output-sxml) 'sxml)
388                               (else #f)))
389            (unified-envs (map (lambda (x) (interpreter x xml: (options 'xml))) operands)))
390        (if (options 'verbose) (begin (repr-verbose 1) (ivp-verbose 1)))
391        (simulation-platform (or (options 'platform) (defopt 'platform) ))
392        (for-each
393         (lambda (operand uenv)
394           (let ((source-defs (car uenv))
395                 (mty         (cadr uenv))
396                 (eval-env    (caddr uenv)))
397             
398             (let ((type-env-opt (options 'print-type-env)))
399               (if type-env-opt
400                   (if (and (string? type-env-opt) (string=?  type-env-opt "all"))
401                     (print-type-env mty output-type)
402                     (let ((fc (lambda (x) (and (member (ident-name (car x)) type-env-opt) x))))
403                       (print-type-env mty output-type fc)))
404                   ))
405
406             (let ((eval-env-opt (options 'print-eval-env)))
407               (if eval-env-opt
408                   (if (and (string? eval-env-opt) (string=? eval-env-opt "all"))
409                     (print-eval-env eval-env)
410                     (let ((fc (lambda (x) (and (member (ident-name (car x)) eval-env-opt) x))))
411                       (print-eval-env eval-env output-type fc)))
412                   ))
413
414             (if (options 'print-source-defs)
415                 (print-source-defs source-defs output-type))
416
417             (if (options 'html-report)
418                 (html-report operand uenv value-hook: (make-ivp-hook diagram: #t ivp: #t ivp-plot: #t)))
419
420             (if (options 'data)
421                 (traverse-definitions operand uenv value-hook: (make-ivp-hook ivp: #t)))
422
423             (if (options 'plot)
424                    (traverse-definitions operand uenv value-hook: (make-ivp-hook ivp: #t ivp-plot: #t plot-index: (options 'plot-index))))
425
426             (if (options 'plot-eps)
427                 (traverse-definitions operand uenv value-hook: (make-ivp-hook ivp: #t ivp-plot: #t plot-format: 'eps plot-index: (options 'plot-index))))
428
429             ))
430         operands unified-envs))))
431
432(width 40)
433(main opt (opt '@))
Note: See TracBrowser for help on using the repository browser.