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

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

picnic: some more bug fixes in point/segment projections and in GLH example model

File size: 37.5 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 fold-right 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 (sprintf "~A.pointset.dat" name)
137              (lambda (out)
138                (for-each (match-lambda
139                           ((gid p)
140                            (fprintf out "~A ~A ~A~%" 
141                                     (coord 0 p)
142                                     (coord 1 p)
143                                     (coord 2 p))))
144                          pts))
145              ))
146
147       
148        (define (write-layout name pts #!optional rank)
149            (call-with-output-file (if rank 
150                                       (sprintf "~A.~A.layout.dat" name rank)
151                                       (sprintf "~A.layout.dat" name))
152              (lambda (out)
153                (for-each (match-lambda
154                           ((gid p)
155                            (fprintf out "~A ~A ~A ~A~%" 
156                                     gid
157                                     (coord 0 p)
158                                     (coord 1 p)
159                                     (coord 2 p))))
160                          pts))
161              ))
162
163        (define (write-sections forest-name section-name layout sections #!optional rank)
164          (call-with-output-file (if rank
165                                     (sprintf "~A.~A.~A.section.dat" forest-name section-name rank)
166                                     (sprintf "~A.~A.section.dat" forest-name section-name))
167            (lambda (out)
168              (for-each
169               (match-lambda*
170                (((gid p) section)
171                 (fprintf out "~A " gid)
172                 (for-each
173                  (lambda (neurites)
174                    (for-each
175                     (lambda (gd)
176                       (let ((p (genpoint-coords gd)))
177                         (fprintf out "~A ~A ~A "
178                                  (coord 0 p)
179                                  (coord 1 p)
180                                  (coord 2 p))))
181                     (cdr neurites)))
182                  (cdr section))
183                 (fprintf out "~%")
184                 ))
185               layout 
186               sections))))
187
188
189        (define (cells-sections->kd-tree cells section-name #!key
190                                         (make-value
191                                          (lambda (i v) 
192                                            (list (genpoint-parent-index v)
193                                                  (genpoint-parent-distance v))))
194                                         (make-point
195                                          (lambda (v) (genpoint-coords v))))
196          (let ((t 
197                 (let recur ((cells cells) (points '()))
198                   (if (null? cells)
199                       (list->kd-tree
200                        points
201                        make-value: make-value
202                        make-point: make-point)
203                       (let ((cell (car cells)))
204                         (recur (cdr cells) 
205                                (let inner ((sections (append (cell-section-ref section-name cell)))
206                                            (points points))
207                                  (if (null? sections) points
208                                      (inner
209                                       (cdr sections)
210                                       (append (cdr (car sections)) points))
211                                      ))
212                                ))
213                       ))
214                 ))
215            t))
216
217        (define (sections->kd-tree cells #!key
218                                   (make-value
219                                    (lambda (i v) 
220                                      (list (genpoint-parent-index v)
221                                            (genpoint-parent-distance v))))
222                                   (make-point
223                                    (lambda (v) (genpoint-coords v))))
224          (let ((t 
225                 (let recur ((cells cells) (points '()))
226                   (if (null? cells)
227                       (list->kd-tree
228                        points
229                        make-value: make-value
230                        make-point: make-point)
231                       (let ((sections (car cells)))
232                         (let inner ((sections sections) (points points))
233                           (if (null? sections)
234                               (recur (cdr cells) points)
235                               (let ((section (car sections)))
236                                 (inner (cdr sections) 
237                                        (append (cdr (cadr section)) points))
238                                 ))
239                           ))
240                       ))
241                 ))
242            t))
243
244
245        (define (cells-origins->kd-tree cells)
246          (let ((t 
247                 (let recur ((cells cells) (points '()))
248                   (if (null? cells)
249                       (list->kd-tree
250                        points
251                        make-value: (lambda (i v) 
252                                      (list (genpoint-parent-index v)
253                                            (genpoint-parent-distance v)))
254                        make-point: (lambda (v) (genpoint-coords v))
255                        )
256                       (let ((cell (car cells)))
257                         (recur (cdr cells) 
258                                (cons (make-genpoint (cell-origin cell) (cell-index cell) 0. 0)
259                                      points)))
260                       ))
261                 ))
262            t))
263
264
265
266        (define (make-line-segment p dx dy dz) 
267          (let ((c (pcurve:line-segment 3 (list dx dy dz))))
268            (pcurve:translate-curve (point->list p)
269             (pcurve:line-segment 3 (list dx dy dz)))
270            ))
271       
272       
273        ;; auxiliary function to create segment points
274        (define (make-segment si np sp xyz)
275          (list-tabulate 
276           np
277           (lambda (i) (let* ((pi (+ sp i))
278                              (p (make-point 
279                                  (f64vector-ref (car xyz) pi)
280                                  (f64vector-ref (cadr xyz) pi)
281                                  (f64vector-ref (caddr xyz) pi))))
282                         (list si p)
283                         ))
284           ))
285
286
287         
288
289        ;;
290        ;; Creates a process of the given segments and number of points per
291        ;; segment from the given curve.
292        ;;
293        (define (make-segmented-process c ns np)
294          (let ((xyz (pcurve:iterate-curve c (* ns np))))
295            (let recur ((si ns) (ax '()))
296              (if (positive? si)
297                  (let ((sp (* (- si 1) np)))
298                    (recur (- si 1) (append (make-segment si np sp xyz) ax)))
299                  ax)
300              ))
301          )
302
303        ;;
304        ;; Non-segmented process
305        ;;
306        (define (make-process c np)
307          (let ((xyz (pcurve:iterate-curve c np)))
308            (list-tabulate 
309             np
310             (lambda (i) 
311               (make-point 
312                (f64vector-ref (car xyz) i)
313                (f64vector-ref (cadr xyz) i)
314                (f64vector-ref (caddr xyz) i)))
315             ))
316          )
317       
318       
319        ;;
320        ;; Creates a named section containing points from the given segmented processes.
321        ;;
322        (define (make-segmented-section gid p label ps)
323          `(,label . 
324                   ,(fold (lambda (i+proc ax)
325                            (let ((seci (car i+proc)) 
326                                  (proc (cadr i+proc)))
327                              (cons
328                               `(,seci . 
329                                       ,(map (lambda (sp)
330                                               (let ((segi (car sp))
331                                                     (dpoint (cadr sp)))
332                                                 (let ((soma-distance (sqrt (dist2 p dpoint))))
333                                                   (make-genpoint dpoint gid soma-distance segi))
334                                                 ))
335                                             proc))
336                               ax)
337                              ))
338                          '() ps)
339                   ))
340       
341        ;;
342        ;; Non-segmented named section
343        ;;
344        (define (make-section gid p label ps)
345          `(,label . 
346                   ,(fold (lambda (i+proc ax)
347                            (let* ((seci (car i+proc)) 
348                                   (proc (cadr i+proc))
349                                   (pts (map (lambda (dpoint)
350                                               (let ((soma-distance (sqrt (dist2 p dpoint))))
351                                                 (make-genpoint dpoint gid soma-distance #f)))
352                                             proc)))
353                              (d "make-section: label = ~A gid = ~A p = ~A pts = ~A~%" 
354                                 label gid p proc)
355                              (cons `(,seci . ,pts) ax)
356                              ))
357                          '() ps)
358                   ))
359       
360       
361       
362        (define (make-gen-kd-tree data #!key (threshold 0.0))
363         
364          (define (update-bb p x-min y-min z-min x-max y-max z-max)
365            (let ((x (coord 0 p)) (y (coord 1 p)) (z (coord 2 p)))
366              (if (< x (x-min)) (x-min x))
367              (if (< y (y-min)) (y-min y))
368              (if (< z (z-min)) (z-min z))
369              (if (> x (x-max)) (x-max x))
370              (if (> y (y-max)) (y-max y))
371              (if (> z (z-max)) (z-max z))
372              ))
373
374
375          (let* (
376                 (t (list->kd-tree
377                     (filter (lambda (x) (<= threshold (cdr x))) data)
378                     make-value: (lambda (i v) (cdr v))
379                     make-point: (lambda (v) (apply make-point (car v)))
380                     offset: 0
381                     ))
382                 (n (kd-tree-size t))
383                 (x-min (make-parameter +inf.0))
384                 (y-min (make-parameter +inf.0))
385                 (z-min (make-parameter +inf.0))
386                 (x-max (make-parameter -inf.0))
387                 (y-max (make-parameter -inf.0))
388                 (z-max (make-parameter -inf.0))
389                 (bb (begin (kd-tree-for-each
390                             (lambda (p) (update-bb p x-min y-min z-min
391                                                    x-max y-max z-max))
392                             t)
393                            (list (x-min) (y-min) (z-min) (x-max) (y-max) (z-max))))
394                 )
395
396            (cons bb t)
397
398            ))
399
400
401
402
403
404        (define (genpoint-projection prefix my-comm myrank size cells fibers zone cell-start fiber-start)
405
406          (d "rank ~A: zone = ~A~%" myrank zone)
407
408          (fold (lambda (cell ax)
409
410                  (let* ((i (+ cell-start (car cell)))
411                         (root (modulo i size))
412                         (sections (cadr cell)))
413                   
414                    (fold 
415                     
416                     (lambda (sec ax)
417                       
418                       (let ((seci (car sec))
419                             (gxs  (cdr sec)))
420                         
421                         (let ((query-data
422                                (fold 
423                                 (lambda (gd ax)
424                                   (d "rank ~A: querying point ~A (~A)~%" 
425                                      myrank (genpoint-coords gd) 
426                                      (genpoint-parent-index gd))
427                                   (fold
428                                    (lambda (x ax) 
429                                      (let (
430                                            (source (car x))
431                                            (target i)
432                                            (distance (cadr x))
433                                            (segi (genpoint-segment gd))
434                                            )
435                                        (if (not segi)
436                                            (error 'genpoint-projection "point does not have segment index" gd))
437                                        (append (list source target distance segi seci) ax)
438                                        ))
439                                    ax
440                                    (delete-duplicates
441                                     (map (lambda (x) 
442                                            (d "rank ~A: query result = ~A (~A) (~A) ~%" 
443                                               myrank (kdnn-point x) (kdnn-distance x) (kdnn-parent-index x))
444                                            (list (+ fiber-start (kdnn-parent-index x))
445                                                  (+ (kdnn-distance x) (genpoint-parent-distance gd)  (kdnn-parent-distance x))
446                                                  ))
447                                          (kd-tree-near-neighbors* fibers zone (genpoint-coords gd)))
448                                     (lambda (u v)
449                                       (= (car u) (car v)))
450                                     )
451                                    ))
452                                 '() gxs)))
453                           (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
454
455                                  (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
456                             
457                             (append res1 ax))
458                           
459                           ))
460                       )
461                     ax sections)
462                    ))
463                '() cells ))
464
465
466       
467
468        (define (point-projection prefix my-comm myrank size pts fibers zone point-start nn-filter)
469          (fold (lambda (px ax)
470
471                  (d "~A: rank ~A: px = ~A~%"  prefix myrank px)
472
473                  (let* ((i (+ point-start (car px)))
474                         (root (modulo i size))
475                         (dd (d "~A: rank ~A: querying point ~A (root ~A)~%" prefix myrank px root))
476                         (query-data
477                          (fold 
478                           (lambda (pd ax)
479                             (fold
480                              (lambda (x ax) 
481                                (let ((source (car x))
482                                      (target i)
483                                      (distance (cadr x)))
484                                  (if (and (> distance  0.) (not (= source target)))
485                                      (append (list source target distance) ax)
486                                      ax)
487                                  ))
488                              ax
489                              (delete-duplicates
490                               (map (lambda (x) 
491                                      (let ((res (list (car (cadar x)) 
492                                                       (+ (caddr x) (cadr (cadar x))))))
493                                        (d "~A: x = ~A res = ~A~%" prefix x res)
494                                        res))
495                                    (nn-filter pd (kd-tree-near-neighbors* fibers zone pd))
496                                    )
497                               (lambda (u v) (= (car u) (car v)))
498                               )
499                              ))
500                           '() (cadr px)))
501                         )
502
503                    (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
504                           (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
505                      (append res1 ax))
506                    ))
507
508                '() pts))
509       
510
511
512
513
514        (define-datatype layer-boundary layer-boundary?
515          (Bounds (top number?) (left number?) (bottom number?) (right number?))
516          (BoundsXZ (top number?) (left number?) (bottom number?) (right number?)
517                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
518          (BoundsYZ (top number?) (left number?) (bottom number?) (right number?)
519                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
520          )
521
522
523
524        (define (boundary-z-extent-function boundary)
525          (cases layer-boundary boundary
526                 (Bounds (top left bottom right) 
527                         (lambda (x y) 0.))
528                 (BoundsXZ (top left bottom right n k x y d d2) 
529                           (lambda (xp yp) 
530                             (let-values (((y0tab y1tab y2tab res)
531                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
532                               (f64vector-ref y0tab 0))))
533                 (BoundsYZ (top left bottom right n k x y d d2) 
534                           (lambda (xp yp) 
535                             (let-values (((y0tab y1tab y2tab res)
536                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
537                               (f64vector-ref y0tab 0))))
538                 ))
539
540
541        (define (point2d-rejection top left bottom right)
542            (lambda (p)
543              (let ((x (coord 0 p)) (y (coord 1 p)))
544                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
545            )
546
547
548        (define (compute-point3d p zu z-extent-function)
549          (let* ((x (coord 0 p))
550                 (y (coord 1 p))
551                 (z-extent (z-extent-function x y)))
552            (if (zero? zu)
553                (make-point x y zu)
554                (make-point x y (fp* zu z-extent)))
555            ))
556
557
558        (define (cluster-point-rejection p cluster-pts mean-distance)
559          (let ((D (* 4 mean-distance mean-distance))
560                (nn (kd-tree-nearest-neighbor cluster-pts p)))
561            (and (< (dist2 p nn) D) p)))
562
563
564
565        (define (XZAxis n k x-points z-points boundary)
566         
567          (if (not (>= n 3)) 
568              (error 'generate-boundary "XZAxis requires at least 3 interpolation points (n >= 3)"))
569         
570          (cases layer-boundary boundary
571                 (Bounds (top left bottom right) 
572
573                         (let-values (((d d2 constr errc diagn)
574                                       (bvsp-spline:compute n k x-points z-points)))
575                           
576                           (if (not (zero? errc)) 
577                               (error 'generate-boundary "error in constructing spline from boundary points" errc))
578                           
579                           (BoundsXZ top left bottom right n k x-points z-points d d2)))
580                 
581                 (else (error 'generate-boundary "boundary argument to XZAxis is already a pseudo-3D boundary")))
582          )
583
584
585        (define (Grid x-spacing y-spacing z-spacing boundary)
586
587          (match-let (((top left bottom right)
588                       (cases layer-boundary boundary
589                                   (Bounds (top left bottom right) 
590                                           (list top left bottom right))
591                                   (BoundsXZ (top left bottom right n k x y d d2) 
592                                             (list top left bottom right))
593                                   (BoundsYZ (top left bottom right n k x y d d2) 
594                                             (list top left bottom right))
595                                   )))
596
597          (let* (
598                 (x-extent   (- right left))
599                 (y-extent   (- top bottom))
600                 (z-extent-function
601                  (boundary-z-extent-function boundary))
602                 (compute-grid-points3d
603                  (lambda (p z-spacing z-extent-function)
604                    (let* ((x (coord 0 p))
605                           (y (coord 1 p))
606                           (z-extent (z-extent-function x y)))
607                      (let recur ((points '()) (z (/ z-spacing 2.)))
608                        (if (> z z-extent)
609                            points
610                            (recur (cons (make-point x y z) points) (+ z z-spacing)))
611                        ))
612                    ))
613                 )
614           
615            (d "Grid: x-spacing = ~A~%" x-spacing)
616            (d "Grid: y-spacing = ~A~%" y-spacing)
617            (d "Grid: z-spacing = ~A~%" z-spacing)
618           
619            (d "Grid: x-extent = ~A~%" x-extent)
620            (d "Grid: y-extent = ~A~%" y-extent)
621           
622            (let ((x-points (let recur ((points '()) (x (/ x-spacing 2.)))
623                              (if (> x x-extent)
624                                  (list->f64vector (reverse points))
625                                  (recur (cons x points) (+ x x-spacing)))))
626                  (y-points (let recur ((points '()) (y (/ y-spacing 2.)))
627                              (if (> y y-extent)
628                                  (list->f64vector (reverse points))
629                                  (recur (cons y points) (+ y y-spacing)))))
630                  )
631             
632              (let ((nx (f64vector-length x-points))
633                    (ny (f64vector-length y-points))
634                    )
635               
636                (let recur ((i 0) (j 0) (ax '()))
637                  (if (< i nx)
638                      (let ((x (f64vector-ref x-points i)))
639                        (if (< j ny)
640                            (let* ((y   (f64vector-ref y-points j))
641                                   (p   (make-point x y))
642                                   (p3ds (if (zero? z-spacing)
643                                             (list (make-point x y 0.0))
644                                             (compute-grid-points3d p z-spacing z-extent-function)))
645                                   )
646                              (recur i (+ 1 j) (if p3ds (append p3ds ax) ax)))
647                            (recur (+ 1 i) 0 ax)))
648                      (let* ((t (list->kd-tree ax))
649                             (n (kd-tree-size t)))
650                        (list t boundary)
651                        ))
652                  )))
653            ))
654          )
655
656        (define (UniformRandomPointProcess n x-seed y-seed boundary)
657
658          (match-let (((top left bottom right)
659                       (cases layer-boundary boundary
660                                   (Bounds (top left bottom right) 
661                                           (list top left bottom right))
662                                   (BoundsXZ (top left bottom right n k x y d d2) 
663                                             (list top left bottom right))
664                                   (BoundsYZ (top left bottom right n k x y d d2) 
665                                             (list top left bottom right))
666                                   )))
667
668          (let* (
669                 (x-extent   (- right left))
670                 (y-extent   (- top bottom))
671                 (z-extent-function (boundary-z-extent-function boundary))
672                 )
673
674            (let ((x-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init x-seed)))
675                  (y-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init y-seed)))
676                  (z-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init (current-milliseconds)))))
677             
678              (let ((point-rejection1 (point2d-rejection top left bottom right)))
679               
680                (let recur ((i 0) (ax '()))
681                  (if (< i n)
682                      (let ((x (f64vector-ref x-points i))
683                            (y (f64vector-ref y-points i))
684                            (z (f64vector-ref z-points i)))
685                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
686                          (let ((p3d 
687                                 (and (point-rejection1 p)
688                                      (compute-point3d p z z-extent-function))))
689                            (recur (+ 1 i) (if p3d (cons p3d ax) ax)))))
690                      (let* ((t (list->kd-tree ax))
691                             (n (kd-tree-size t)))
692
693                        (list t boundary))))
694                ))
695            ))
696          )
697
698
699        (define (ClusteredRandomPointProcess cluster-pts n mean-distance x-seed y-seed boundary)
700
701          (match-let (((top left bottom right)
702                       (cases layer-boundary boundary
703                                   (Bounds (top left bottom right) 
704                                           (list top left bottom right))
705                                   (BoundsXZ (top left bottom right n k x y d d2) 
706                                             (list top left bottom right))
707                                   (BoundsYZ (top left bottom right n k x y d d2) 
708                                             (list top left bottom right))
709                                   )))
710
711
712          (let* (
713                 (x-extent   (- right left))
714                 (y-extent   (- top bottom))
715                 (z-extent-function (boundary-z-extent-function boundary))
716                 )
717
718            (let recur ((pts '()) (x-seed x-seed) (y-seed y-seed))
719
720              (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
721                    (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
722                    (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
723               
724                (let ((point-rejection1 (point2d-rejection top left bottom right)))
725                 
726                  (let inner-recur ((j 0) (ax pts))
727                    (if (< j n)
728                        (let ((x (f64vector-ref x-points j))
729                              (y (f64vector-ref y-points j))
730                              (z (f64vector-ref z-points j)))
731                          (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
732                            (let ((p3d 
733                                   (and (point-rejection1 p)
734                                        (compute-point3d p z z-extent-function))))
735                              (let ((pp (cluster-point-rejection p3d cluster-pts mean-distance)))
736                                (inner-recur (+ 1 j)  (if pp (cons pp ax) ax))))))
737
738                        (if (< (length ax) n)
739                            (recur ax (+ 1 x-seed) (+ 1 y-seed))
740                            (let* ((t (list->kd-tree (take ax n)))
741                                   (n (kd-tree-size t)))
742                             
743                              (list t boundary))))
744                    ))
745                ))
746            ))
747          )
748
749
750
751       
752        (define comment-pat (string->irregex "^#.*"))
753
754
755        (define (load-points-from-file filename)
756
757          (let ((in (open-input-file filename)))
758           
759            (if (not in) (error 'load-points-from-file "file not found" filename))
760
761            (let* ((lines
762                    (filter (lambda (line) (not (irregex-match comment-pat line)))
763                            (read-lines in)))
764
765                   (point-data
766                    (filter-map
767                     (lambda (line) 
768                       (let ((lst (map string->number (string-split line " \t"))))
769                         (and (not (null? lst)) (apply make-point lst))))
770                     lines))
771
772                   (point-tree (list->kd-tree point-data))
773                   )
774             
775              (list point-tree)
776             
777              ))
778          )
779
780
781        (define (segment-projection label source-tree target-sections zone my-comm myrank size)
782
783          (MPI:barrier my-comm)
784         
785          (let ((my-results
786                 (genpoint-projection label my-comm myrank size target-sections source-tree zone 0 0)))
787
788            (MPI:barrier my-comm)
789
790            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
791              (lambda (out-sources)
792                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
793                  (lambda (out-targets)
794                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
795                      (lambda (out-distances)
796                        (call-with-output-file (sprintf "~Asegments~A.dat"  label (if (> size 1) myrank ""))
797                          (lambda (out-segments)
798                            (for-each
799                             (lambda (my-data)
800                               (let* ((my-entry-len 5)
801                                      (my-data-len (/ (f64vector-length my-data) my-entry-len)))
802                                 (d "rank ~A: length my-data = ~A~%" myrank my-data-len)
803                                 (let recur ((m 0))
804                                   (if (< m my-data-len)
805                                       (let* (
806                                              (my-entry-offset (* m my-entry-len))
807                                              (source (inexact->exact (f64vector-ref my-data my-entry-offset)))
808                                              (target (inexact->exact (f64vector-ref my-data (+ 1 my-entry-offset))))
809                                              (distance (f64vector-ref my-data (+ 2 my-entry-offset)))
810                                              (segment (f64vector-ref my-data (+ 3 my-entry-offset)))
811                                              (section (f64vector-ref my-data (+ 4 my-entry-offset)))
812                                              )
813                                         (fprintf out-sources   "~A~%" source)
814                                         (fprintf out-targets   "~A~%" target)
815                                         (fprintf out-distances "~A~%" distance)
816                                         (fprintf out-segments  "~A ~A~%" segment section)
817                                         (recur (+ 1 m)))))
818                                 ))
819                             my-results))
820                          ))
821                      ))
822                  ))
823              ))
824          )
825       
826
827        (define (projection label source-tree target zone my-comm myrank size) 
828
829          (MPI:barrier my-comm)
830         
831          (let ((my-results (point-projection label my-comm myrank size target source-tree zone 0 (lambda (x nn) nn))))
832
833          (MPI:barrier my-comm)
834           
835            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
836              (lambda (out-sources)
837                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
838                  (lambda (out-targets)
839                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
840                      (lambda (out-distances)
841                        (for-each
842                         (lambda (my-data)
843                           (let* ((my-entry-len 3)
844                                  (my-data-len (/ (f64vector-length my-data) my-entry-len)))
845                             (d "~A: rank ~A: length my-data = ~A~%" label myrank my-data-len)
846                             (let recur ((m 0))
847                               (if (< m my-data-len)
848                                   (let ((source (inexact->exact (f64vector-ref my-data (* m my-entry-len))))
849                                         (target (inexact->exact (f64vector-ref my-data (+ 1 (* m my-entry-len)))))
850                                         (distance (f64vector-ref my-data (+ 2 (* m my-entry-len)))))
851                                     (fprintf out-sources "~A~%" source)
852                                     (fprintf out-targets "~A~%" target)
853                                     (fprintf out-distances "~A~%" distance)
854                                     (recur (+ 1 m)))))
855                             ))
856                         my-results)
857                        ))
858                    ))
859                ))
860            ))
861       
862
863        (include "calc-parser.scm")
864
865       
866        (define (load-config-file filename)
867          (let ((in (open-input-file filename)))
868            (if (not in) (error 'load-config-file "file not found" filename))
869            (init-bindings)
870            (let* ((lines (reverse (filter (lambda (x) (not (string-null? x))) (read-lines in))))
871                   (properties (filter-map
872                                (lambda (line) 
873                                  (and (not (string-prefix? "//" line))
874                                       (let ((lst (string-split line "=")))
875                                         (and (> (length lst) 1)
876                                              (let ((key (string->symbol (string-trim-both (car lst) #\space)))
877                                                    (val (calc-eval-string (cadr lst))))
878                                                (add-binding key val)
879                                                (cons key val))))))
880                                lines)))
881              properties
882              ))
883          )
884
885        (define (make-output-fname dirname sysname suffix . rest) 
886          (let-optionals 
887           rest ((x #t))
888           (and x
889                (let ((dirname (if (string? x) x dirname)))
890                  (let ((fname (sprintf "~A~A" sysname suffix)))
891                    (or (and dirname (make-pathname dirname fname)) fname)))
892                )))
893
894       
895)
896
897       
Note: See TracBrowser for help on using the repository browser.