source: project/release/4/picnic/trunk/picnic-utils.scm @ 30645

Last change on this file since 30645 was 30645, checked in by Ivan Raikov, 7 years ago

picnic: reintegrated calculator language parser

File size: 35.8 KB
Line 
1;;
2;; Neural Parametric Curve Connectivity spatial and geometric utility procedures.
3;;
4;; Copyright 2012-2014 Ivan Raikov and the Okinawa Institute of Science and Technology
5;;
6;; This program is free software: you can redistribute it and/or
7;; modify it under the terms of the GNU General Public License as
8;; published by the Free Software Foundation, either version 3 of the
9;; License, or (at your option) any later version.
10;;
11;; This program is distributed in the hope that it will be useful, but
12;; WITHOUT ANY WARRANTY; without even the implied warranty of
13;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14;; General Public License for more details.
15;;
16;; A full copy of the GPL license can be found at
17;; <http://www.gnu.org/licenses/>.
18;;
19
20(module picnic-utils
21
22        *
23
24
25        (import scheme chicken)
26
27       
28        (require-extension srfi-69 datatype matchable vector-lib 
29                           mpi mathh typeclass kd-tree random-mtzig
30                           lalr-driver)
31
32
33        (require-library srfi-1 srfi-4 srfi-13 irregex files posix data-structures
34                         bvsp-spline parametric-curve)
35
36
37        (import 
38                (only srfi-1 
39                      fold filter-map filter every zip list-tabulate delete-duplicates partition 
40                      first second third take)
41                (only srfi-4 
42                      s32vector s32vector-length s32vector-ref
43                      f64vector f64vector? f64vector-ref f64vector-length f64vector->list list->f64vector)
44                (only srfi-13 string= string< string-null? string-prefix? string-trim-both)
45                (only irregex string->irregex irregex-match)
46                (only files make-pathname)
47                (only posix glob)
48                (only extras read-lines pp fprintf )
49                (only ports with-output-to-port )
50                (only data-structures ->string alist-ref compose identity string-split sort atom?)
51                (only lolevel extend-procedure procedure-data extended-procedure?)
52                (prefix bvsp-spline bvsp-spline:)
53                (prefix parametric-curve pcurve:)
54                )
55
56        (define picnic-verbose (make-parameter 0))
57
58
59        (define (d fstr . args)
60          (let ([port (current-error-port)])
61            (if (positive? (picnic-verbose)) 
62                (begin (apply fprintf port fstr args)
63                       (flush-output port) ) )))
64
65
66        (include "mathh-constants")
67
68        (define (sign x) (if (negative? x) -1.0 1.0))
69
70        (define (f64vector-empty? x) (zero? (f64vector-length x)))
71
72        (define (random-init seed) (random-mtzig:init seed))
73
74        (define (random-uniform low high st)
75          (let ((rlo (if (< low high) low high))
76                (rhi (if (< low high) high low))) 
77            (let ((delta (+ 1 (- rhi rlo)))
78                  (v (random-mtzig:randu! st)))
79              (+ rlo (floor (* delta v)))
80              ))
81          )
82
83
84        (define (random-normal mean sdev st)
85          (+ mean (* sdev (random-mtzig:randn! st))))
86
87        (import-instance (<KdTree> KdTree3d)
88                         (<Point> Point3d))
89       
90        ;; convenience procedure to access to results of kd-tree-nearest-neighbor queries
91        (define (kdnn-point x) (cadr x))
92        (define (kdnn-distance x) (caddr x))
93        (define (kdnn-index x) (caar x))
94        (define (kdnn-parent-index x) (car (cadar x)))
95        (define (kdnn-parent-distance x) (cadr (cadar x)))
96
97
98        (define point->list f64vector->list)
99       
100
101       
102        (define-record-type pointset (make-pointset prefix id points boundary)
103          pointset? 
104          (prefix pointset-prefix)
105          (id pointset-id)
106          (points pointset-points)
107          (boundary pointset-boundary)
108          )
109
110
111       
112
113        (define-record-type genpoint (make-genpoint coords parent-index parent-distance segment)
114          genpoint? 
115          (coords genpoint-coords)
116          (parent-index genpoint-parent-index)
117          (parent-distance genpoint-parent-distance)
118          (segment genpoint-segment)
119          )
120
121       
122        (define-record-type cell (make-cell ty index origin sections)
123          cell? 
124          (ty cell-type)
125          (index cell-index)
126          (origin cell-origin)
127          (sections cell-sections)
128          )
129       
130       
131        (define (cell-section-ref name cell)
132          (alist-ref name (cell-sections cell)))
133
134       
135        (define (write-pointset name pts)
136            (call-with-output-file (string-append (->string name) ".sorted.dat")
137              (lambda (out)
138                (for-each (lambda (x) 
139                            (fprintf out "~A ~A ~A~%" 
140                                     (coord 0 x)
141                                     (coord 1 x)
142                                     (coord 2 x)))
143                          pts))
144              ))
145
146        (define (write-sections section-name cells)
147          (lambda (out)
148            (for-each
149             (lambda (gx) 
150               (for-each
151                (lambda (section)
152                  (for-each
153                   (lambda (gd)
154                     (let ((p (genpoint-coords gd)))
155                       (fprintf out "~A ~A ~A "
156                                (coord 0 p)
157                                (coord 1 p)
158                                (coord 2 p))))
159                   (cdr section)))
160                (cell-section-ref section-name gx))
161               (fprintf out "~%"))
162             cells)))
163
164
165        (define (cells-sections->kd-tree cells section-name #!key
166                                         (make-value
167                                          (lambda (i v) 
168                                            (list (genpoint-parent-index v)
169                                                  (genpoint-parent-distance v))))
170                                         (make-point
171                                          (lambda (v) (genpoint-coords v))))
172          (let ((t 
173                 (let recur ((cells cells) (points '()))
174                   (if (null? cells)
175                       (list->kd-tree
176                        points
177                        make-value: make-value
178                        make-point: make-point)
179                       (let ((cell (car cells)))
180                         (recur (cdr cells) 
181                                (let inner ((sections (append (cell-section-ref section-name cell)))
182                                            (points points))
183                                  (if (null? sections) points
184                                      (inner
185                                       (cdr sections)
186                                       (append (cdr (car sections)) points))
187                                      ))
188                                ))
189                       ))
190                 ))
191            t))
192
193        (define (sections->kd-tree cells #!key
194                                   (make-value
195                                    (lambda (i v) 
196                                      (list (genpoint-parent-index v)
197                                            (genpoint-parent-distance v))))
198                                   (make-point
199                                    (lambda (v) (genpoint-coords v))))
200          (let ((t 
201                 (let recur ((cells cells) (points '()))
202                   (if (null? cells)
203                       (list->kd-tree
204                        points
205                        make-value: make-value
206                        make-point: make-point)
207                       (let ((sections (car cells)))
208                         (let inner ((sections sections) (points points))
209                           (if (null? sections)
210                               (recur (cdr cells) points)
211                               (let ((section (car sections)))
212                                 (inner (cdr sections) 
213                                        (append (cdr (cadr section)) points))
214                                 ))
215                           ))
216                       ))
217                 ))
218            t))
219
220
221        (define (cells-origins->kd-tree cells)
222          (let ((t 
223                 (let recur ((cells cells) (points '()))
224                   (if (null? cells)
225                       (list->kd-tree
226                        points
227                        make-value: (lambda (i v) 
228                                      (list (genpoint-parent-index v)
229                                            (genpoint-parent-distance v)))
230                        make-point: (lambda (v) (genpoint-coords v))
231                        )
232                       (let ((cell (car cells)))
233                         (recur (cdr cells) 
234                                (cons (make-genpoint (cell-origin cell) (cell-index cell) 0. 0)
235                                      points)))
236                       ))
237                 ))
238            t))
239
240
241
242        (define (make-line-segment p dx dy dz) 
243          (let ((c (pcurve:line-segment 3 (list dx dy dz))))
244            (pcurve:translate-curve (point->list p)
245             (pcurve:line-segment 3 (list dx dy dz)))
246            ))
247       
248       
249        ;; auxiliary function to create segment points
250        (define (make-segment si np sp xyz)
251          (list-tabulate 
252           np
253           (lambda (i) (let* ((pi (+ sp i))
254                              (p (make-point 
255                                  (f64vector-ref (car xyz) pi)
256                                  (f64vector-ref (cadr xyz) pi)
257                                  (f64vector-ref (caddr xyz) pi))))
258                         (list si p)
259                         ))
260           ))
261
262
263         
264
265        ;;
266        ;; Creates a process of the given segments and number of points per
267        ;; segment from the given curve.
268        ;;
269        (define (make-segmented-process c ns np)
270          (let ((xyz (pcurve:iterate-curve c (* ns np))))
271            (let recur ((si ns) (ax '()))
272              (if (positive? si)
273                  (let ((sp (* (- si 1) np)))
274                    (recur (- si 1) (append (make-segment si np sp xyz) ax)))
275                  ax)
276              ))
277          )
278
279        ;;
280        ;; Non-segmented process
281        ;;
282        (define (make-process c np)
283          (let ((xyz (pcurve:iterate-curve c np)))
284            (list-tabulate 
285             np
286             (lambda (i) 
287               (make-point 
288                (f64vector-ref (car xyz) i)
289                (f64vector-ref (cadr xyz) i)
290                (f64vector-ref (caddr xyz) i)))
291             ))
292          )
293       
294       
295        ;;
296        ;; Creates a named section containing points from the given segmented processes.
297        ;;
298        (define (make-segmented-section gid p label ps)
299          `(,label . 
300                   ,(fold (lambda (i+proc ax)
301                            (let ((seci (car i+proc)) 
302                                  (proc (cadr i+proc)))
303                              (cons
304                               `(,seci . 
305                                       ,(map (lambda (sp)
306                                               (let ((segi (car sp))
307                                                     (dpoint (cadr sp)))
308                                                 (let ((soma-distance (sqrt (dist2 p dpoint))))
309                                                   (make-genpoint dpoint gid soma-distance segi))
310                                                 ))
311                                             proc))
312                               ax)
313                              ))
314                          '() ps)
315                   ))
316       
317        ;;
318        ;; Non-segmented named section
319        ;;
320        (define (make-section gid p label ps)
321          `(,label . 
322                   ,(fold (lambda (i+proc ax)
323                            (let* ((seci (car i+proc)) 
324                                   (proc (cadr i+proc))
325                                   (pts (map (lambda (dpoint)
326                                               (let ((soma-distance (sqrt (dist2 p dpoint))))
327                                                 (make-genpoint dpoint gid soma-distance #f)))
328                                             proc)))
329                              (d "make-section: gid = ~A pts = ~A~%" gid proc)
330                              (cons `(,seci . ,pts) ax)
331                              ))
332                          '() ps)
333                   ))
334       
335       
336       
337        (define (make-gen-kd-tree data #!key (threshold 0.0))
338         
339          (define (update-bb p x-min y-min z-min x-max y-max z-max)
340            (let ((x (coord 0 p)) (y (coord 1 p)) (z (coord 2 p)))
341              (if (< x (x-min)) (x-min x))
342              (if (< y (y-min)) (y-min y))
343              (if (< z (z-min)) (z-min z))
344              (if (> x (x-max)) (x-max x))
345              (if (> y (y-max)) (y-max y))
346              (if (> z (z-max)) (z-max z))
347              ))
348
349
350          (let* (
351                 (t (list->kd-tree
352                     (filter (lambda (x) (<= threshold (cdr x))) data)
353                     make-value: (lambda (i v) (cdr v))
354                     make-point: (lambda (v) (apply make-point (car v)))
355                     offset: 0
356                     ))
357                 (n (kd-tree-size t))
358                 (x-min (make-parameter +inf.0))
359                 (y-min (make-parameter +inf.0))
360                 (z-min (make-parameter +inf.0))
361                 (x-max (make-parameter -inf.0))
362                 (y-max (make-parameter -inf.0))
363                 (z-max (make-parameter -inf.0))
364                 (bb (begin (kd-tree-for-each
365                             (lambda (p) (update-bb p x-min y-min z-min
366                                                    x-max y-max z-max))
367                             t)
368                            (list (x-min) (y-min) (z-min) (x-max) (y-max) (z-max))))
369                 )
370
371            (cons bb t)
372
373            ))
374
375
376
377
378
379        (define (genpoint-projection prefix my-comm myrank size cells fibers zone cell-start fiber-start)
380
381          (d "rank ~A: zone = ~A~%" myrank zone)
382
383          (fold (lambda (cell ax)
384
385                  (let* ((i (+ cell-start (car cell)))
386                         (root (modulo i size))
387                         (sections (cadr cell)))
388                   
389                    (fold 
390                     
391                     (lambda (sec ax)
392                       
393                       (let ((seci (car sec))
394                             (gxs  (cdr sec)))
395                         
396                         (let ((query-data
397                                (fold 
398                                 (lambda (gd ax)
399                                   (d "rank ~A: querying point ~A (~A)~%" 
400                                      myrank (genpoint-coords gd) 
401                                      (genpoint-parent-index gd))
402                                   (fold
403                                    (lambda (x ax) 
404                                      (let (
405                                            (source (car x))
406                                            (target i)
407                                            (distance (cadr x))
408                                            (segi (genpoint-segment gd))
409                                            )
410                                        (append (list source target distance segi seci) ax)
411                                        ))
412                                    ax
413                                    (delete-duplicates
414                                     (map (lambda (x) 
415                                            (d "rank ~A: query result = ~A (~A) (~A) ~%" 
416                                               myrank (kdnn-point x) (kdnn-distance x) (kdnn-parent-index x))
417                                            (list (+ fiber-start (kdnn-parent-index x))
418                                                  (+ (kdnn-distance x) (genpoint-parent-distance gd)  (kdnn-parent-distance x))
419                                                  ))
420                                          (kd-tree-near-neighbors* fibers zone (genpoint-coords gd)))
421                                     (lambda (u v)
422                                       (= (car u) (car v)))
423                                     )
424                                    ))
425                                 '() gxs)))
426                           
427                           (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
428
429                                  (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
430                             
431                             (append res1 ax))
432                           
433                           ))
434                       )
435                     ax sections)
436                    ))
437                '() cells ))
438
439
440       
441
442        (define (point-projection prefix my-comm myrank size pts fibers zone point-start nn-filter)
443          (fold (lambda (px ax)
444
445                  (d "~A: rank ~A: px = ~A~%"  prefix myrank px)
446
447                  (let* ((i (+ point-start (car px)))
448                         (root (modulo i size))
449                         (dd (d "~A: rank ~A: querying point ~A (root ~A)~%" prefix myrank px root))
450                         (query-data
451                          (fold 
452                           (lambda (pd ax)
453                             (fold
454                              (lambda (x ax) 
455                                (let ((source (car x))
456                                      (target i)
457                                      (distance (cadr x)))
458                                  (if (and (> distance  0.) (not (= source target)))
459                                      (append (list source target distance) ax)
460                                      ax)
461                                  ))
462                              ax
463                              (delete-duplicates
464                               (map (lambda (x) 
465                                      (let ((res (list (car (cadar x))  (+ (caddr x) (cadr (cadar x))))))
466                                        (d "~A: axon: x = ~A res = ~A~%" prefix x res)
467                                        res))
468                                    (nn-filter pd (kd-tree-near-neighbors* fibers zone pd))
469                                    )
470                               (lambda (u v) (= (car u) (car v)))
471                               )
472                              ))
473                           '() (list (cadr px))))
474                         )
475
476                    (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
477                           (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
478                      (append res1 ax))
479                    ))
480
481                '() pts))
482       
483
484
485
486
487        (define-datatype layer-boundary layer-boundary?
488          (Bounds (top number?) (left number?) (bottom number?) (right number?))
489          (BoundsXZ (top number?) (left number?) (bottom number?) (right number?)
490                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
491          (BoundsYZ (top number?) (left number?) (bottom number?) (right number?)
492                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
493          )
494
495
496
497        (define (boundary-z-extent-function boundary)
498          (cases layer-boundary boundary
499                 (Bounds (top left bottom right) 
500                         (lambda (x y) 0.))
501                 (BoundsXZ (top left bottom right n k x y d d2) 
502                           (lambda (xp yp) 
503                             (let-values (((y0tab y1tab y2tab res)
504                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
505                               (f64vector-ref y0tab 0))))
506                 (BoundsYZ (top left bottom right n k x y d d2) 
507                           (lambda (xp yp) 
508                             (let-values (((y0tab y1tab y2tab res)
509                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
510                               (f64vector-ref y0tab 0))))
511                 ))
512
513
514        (define (point2d-rejection top left bottom right)
515            (lambda (p)
516              (let ((x (coord 0 p)) (y (coord 1 p)))
517                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
518            )
519
520
521        (define (compute-point3d p zu z-extent-function)
522          (let* ((x (coord 0 p))
523                 (y (coord 1 p))
524                 (z-extent (z-extent-function x y)))
525            (if (zero? zu)
526                (make-point x y zu)
527                (make-point x y (fp* zu z-extent)))
528            ))
529
530
531        (define (cluster-point-rejection p cluster-pts mean-distance)
532          (let ((D (* 4 mean-distance mean-distance))
533                (nn (kd-tree-nearest-neighbor cluster-pts p)))
534            (and (< (dist2 p nn) D) p)))
535
536
537
538        (define (XZAxis n k x-points z-points boundary)
539         
540          (if (not (>= n 3)) 
541              (error 'generate-boundary "XZAxis requires at least 3 interpolation points (n >= 3)"))
542         
543          (cases layer-boundary boundary
544                 (Bounds (top left bottom right) 
545
546                         (let-values (((d d2 constr errc diagn)
547                                       (bvsp-spline:compute n k x-points z-points)))
548                           
549                           (if (not (zero? errc)) 
550                               (error 'generate-boundary "error in constructing spline from boundary points" errc))
551                           
552                           (BoundsXZ top left bottom right n k x-points z-points d d2)))
553                 
554                 (else (error 'generate-boundary "boundary argument to XZAxis is already a pseudo-3D boundary")))
555          )
556
557
558        (define (Grid x-spacing y-spacing z-spacing boundary)
559
560          (match-let (((top left bottom right)
561                       (cases layer-boundary boundary
562                                   (Bounds (top left bottom right) 
563                                           (list top left bottom right))
564                                   (BoundsXZ (top left bottom right n k x y d d2) 
565                                             (list top left bottom right))
566                                   (BoundsYZ (top left bottom right n k x y d d2) 
567                                             (list top left bottom right))
568                                   )))
569
570          (let* (
571                 (x-extent   (- right left))
572                 (y-extent   (- top bottom))
573                 (z-extent-function
574                  (boundary-z-extent-function boundary))
575                 (compute-grid-points3d
576                  (lambda (p z-spacing z-extent-function)
577                    (let* ((x (coord 0 p))
578                           (y (coord 1 p))
579                           (z-extent (z-extent-function x y)))
580                      (let recur ((points '()) (z (/ z-spacing 2.)))
581                        (if (> z z-extent)
582                            points
583                            (recur (cons (make-point x y z) points) (+ z z-spacing)))
584                        ))
585                    ))
586                 )
587           
588            (d "Grid: x-spacing = ~A~%" x-spacing)
589            (d "Grid: y-spacing = ~A~%" y-spacing)
590            (d "Grid: z-spacing = ~A~%" z-spacing)
591           
592            (d "Grid: x-extent = ~A~%" x-extent)
593            (d "Grid: y-extent = ~A~%" y-extent)
594           
595            (let ((x-points (let recur ((points '()) (x (/ x-spacing 2.)))
596                              (if (> x x-extent)
597                                  (list->f64vector (reverse points))
598                                  (recur (cons x points) (+ x x-spacing)))))
599                  (y-points (let recur ((points '()) (y (/ y-spacing 2.)))
600                              (if (> y y-extent)
601                                  (list->f64vector (reverse points))
602                                  (recur (cons y points) (+ y y-spacing)))))
603                  )
604             
605              (let ((nx (f64vector-length x-points))
606                    (ny (f64vector-length y-points))
607                    )
608               
609                (let recur ((i 0) (j 0) (ax '()))
610                  (if (< i nx)
611                      (let ((x (f64vector-ref x-points i)))
612                        (if (< j ny)
613                            (let* ((y   (f64vector-ref y-points j))
614                                   (p   (make-point x y))
615                                   (p3ds (if (zero? z-spacing)
616                                             (list (make-point x y 0.0))
617                                             (compute-grid-points3d p z-spacing z-extent-function)))
618                                   )
619                              (recur i (+ 1 j) (if p3ds (append p3ds ax) ax)))
620                            (recur (+ 1 i) 0 ax)))
621                      (let* ((t (list->kd-tree ax))
622                             (n (kd-tree-size t)))
623                        (list t boundary)
624                        ))
625                  )))
626            ))
627          )
628
629        (define (UniformRandomPointProcess n x-seed y-seed boundary)
630
631          (match-let (((top left bottom right)
632                       (cases layer-boundary boundary
633                                   (Bounds (top left bottom right) 
634                                           (list top left bottom right))
635                                   (BoundsXZ (top left bottom right n k x y d d2) 
636                                             (list top left bottom right))
637                                   (BoundsYZ (top left bottom right n k x y d d2) 
638                                             (list top left bottom right))
639                                   )))
640
641          (let* (
642                 (x-extent   (- right left))
643                 (y-extent   (- top bottom))
644                 (z-extent-function (boundary-z-extent-function boundary))
645                 )
646
647            (let ((x-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init x-seed)))
648                  (y-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init y-seed)))
649                  (z-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init (current-milliseconds)))))
650             
651              (let ((point-rejection1 (point2d-rejection top left bottom right)))
652               
653                (let recur ((i 0) (ax '()))
654                  (if (< i n)
655                      (let ((x (f64vector-ref x-points i))
656                            (y (f64vector-ref y-points i))
657                            (z (f64vector-ref z-points i)))
658                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
659                          (let ((p3d 
660                                 (and (point-rejection1 p)
661                                      (compute-point3d p z z-extent-function))))
662                            (recur (+ 1 i) (if p3d (cons p3d ax) ax)))))
663                      (let* ((t (list->kd-tree ax))
664                             (n (kd-tree-size t)))
665
666                        (list t boundary))))
667                ))
668            ))
669          )
670
671
672        (define (ClusteredRandomPointProcess cluster-pts n mean-distance x-seed y-seed boundary)
673
674          (match-let (((top left bottom right)
675                       (cases layer-boundary boundary
676                                   (Bounds (top left bottom right) 
677                                           (list top left bottom right))
678                                   (BoundsXZ (top left bottom right n k x y d d2) 
679                                             (list top left bottom right))
680                                   (BoundsYZ (top left bottom right n k x y d d2) 
681                                             (list top left bottom right))
682                                   )))
683
684
685          (let* (
686                 (x-extent   (- right left))
687                 (y-extent   (- top bottom))
688                 (z-extent-function (boundary-z-extent-function boundary))
689                 )
690
691            (let recur ((pts '()) (x-seed x-seed) (y-seed y-seed))
692
693              (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
694                    (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
695                    (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
696               
697                (let ((point-rejection1 (point2d-rejection top left bottom right)))
698                 
699                  (let inner-recur ((j 0) (ax pts))
700                    (if (< j n)
701                        (let ((x (f64vector-ref x-points j))
702                              (y (f64vector-ref y-points j))
703                              (z (f64vector-ref z-points j)))
704                          (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
705                            (let ((p3d 
706                                   (and (point-rejection1 p)
707                                        (compute-point3d p z z-extent-function))))
708                              (let ((pp (cluster-point-rejection p3d cluster-pts mean-distance)))
709                                (inner-recur (+ 1 j)  (if pp (cons pp ax) ax))))))
710
711                        (if (< (length ax) n)
712                            (recur ax (+ 1 x-seed) (+ 1 y-seed))
713                            (let* ((t (list->kd-tree (take ax n)))
714                                   (n (kd-tree-size t)))
715                             
716                              (list t boundary))))
717                    ))
718                ))
719            ))
720          )
721
722
723
724       
725        (define comment-pat (string->irregex "^#.*"))
726
727
728        (define (load-points-from-file filename)
729
730          (let ((in (open-input-file filename)))
731           
732            (if (not in) (error 'load-points-from-file "file not found" filename))
733
734            (let* ((lines
735                    (filter (lambda (line) (not (irregex-match comment-pat line)))
736                            (read-lines in)))
737
738                   (point-data
739                    (filter-map
740                     (lambda (line) 
741                       (let ((lst (map string->number (string-split line " \t"))))
742                         (and (not (null? lst)) (apply make-point lst))))
743                     lines))
744
745                   (point-tree (list->kd-tree point-data))
746                   )
747             
748              (list point-tree)
749             
750              ))
751          )
752
753
754        (define (segment-projection label source-tree target-sections zone my-comm myrank size)
755
756          (MPI:barrier my-comm)
757         
758          (let ((my-results
759                 (genpoint-projection label my-comm myrank size target-sections source-tree zone 0 0)))
760
761            (MPI:barrier my-comm)
762
763            (print "my-results = " my-results)
764     
765            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
766              (lambda (out-sources)
767                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
768                  (lambda (out-targets)
769                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
770                      (lambda (out-distances)
771                        (call-with-output-file (sprintf "~Asegments~A.dat"  label (if (> size 1) myrank ""))
772                          (lambda (out-segments)
773                            (for-each
774                             (lambda (my-data)
775                               (let* ((my-entry-len 5)
776                                      (my-data-len (/ (f64vector-length my-data) my-entry-len)))
777                                 (d "rank ~A: length my-data = ~A~%" myrank my-data-len)
778                                 (let recur ((m 0))
779                                   (if (< m my-data-len)
780                                       (let* (
781                                              (my-entry-offset (* m my-entry-len))
782                                              (source (inexact->exact (f64vector-ref my-data my-entry-offset)))
783                                              (target (inexact->exact (f64vector-ref my-data (+ 1 my-entry-offset))))
784                                              (distance (f64vector-ref my-data (+ 2 my-entry-offset)))
785                                              (segment (f64vector-ref my-data (+ 3 my-entry-offset)))
786                                              (section (f64vector-ref my-data (+ 4 my-entry-offset)))
787                                              )
788                                         (fprintf out-sources   "~A~%" source)
789                                         (fprintf out-targets   "~A~%" target)
790                                         (fprintf out-distances "~A~%" distance)
791                                         (fprintf out-segments  "~A ~A~%" segment section)
792                                         (recur (+ 1 m)))))
793                                 ))
794                             my-results))
795                          ))
796                      ))
797                  ))
798              ))
799          )
800       
801
802        (define (projection label source-tree target zone my-comm myrank size) 
803
804          (MPI:barrier my-comm)
805         
806          (let ((my-results (point-projection label my-comm myrank size target source-tree zone 0 (lambda (x nn) nn))))
807
808          (MPI:barrier my-comm)
809           
810            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
811              (lambda (out-sources)
812                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
813                  (lambda (out-targets)
814                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
815                      (lambda (out-distances)
816                        (for-each
817                         (lambda (my-data)
818                           (let* ((my-entry-len 3)
819                                  (my-data-len (/ (f64vector-length my-data) my-entry-len)))
820                             (d "~A: rank ~A: length my-data = ~A~%" label myrank my-data-len)
821                             (let recur ((m 0))
822                               (if (< m my-data-len)
823                                   (let ((source (inexact->exact (f64vector-ref my-data (* m my-entry-len))))
824                                         (target (inexact->exact (f64vector-ref my-data (+ 1 (* m my-entry-len)))))
825                                         (distance (f64vector-ref my-data (+ 2 (* m my-entry-len)))))
826                                     (fprintf out-sources "~A~%" source)
827                                     (fprintf out-targets "~A~%" target)
828                                     (fprintf out-distances "~A~%" distance)
829                                     (recur (+ 1 m)))))
830                             ))
831                         my-results)
832                        ))
833                    ))
834                ))
835            ))
836       
837
838        (include "calc-parser.scm")
839
840       
841        (define (load-config-file filename)
842          (let ((in (open-input-file filename)))
843            (if (not in) (error 'load-config-file "file not found" filename))
844            (init-bindings)
845            (let* ((lines (reverse (filter (lambda (x) (not (string-null? x))) (read-lines in))))
846                   (properties (filter-map
847                                (lambda (line) 
848                                  (and (not (string-prefix? "//" line))
849                                       (let ((lst (string-split line "=")))
850                                         (and (> (length lst) 1)
851                                              (let ((key (string->symbol (string-trim-both (car lst) #\space)))
852                                                    (val (calc-eval-string (cadr lst))))
853                                                (add-binding key val)
854                                                (cons key val))))))
855                                lines)))
856              properties
857              ))
858          )
859
860       
861)
862
863       
Note: See TracBrowser for help on using the repository browser.