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

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

picnic: bug fix in ordering of layouts and local sections

File size: 37.3 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.sorted.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                                        (append (list source target distance segi seci) ax)
436                                        ))
437                                    ax
438                                    (delete-duplicates
439                                     (map (lambda (x) 
440                                            (d "rank ~A: query result = ~A (~A) (~A) ~%" 
441                                               myrank (kdnn-point x) (kdnn-distance x) (kdnn-parent-index x))
442                                            (list (+ fiber-start (kdnn-parent-index x))
443                                                  (+ (kdnn-distance x) (genpoint-parent-distance gd)  (kdnn-parent-distance x))
444                                                  ))
445                                          (kd-tree-near-neighbors* fibers zone (genpoint-coords gd)))
446                                     (lambda (u v)
447                                       (= (car u) (car v)))
448                                     )
449                                    ))
450                                 '() gxs)))
451                           
452                           (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
453
454                                  (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
455                             
456                             (append res1 ax))
457                           
458                           ))
459                       )
460                     ax sections)
461                    ))
462                '() cells ))
463
464
465       
466
467        (define (point-projection prefix my-comm myrank size pts fibers zone point-start nn-filter)
468          (fold (lambda (px ax)
469
470                  (d "~A: rank ~A: px = ~A~%"  prefix myrank px)
471
472                  (let* ((i (+ point-start (car px)))
473                         (root (modulo i size))
474                         (dd (d "~A: rank ~A: querying point ~A (root ~A)~%" prefix myrank px root))
475                         (query-data
476                          (fold 
477                           (lambda (pd ax)
478                             (fold
479                              (lambda (x ax) 
480                                (let ((source (car x))
481                                      (target i)
482                                      (distance (cadr x)))
483                                  (if (and (> distance  0.) (not (= source target)))
484                                      (append (list source target distance) ax)
485                                      ax)
486                                  ))
487                              ax
488                              (delete-duplicates
489                               (map (lambda (x) 
490                                      (let ((res (list (car (cadar x))  (+ (caddr x) (cadr (cadar x))))))
491                                        (d "~A: axon: x = ~A res = ~A~%" prefix x res)
492                                        res))
493                                    (nn-filter pd (kd-tree-near-neighbors* fibers zone pd))
494                                    )
495                               (lambda (u v) (= (car u) (car v)))
496                               )
497                              ))
498                           '() (list (cadr px))))
499                         )
500
501                    (let* ((res0 (MPI:gatherv-f64vector (list->f64vector query-data) root my-comm))
502                           (res1 (or (and (= myrank root) (filter (lambda (x) (not (f64vector-empty? x))) res0)) '())))
503                      (append res1 ax))
504                    ))
505
506                '() pts))
507       
508
509
510
511
512        (define-datatype layer-boundary layer-boundary?
513          (Bounds (top number?) (left number?) (bottom number?) (right number?))
514          (BoundsXZ (top number?) (left number?) (bottom number?) (right number?)
515                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
516          (BoundsYZ (top number?) (left number?) (bottom number?) (right number?)
517                    (n integer?) (k integer?) (x f64vector?) (y f64vector?) (d f64vector?) (d2 f64vector?))
518          )
519
520
521
522        (define (boundary-z-extent-function boundary)
523          (cases layer-boundary boundary
524                 (Bounds (top left bottom right) 
525                         (lambda (x y) 0.))
526                 (BoundsXZ (top left bottom right n k x y d d2) 
527                           (lambda (xp yp) 
528                             (let-values (((y0tab y1tab y2tab res)
529                                           (bvsp-spline:evaluate n k x y d d2 (f64vector xp) 0)))
530                               (f64vector-ref y0tab 0))))
531                 (BoundsYZ (top left bottom right n k x y d d2) 
532                           (lambda (xp yp) 
533                             (let-values (((y0tab y1tab y2tab res)
534                                           (bvsp-spline:evaluate n k x y d d2 (f64vector yp) 0)))
535                               (f64vector-ref y0tab 0))))
536                 ))
537
538
539        (define (point2d-rejection top left bottom right)
540            (lambda (p)
541              (let ((x (coord 0 p)) (y (coord 1 p)))
542                (and (fp> x left)  (fp< x right) (fp> y bottom) (fp< y top) p)))
543            )
544
545
546        (define (compute-point3d p zu z-extent-function)
547          (let* ((x (coord 0 p))
548                 (y (coord 1 p))
549                 (z-extent (z-extent-function x y)))
550            (if (zero? zu)
551                (make-point x y zu)
552                (make-point x y (fp* zu z-extent)))
553            ))
554
555
556        (define (cluster-point-rejection p cluster-pts mean-distance)
557          (let ((D (* 4 mean-distance mean-distance))
558                (nn (kd-tree-nearest-neighbor cluster-pts p)))
559            (and (< (dist2 p nn) D) p)))
560
561
562
563        (define (XZAxis n k x-points z-points boundary)
564         
565          (if (not (>= n 3)) 
566              (error 'generate-boundary "XZAxis requires at least 3 interpolation points (n >= 3)"))
567         
568          (cases layer-boundary boundary
569                 (Bounds (top left bottom right) 
570
571                         (let-values (((d d2 constr errc diagn)
572                                       (bvsp-spline:compute n k x-points z-points)))
573                           
574                           (if (not (zero? errc)) 
575                               (error 'generate-boundary "error in constructing spline from boundary points" errc))
576                           
577                           (BoundsXZ top left bottom right n k x-points z-points d d2)))
578                 
579                 (else (error 'generate-boundary "boundary argument to XZAxis is already a pseudo-3D boundary")))
580          )
581
582
583        (define (Grid x-spacing y-spacing z-spacing boundary)
584
585          (match-let (((top left bottom right)
586                       (cases layer-boundary boundary
587                                   (Bounds (top left bottom right) 
588                                           (list top left bottom right))
589                                   (BoundsXZ (top left bottom right n k x y d d2) 
590                                             (list top left bottom right))
591                                   (BoundsYZ (top left bottom right n k x y d d2) 
592                                             (list top left bottom right))
593                                   )))
594
595          (let* (
596                 (x-extent   (- right left))
597                 (y-extent   (- top bottom))
598                 (z-extent-function
599                  (boundary-z-extent-function boundary))
600                 (compute-grid-points3d
601                  (lambda (p z-spacing z-extent-function)
602                    (let* ((x (coord 0 p))
603                           (y (coord 1 p))
604                           (z-extent (z-extent-function x y)))
605                      (let recur ((points '()) (z (/ z-spacing 2.)))
606                        (if (> z z-extent)
607                            points
608                            (recur (cons (make-point x y z) points) (+ z z-spacing)))
609                        ))
610                    ))
611                 )
612           
613            (d "Grid: x-spacing = ~A~%" x-spacing)
614            (d "Grid: y-spacing = ~A~%" y-spacing)
615            (d "Grid: z-spacing = ~A~%" z-spacing)
616           
617            (d "Grid: x-extent = ~A~%" x-extent)
618            (d "Grid: y-extent = ~A~%" y-extent)
619           
620            (let ((x-points (let recur ((points '()) (x (/ x-spacing 2.)))
621                              (if (> x x-extent)
622                                  (list->f64vector (reverse points))
623                                  (recur (cons x points) (+ x x-spacing)))))
624                  (y-points (let recur ((points '()) (y (/ y-spacing 2.)))
625                              (if (> y y-extent)
626                                  (list->f64vector (reverse points))
627                                  (recur (cons y points) (+ y y-spacing)))))
628                  )
629             
630              (let ((nx (f64vector-length x-points))
631                    (ny (f64vector-length y-points))
632                    )
633               
634                (let recur ((i 0) (j 0) (ax '()))
635                  (if (< i nx)
636                      (let ((x (f64vector-ref x-points i)))
637                        (if (< j ny)
638                            (let* ((y   (f64vector-ref y-points j))
639                                   (p   (make-point x y))
640                                   (p3ds (if (zero? z-spacing)
641                                             (list (make-point x y 0.0))
642                                             (compute-grid-points3d p z-spacing z-extent-function)))
643                                   )
644                              (recur i (+ 1 j) (if p3ds (append p3ds ax) ax)))
645                            (recur (+ 1 i) 0 ax)))
646                      (let* ((t (list->kd-tree ax))
647                             (n (kd-tree-size t)))
648                        (list t boundary)
649                        ))
650                  )))
651            ))
652          )
653
654        (define (UniformRandomPointProcess n x-seed y-seed boundary)
655
656          (match-let (((top left bottom right)
657                       (cases layer-boundary boundary
658                                   (Bounds (top left bottom right) 
659                                           (list top left bottom right))
660                                   (BoundsXZ (top left bottom right n k x y d d2) 
661                                             (list top left bottom right))
662                                   (BoundsYZ (top left bottom right n k x y d d2) 
663                                             (list top left bottom right))
664                                   )))
665
666          (let* (
667                 (x-extent   (- right left))
668                 (y-extent   (- top bottom))
669                 (z-extent-function (boundary-z-extent-function boundary))
670                 )
671
672            (let ((x-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init x-seed)))
673                  (y-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init y-seed)))
674                  (z-points (random-mtzig:f64vector-randu! (inexact->exact n) (random-mtzig:init (current-milliseconds)))))
675             
676              (let ((point-rejection1 (point2d-rejection top left bottom right)))
677               
678                (let recur ((i 0) (ax '()))
679                  (if (< i n)
680                      (let ((x (f64vector-ref x-points i))
681                            (y (f64vector-ref y-points i))
682                            (z (f64vector-ref z-points i)))
683                        (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
684                          (let ((p3d 
685                                 (and (point-rejection1 p)
686                                      (compute-point3d p z z-extent-function))))
687                            (recur (+ 1 i) (if p3d (cons p3d ax) ax)))))
688                      (let* ((t (list->kd-tree ax))
689                             (n (kd-tree-size t)))
690
691                        (list t boundary))))
692                ))
693            ))
694          )
695
696
697        (define (ClusteredRandomPointProcess cluster-pts n mean-distance x-seed y-seed boundary)
698
699          (match-let (((top left bottom right)
700                       (cases layer-boundary boundary
701                                   (Bounds (top left bottom right) 
702                                           (list top left bottom right))
703                                   (BoundsXZ (top left bottom right n k x y d d2) 
704                                             (list top left bottom right))
705                                   (BoundsYZ (top left bottom right n k x y d d2) 
706                                             (list top left bottom right))
707                                   )))
708
709
710          (let* (
711                 (x-extent   (- right left))
712                 (y-extent   (- top bottom))
713                 (z-extent-function (boundary-z-extent-function boundary))
714                 )
715
716            (let recur ((pts '()) (x-seed x-seed) (y-seed y-seed))
717
718              (let ((x-points (random-mtzig:f64vector-randu! n (random-mtzig:init x-seed)))
719                    (y-points (random-mtzig:f64vector-randu! n (random-mtzig:init y-seed)))
720                    (z-points (random-mtzig:f64vector-randu! n (random-mtzig:init (current-milliseconds)))))
721               
722                (let ((point-rejection1 (point2d-rejection top left bottom right)))
723                 
724                  (let inner-recur ((j 0) (ax pts))
725                    (if (< j n)
726                        (let ((x (f64vector-ref x-points j))
727                              (y (f64vector-ref y-points j))
728                              (z (f64vector-ref z-points j)))
729                          (let ((p (make-point (fp* x x-extent) (fp* y y-extent))))
730                            (let ((p3d 
731                                   (and (point-rejection1 p)
732                                        (compute-point3d p z z-extent-function))))
733                              (let ((pp (cluster-point-rejection p3d cluster-pts mean-distance)))
734                                (inner-recur (+ 1 j)  (if pp (cons pp ax) ax))))))
735
736                        (if (< (length ax) n)
737                            (recur ax (+ 1 x-seed) (+ 1 y-seed))
738                            (let* ((t (list->kd-tree (take ax n)))
739                                   (n (kd-tree-size t)))
740                             
741                              (list t boundary))))
742                    ))
743                ))
744            ))
745          )
746
747
748
749       
750        (define comment-pat (string->irregex "^#.*"))
751
752
753        (define (load-points-from-file filename)
754
755          (let ((in (open-input-file filename)))
756           
757            (if (not in) (error 'load-points-from-file "file not found" filename))
758
759            (let* ((lines
760                    (filter (lambda (line) (not (irregex-match comment-pat line)))
761                            (read-lines in)))
762
763                   (point-data
764                    (filter-map
765                     (lambda (line) 
766                       (let ((lst (map string->number (string-split line " \t"))))
767                         (and (not (null? lst)) (apply make-point lst))))
768                     lines))
769
770                   (point-tree (list->kd-tree point-data))
771                   )
772             
773              (list point-tree)
774             
775              ))
776          )
777
778
779        (define (segment-projection label source-tree target-sections zone my-comm myrank size)
780
781          (MPI:barrier my-comm)
782         
783          (let ((my-results
784                 (genpoint-projection label my-comm myrank size target-sections source-tree zone 0 0)))
785
786            (MPI:barrier my-comm)
787
788            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
789              (lambda (out-sources)
790                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
791                  (lambda (out-targets)
792                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
793                      (lambda (out-distances)
794                        (call-with-output-file (sprintf "~Asegments~A.dat"  label (if (> size 1) myrank ""))
795                          (lambda (out-segments)
796                            (for-each
797                             (lambda (my-data)
798                               (let* ((my-entry-len 5)
799                                      (my-data-len (/ (f64vector-length my-data) my-entry-len)))
800                                 (d "rank ~A: length my-data = ~A~%" myrank my-data-len)
801                                 (let recur ((m 0))
802                                   (if (< m my-data-len)
803                                       (let* (
804                                              (my-entry-offset (* m my-entry-len))
805                                              (source (inexact->exact (f64vector-ref my-data my-entry-offset)))
806                                              (target (inexact->exact (f64vector-ref my-data (+ 1 my-entry-offset))))
807                                              (distance (f64vector-ref my-data (+ 2 my-entry-offset)))
808                                              (segment (f64vector-ref my-data (+ 3 my-entry-offset)))
809                                              (section (f64vector-ref my-data (+ 4 my-entry-offset)))
810                                              )
811                                         (fprintf out-sources   "~A~%" source)
812                                         (fprintf out-targets   "~A~%" target)
813                                         (fprintf out-distances "~A~%" distance)
814                                         (fprintf out-segments  "~A ~A~%" segment section)
815                                         (recur (+ 1 m)))))
816                                 ))
817                             my-results))
818                          ))
819                      ))
820                  ))
821              ))
822          )
823       
824
825        (define (projection label source-tree target zone my-comm myrank size) 
826
827          (MPI:barrier my-comm)
828         
829          (let ((my-results (point-projection label my-comm myrank size target source-tree zone 0 (lambda (x nn) nn))))
830
831          (MPI:barrier my-comm)
832           
833            (call-with-output-file (sprintf "~Asources~A.dat"  label (if (> size 1) myrank ""))
834              (lambda (out-sources)
835                (call-with-output-file (sprintf "~Atargets~A.dat"  label (if (> size 1) myrank ""))
836                  (lambda (out-targets)
837                    (call-with-output-file (sprintf "~Adistances~A.dat"  label (if (> size 1) myrank ""))
838                      (lambda (out-distances)
839                        (for-each
840                         (lambda (my-data)
841                           (let* ((my-entry-len 3)
842                                  (my-data-len (/ (f64vector-length my-data) my-entry-len)))
843                             (d "~A: rank ~A: length my-data = ~A~%" label myrank my-data-len)
844                             (let recur ((m 0))
845                               (if (< m my-data-len)
846                                   (let ((source (inexact->exact (f64vector-ref my-data (* m my-entry-len))))
847                                         (target (inexact->exact (f64vector-ref my-data (+ 1 (* m my-entry-len)))))
848                                         (distance (f64vector-ref my-data (+ 2 (* m my-entry-len)))))
849                                     (fprintf out-sources "~A~%" source)
850                                     (fprintf out-targets "~A~%" target)
851                                     (fprintf out-distances "~A~%" distance)
852                                     (recur (+ 1 m)))))
853                             ))
854                         my-results)
855                        ))
856                    ))
857                ))
858            ))
859       
860
861        (include "calc-parser.scm")
862
863       
864        (define (load-config-file filename)
865          (let ((in (open-input-file filename)))
866            (if (not in) (error 'load-config-file "file not found" filename))
867            (init-bindings)
868            (let* ((lines (reverse (filter (lambda (x) (not (string-null? x))) (read-lines in))))
869                   (properties (filter-map
870                                (lambda (line) 
871                                  (and (not (string-prefix? "//" line))
872                                       (let ((lst (string-split line "=")))
873                                         (and (> (length lst) 1)
874                                              (let ((key (string->symbol (string-trim-both (car lst) #\space)))
875                                                    (val (calc-eval-string (cadr lst))))
876                                                (add-binding key val)
877                                                (cons key val))))))
878                                lines)))
879              properties
880              ))
881          )
882
883        (define (make-output-fname dirname sysname suffix . rest) 
884          (let-optionals 
885           rest ((x #t))
886           (and x
887                (let ((dirname (if (string? x) x dirname)))
888                  (let ((fname (sprintf "~A~A" sysname suffix)))
889                    (or (and dirname (make-pathname dirname fname)) fname)))
890                )))
891
892       
893)
894
895       
Note: See TracBrowser for help on using the repository browser.