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

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

picnic: bug fixes in write-sections and added some plotting scripts

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